Defining Penalized Spline Smooths in VGAM Formulas
sm.ps.RdThis function represents a P-spline smooth term
in a vgam formula
and confers automatic smoothing parameter selection.
Usage
sm.ps(x, ..., ps.int = NULL, spar = -1, degree = 3, p.order = 2,
ridge.adj = 1e-5, spillover = 0.01, maxspar = 1e12,
outer.ok = FALSE, mux = NULL, fixspar = FALSE)Arguments
- x, ...
See
sm.os.- ps.int
the number of equally-spaced B-spline intervals. Note that the number of knots is equal to
ps.int + 2*degree + 1. The default, signified byNULL, means that the maximum of the value 7 anddegreeis chosen. This usually means 6 interior knots for big data sets. However, if this is too high compared to the length ofx, then some adjustment is made. In the case wheremuxis assigned a numerical value (suggestions: some value between 1 and 2) thenceiling(mux * log(length(unique(x.index))))is used, wherex.indexis the combined data. No matter what, the above is not guaranteed to work on every data set. This argument may change in the future. See also argumentmux.- spar, maxspar
See
sm.os.- mux
numeric. If given, then this argument multiplies
log(length(unique(x)))to obtainps.int. Ifps.intis given then this argument is ignored.- degree
degree of B-spline basis. Usually this will be 2 or 3; and the values 1 or 4 might possibly be used.
- p.order
order of difference penalty (0 is the ridge penalty).
- ridge.adj, spillover
See
sm.os.- outer.ok, fixspar
See
sm.os.
Details
This function can be used by vgam to
allow automatic smoothing parameter selection based on
P-splines and minimizing an UBRE quantity.
This function should only be used with vgam
and is an alternative to sm.os;
see that function for some details that also apply here.
Value
A matrix with attributes that are (only) used by vgam.
The number of rows of the matrix is length(x) and
the number of columns is ps.int + degree - 1.
The latter is because the function is centred.
References
Eilers, P. H. C. and Marx, B. D. (1996). Flexible smoothing with B-splines and penalties (with comments and rejoinder). Statistical Science, 11(2): 89–121.
Author
B. D. Marx wrote the original function. Subsequent edits were made by T. W. Yee and C. Somchit.
Note
This function is currently under development and
may change in the future.
In particular, the default for ps.int is
subject to change.
Warning
See sm.os.
See also
sm.os,
s,
vgam,
smartpred,
is.smart,
summarypvgam,
splineDesign,
bs,
magic.
Examples
sm.ps(runif(20))
#> 2 3 4 5 6 7
#> 1 0.56848499 0.04407027 -0.10896053 -0.08828988 -0.15682395 -0.16443919
#> 2 0.13552785 0.50570134 0.01635008 -0.08352490 -0.14836022 -0.15556446
#> 3 -0.07129744 0.01244535 0.57878525 0.13732230 -0.11935405 -0.12521252
#> 4 -0.09200077 -0.15632893 -0.10705466 0.11755306 0.50776414 -0.02778718
#> 5 -0.10371988 -0.17624219 -0.12075757 -0.09780058 -0.17319275 0.06694847
#> 6 0.35852401 0.37036452 -0.06949543 -0.07958894 -0.14136900 -0.14823376
#> 7 0.03445295 0.50604843 0.11321354 -0.08434429 -0.15005959 -0.15734636
#> 8 -0.06086662 -0.10342536 -0.07086496 -0.05739296 -0.10194362 -0.10684005
#> 9 -0.06382017 -0.10533492 0.25081240 0.54547728 -0.04077088 -0.11208094
#> 10 0.24035650 0.45766835 -0.03521869 -0.08151531 -0.14479070 -0.15182161
#> 11 -0.09675126 -0.16440103 -0.11264424 0.01683433 0.48775226 0.07181819
#> 12 -0.10610306 -0.18029172 -0.12353222 -0.10004774 -0.17010122 0.20005699
#> 13 -0.07588585 0.12293330 0.55609079 0.02927170 -0.12809255 -0.13431262
#> 14 -0.09806671 -0.16663625 -0.11417577 -0.09247003 -0.16424891 -0.05042176
#> 15 -0.10218683 -0.17363721 -0.11897269 -0.07852650 0.29205472 0.31534863
#> 16 -0.10652897 -0.18101544 -0.12402810 -0.10044935 -0.16456955 0.25145373
#> 17 -0.10664385 -0.18121064 -0.12416185 -0.10055768 -0.12457763 0.39673467
#> 18 -0.06986577 -0.11871684 -0.08134237 -0.06587852 -0.11701603 -0.12031949
#> 19 -0.10123982 -0.17202803 -0.11787011 -0.06506398 0.35026873 0.25845211
#> 20 -0.08236929 -0.13996298 -0.08617287 0.32899199 0.40743080 -0.10643285
#> 8 9 10
#> 1 -0.17310713 -0.088374733 -0.010897812
#> 2 -0.16376460 -0.083605177 -0.010309661
#> 3 -0.13181275 -0.067293103 -0.008298159
#> 4 -0.17008849 -0.086833655 -0.010707776
#> 5 0.45501121 0.005714018 -0.012071738
#> 6 -0.15604748 -0.079665430 -0.009823836
#> 7 -0.16564042 -0.084562824 -0.010427752
#> 8 0.09064512 0.604670433 0.127569798
#> 9 -0.11798897 -0.060235783 -0.007427895
#> 10 -0.15982446 -0.081593654 -0.010061612
#> 11 -0.17846573 -0.091317340 -0.011260676
#> 12 0.36560771 -0.055914637 -0.012349111
#> 13 -0.14139253 -0.072183776 -0.008901246
#> 14 0.47600142 0.128170758 -0.011250820
#> 15 -0.16476208 -0.096447623 -0.011893310
#> 16 0.32082375 -0.070709855 -0.012398682
#> 17 0.15966892 -0.095542860 -0.012412053
#> 18 0.18106762 0.549022148 0.064292217
#> 19 -0.17364905 -0.095553797 -0.011783089
#> 20 -0.15228207 -0.077743110 -0.009586788
#> attr(,"S.arg")
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 0.95808349 -0.527757676 0.273760538 0.10928301 0.13561954
#> [2,] -0.52775768 0.757395026 -0.661944773 0.15030832 -0.10994632
#> [3,] 0.27376054 -0.661944773 0.971426480 -0.58109289 0.16613784
#> [4,] 0.10928301 0.150308324 -0.581092886 0.98636681 -0.59335587
#> [5,] 0.13561954 -0.109946319 0.166137835 -0.59335587 0.89573902
#> [6,] 0.18006973 -0.050945154 0.054452483 0.22513442 -0.60537057
#> [7,] 0.15360546 -0.114727725 0.015460272 0.03861135 0.11669232
#> [8,] 0.09915802 -0.023330379 0.032038839 0.03926758 0.01454059
#> [9,] 0.01017254 -0.006368835 0.001558261 0.00290451 -0.00164880
#> [,6] [,7] [,8] [,9]
#> [1,] 0.180069733 0.15360546 0.09915802 0.010172540
#> [2,] -0.050945154 -0.11472772 -0.02333038 -0.006368835
#> [3,] 0.054452483 0.01546027 0.03203884 0.001558261
#> [4,] 0.225134421 0.03861135 0.03926758 0.002904510
#> [5,] -0.605370568 0.11669232 0.01454059 -0.001648800
#> [6,] 1.024580443 -0.59647550 0.20723466 0.002678119
#> [7,] -0.596475502 0.90105337 -0.60526453 0.154884429
#> [8,] 0.207234662 -0.60526453 0.81089985 -0.310783350
#> [9,] 0.002678119 0.15488443 -0.31078335 0.156222509
#> attr(,"degree")
#> [1] 3
#> attr(,"knots")
#> [1] -0.38118833 -0.24557179 -0.10995526 0.02566128 0.16127782 0.29689435
#> [7] 0.43251089 0.56812742 0.70374396 0.83936050 0.97497703 1.11059357
#> [13] 1.24621010 1.38182664
#> attr(,"spar")
#> [1] -1
#> attr(,"p.order")
#> [1] 2
#> attr(,"ps.int")
#> [1] 7
#> attr(,"ridge.adj")
#> [1] 1e-05
#> attr(,"outer.ok")
#> [1] FALSE
#> attr(,"fixspar")
#> [1] FALSE
sm.ps(runif(20), ps.int = 5)
#> 2 3 4 5 6 7
#> 1 -0.1100373884 -0.16794882 -0.19017801 0.41620128 0.22690429 -0.08672579
#> 2 0.5518215500 0.02060989 -0.25396136 -0.17992137 -0.12101569 -0.09307072
#> 3 0.3931534975 0.31029994 -0.20952529 -0.16117335 -0.10840571 -0.08337263
#> 4 -0.1410237149 -0.20292385 0.10923407 0.30196244 -0.11912309 -0.11666712
#> 5 -0.0937768086 -0.14313048 -0.20815080 0.18320698 0.49950548 -0.01470148
#> 6 -0.0586481705 -0.08951404 -0.13240334 -0.09377554 0.12923210 0.61580376
#> 7 -0.0005573995 0.47001426 -0.05733899 -0.19632390 -0.13215370 -0.10163673
#> 8 -0.1299228362 0.21928444 0.19288862 -0.20021723 -0.15469498 -0.11897278
#> 9 -0.1524297858 -0.13709894 0.29631877 0.01941787 -0.16315819 -0.12610322
#> 10 -0.0602772427 -0.09200047 -0.13608111 -0.09624572 0.15490790 0.60774459
#> 11 -0.0696532735 -0.10631100 -0.15724832 -0.09910999 0.35243235 0.47009571
#> 12 -0.1515106867 -0.15438039 0.27879503 0.05801963 -0.16101843 -0.12534286
#> 13 0.3938167877 0.30963264 -0.20962653 -0.16115936 -0.10839630 -0.08336540
#> 14 -0.1198920772 -0.18298992 -0.13437626 0.47085345 0.07210279 -0.09913969
#> 15 -0.1002601865 -0.15302599 -0.21045749 0.29140895 0.39779914 -0.05623367
#> 16 -0.1533530974 -0.09605299 0.31685617 -0.04636537 -0.16492182 -0.12686706
#> 17 0.4579325463 0.23649697 -0.21948890 -0.16114509 -0.10838670 -0.08335802
#> 18 -0.1489834824 0.06118546 0.28613942 -0.15961153 -0.16199238 -0.12458506
#> 19 -0.1533538842 -0.06388498 0.32040836 -0.08205243 -0.16496047 -0.12686776
#> 20 -0.1530443470 -0.03826170 0.31819598 -0.10396973 -0.16465661 -0.12663407
#> 8
#> 1 -0.01397921
#> 2 -0.01429220
#> 3 -0.01280293
#> 4 -0.01791573
#> 5 -0.01191345
#> 6 0.13588800
#> 7 -0.01560762
#> 8 -0.01826979
#> 9 -0.01936476
#> 10 0.11482915
#> 11 0.02378944
#> 12 -0.01924800
#> 13 -0.01280182
#> 14 -0.01523115
#> 15 -0.01273711
#> 16 -0.01948206
#> 17 -0.01280069
#> 18 -0.01913163
#> 19 -0.01948217
#> 20 -0.01944628
#> attr(,"S.arg")
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 1.25279309 -0.68300904 0.41441780 0.20266289 0.13531510 0.07495159
#> [2,] -0.68300904 1.21704547 -0.90690330 0.24271331 0.02160297 -0.02782579
#> [3,] 0.41441780 -0.90690330 1.10445943 -0.80849381 0.22279055 -0.05461430
#> [4,] 0.20266289 0.24271331 -0.80849381 1.35823060 -0.76213095 0.21652803
#> [5,] 0.13531510 0.02160297 0.22279055 -0.76213095 1.29681892 -0.82864603
#> [6,] 0.07495159 -0.02782579 -0.05461430 0.21652803 -0.82864603 1.02118385
#> [7,] 0.02149145 0.01096192 0.01414778 0.01722195 0.21979027 -0.41155435
#> [,7]
#> [1,] 0.02149145
#> [2,] 0.01096192
#> [3,] 0.01414778
#> [4,] 0.01722195
#> [5,] 0.21979027
#> [6,] -0.41155435
#> [7,] 0.21038647
#> attr(,"degree")
#> [1] 3
#> attr(,"knots")
#> [1] -0.41829825 -0.24520123 -0.07210422 0.10099280 0.27408981 0.44718683
#> [7] 0.62028384 0.79338086 0.96647788 1.13957489 1.31267191 1.48576892
#> attr(,"spar")
#> [1] -1
#> attr(,"p.order")
#> [1] 2
#> attr(,"ps.int")
#> [1] 5
#> attr(,"ridge.adj")
#> [1] 1e-05
#> attr(,"outer.ok")
#> [1] FALSE
#> attr(,"fixspar")
#> [1] FALSE
if (FALSE) { # \dontrun{
data("TravelMode", package = "AER") # Need to install "AER" first
air.df <- subset(TravelMode, mode == "air") # Form 4 smaller data frames
bus.df <- subset(TravelMode, mode == "bus")
trn.df <- subset(TravelMode, mode == "train")
car.df <- subset(TravelMode, mode == "car")
TravelMode2 <- data.frame(income = air.df$income,
wait.air = air.df$wait - car.df$wait,
wait.trn = trn.df$wait - car.df$wait,
wait.bus = bus.df$wait - car.df$wait,
gcost.air = air.df$gcost - car.df$gcost,
gcost.trn = trn.df$gcost - car.df$gcost,
gcost.bus = bus.df$gcost - car.df$gcost,
wait = air.df$wait) # Value is unimportant
TravelMode2$mode <- subset(TravelMode, choice == "yes")$mode # The response
TravelMode2 <- transform(TravelMode2, incom.air = income, incom.trn = 0,
incom.bus = 0)
set.seed(1)
TravelMode2 <- transform(TravelMode2,
junkx2 = runif(nrow(TravelMode2)))
tfit2 <-
vgam(mode ~ sm.ps(gcost.air, gcost.trn, gcost.bus) + ns(junkx2, 4) +
sm.ps(incom.air, incom.trn, incom.bus) + wait ,
crit = "coef",
multinomial(parallel = FALSE ~ 1), data = TravelMode2,
xij = list(sm.ps(gcost.air, gcost.trn, gcost.bus) ~
sm.ps(gcost.air, gcost.trn, gcost.bus) +
sm.ps(gcost.trn, gcost.bus, gcost.air) +
sm.ps(gcost.bus, gcost.air, gcost.trn),
sm.ps(incom.air, incom.trn, incom.bus) ~
sm.ps(incom.air, incom.trn, incom.bus) +
sm.ps(incom.trn, incom.bus, incom.air) +
sm.ps(incom.bus, incom.air, incom.trn),
wait ~ wait.air + wait.trn + wait.bus),
form2 = ~ sm.ps(gcost.air, gcost.trn, gcost.bus) +
sm.ps(gcost.trn, gcost.bus, gcost.air) +
sm.ps(gcost.bus, gcost.air, gcost.trn) +
wait +
sm.ps(incom.air, incom.trn, incom.bus) +
sm.ps(incom.trn, incom.bus, incom.air) +
sm.ps(incom.bus, incom.air, incom.trn) +
junkx2 + ns(junkx2, 4) +
incom.air + incom.trn + incom.bus +
gcost.air + gcost.trn + gcost.bus +
wait.air + wait.trn + wait.bus)
par(mfrow = c(2, 2))
plot(tfit2, se = TRUE, lcol = "orange", scol = "blue", ylim = c(-4, 4))
summary(tfit2)
} # }