sm.ps: Defining Penalized Spline Smooths in VGAM Formulas

Description Usage Arguments Details Value Warning Note Author(s) References See Also Examples

View source: R/sm.ps.R

Description

This function represents a P-spline smooth term in a vgam formula and confers automatic smoothing parameter selection.

Usage

1
2
3
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 by NULL, means that the maximum of the value 7 and degree is chosen. This usually means 6 interior knots for big data sets. However, if this is too high compared to the length of x, then some adjustment is made. In the case where mux is assigned a numerical value (suggestions: some value between 1 and 2) then ceiling(mux * log(length(unique(x.index)))) is used, where x.index is 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 argument mux.

spar, maxspar

See sm.os.

mux

numeric. If given, then this argument multiplies log(length(unique(x))) to obtain ps.int. If ps.int is 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.

Warning

See sm.os.

Note

This function is currently under development and may change in the future. In particular, the default for ps.int is subject to change.

Author(s)

B. D. Marx wrote the original function. Subsequent edits were made by T. W. Yee and C. Somchit.

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.

See Also

sm.os, s, vgam, smartpred, is.smart, summarypvgam, splineDesign, bs, magic.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
sm.ps(runif(20))
sm.ps(runif(20), ps.int = 5)

## Not run: 
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)

## End(Not run)

Example output

Loading required package: stats4
Loading required package: splines
             2           3           4            5           6           7
1  -0.04148624 -0.09919985 -0.15698790 -0.209226808 -0.08270623  0.46065989
2  -0.04677059 -0.11183555 -0.17698441 -0.090994535  0.46466482 -0.01135228
3  -0.04952448 -0.11842049 -0.16000742  0.258650703  0.23703158 -0.19829257
4  -0.05016527 -0.11995274 -0.11549558  0.364559351  0.09143631 -0.21427920
5  -0.02884649 -0.06897631 -0.10915789 -0.145480999 -0.12333361  0.03027080
6  -0.01291211  0.39628562  0.35867070 -0.144408937 -0.14079572 -0.14209974
7  -0.03422122 -0.08182809 -0.12949636 -0.172587260 -0.14186433  0.19974177
8   0.07654279  0.58114621  0.14374827 -0.138697744 -0.11801439 -0.11910742
9  -0.04236036 -0.10129001 -0.16029565 -0.213635242 -0.05350259  0.47693342
10 -0.04754883 -0.10322134  0.23314659  0.300211437 -0.16686030 -0.20517862
11 -0.01415694 -0.03385137 -0.05357121 -0.071397431 -0.06052820 -0.06103493
12 -0.04838644 -0.11102400  0.16762555  0.344344073 -0.15064719 -0.20879300
13 -0.04082895  0.04249811  0.50909872 -0.009667189 -0.17453527 -0.17618161
14 -0.04435026 -0.10604816 -0.16782562 -0.222552985  0.08577021  0.44307808
15 -0.02139974  0.34822230  0.39735013 -0.139201675 -0.14522720 -0.14657227
16 -0.04978560 -0.11894338  0.02394222  0.408703677 -0.08508277 -0.21483055
17 -0.04287043 -0.10250966 -0.16222580 -0.216207666 -0.03040804  0.48088300
18 -0.04739260 -0.11332285 -0.17913876 -0.014234208  0.45336829 -0.08547994
19 -0.04942130 -0.11817378 -0.16317930  0.244126483  0.25339644 -0.19516419
20  0.63588503  0.14044534 -0.09921628 -0.132303045 -0.11216181 -0.11320064
             8           9           10
1   0.15016136 -0.04587758 -0.006993037
2  -0.12923393 -0.05267456 -0.007883783
3  -0.13686006 -0.05577607 -0.008347985
4  -0.13863089 -0.05649776 -0.008456000
5   0.58636024  0.14668598 -0.004860015
6  -0.09100344 -0.03708762 -0.005550892
7   0.49612517  0.01890479 -0.005768425
8  -0.07627871 -0.03108669 -0.004652735
9   0.09550197 -0.04760467 -0.007140381
10 -0.13140038 -0.05355103 -0.008014964
11  0.16405119  0.64617455  0.132267602
12 -0.13371510 -0.05449437 -0.008156154
13 -0.11283013 -0.04598289 -0.006882244
14 -0.03352401 -0.04994870 -0.007475804
15 -0.09386773 -0.03825494 -0.005725604
16 -0.13758166 -0.05607015 -0.008392001
17  0.06276558 -0.04827826 -0.007226360
18 -0.13096865 -0.05337508 -0.007988630
19 -0.13657493 -0.05565987 -0.008330593
20 -0.07249589 -0.02954504 -0.004421996
attr(,"S.arg")
              [,1]         [,2]         [,3]         [,4]         [,5]
 [1,]  0.850856881 -0.556109582  0.294596217  0.172624137  0.161712434
 [2,] -0.556109582  0.868970283 -0.663332392  0.077047462 -0.030398275
 [3,]  0.294596217 -0.663332392  0.987803217 -0.602450772  0.233519799
 [4,]  0.172624137  0.077047462 -0.602450772  0.908254585 -0.572288882
 [5,]  0.161712434 -0.030398275  0.233519799 -0.572288882  1.047892029
 [6,]  0.147934062 -0.067207343  0.020179133  0.132407500 -0.578898850
 [7,]  0.100129376 -0.030153726  0.033317584  0.011911689  0.208817065
 [8,]  0.042707744 -0.007743472  0.020771571  0.014441416  0.029550674
 [9,]  0.004902255 -0.004721296 -0.002528666 -0.005352033 -0.001946817
              [,6]        [,7]         [,8]         [,9]
 [1,]  0.147934062  0.10012938  0.042707744  0.004902255
 [2,] -0.067207343 -0.03015373 -0.007743472 -0.004721296
 [3,]  0.020179133  0.03331758  0.020771571 -0.002528666
 [4,]  0.132407500  0.01191169  0.014441416 -0.005352033
 [5,] -0.578898850  0.20881706  0.029550674 -0.001946817
 [6,]  0.918109898 -0.61416139  0.168869898 -0.004539835
 [7,] -0.614161388  0.95933512 -0.610848115  0.154251075
 [8,]  0.168869898 -0.61084811  0.789158358 -0.312994218
 [9,] -0.004539835  0.15425107 -0.312994218  0.155924906
attr(,"degree")
[1] 3
attr(,"knots")
 [1] -0.36316826 -0.23360946 -0.10405067  0.02550813  0.15506693  0.28462573
 [7]  0.41418453  0.54374333  0.67330213  0.80286092  0.93241972  1.06197852
[13]  1.19153732  1.32109612
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
             2           3           4            5           6           7
1  -0.09658255 -0.24775600 -0.12096211  0.319644390  0.01632236 -0.05699600
2  -0.09944607 -0.24370003  0.12764541  0.176333987 -0.12365017 -0.05868589
3  -0.08901479  0.01927810  0.38307413 -0.217978092 -0.14230780 -0.05283816
4  -0.09342340 -0.23965207 -0.16469050  0.315698287  0.08932404 -0.05478223
5  -0.04730755 -0.12135451 -0.13935533 -0.067943194  0.57041302  0.22389266
6   0.19863754  0.46694627 -0.10367994 -0.241252500 -0.10673073 -0.03962857
7  -0.09564605 -0.17320189  0.33283906 -0.032763653 -0.14961092 -0.05644339
8   0.60347031  0.03621909 -0.17923501 -0.218618371 -0.09671734 -0.03591065
9  -0.01150295  0.40480762  0.08735249 -0.282069680 -0.12606912 -0.04680882
10 -0.07639476 -0.19596978 -0.20916989  0.177173903  0.38434259 -0.01834405
11 -0.08816545 -0.22616424 -0.19712391  0.283144825  0.19379390 -0.04842724
12 -0.02181221 -0.05595324 -0.06425293 -0.078343030  0.15765119  0.65145065
13 -0.05844536  0.28442607  0.21356756 -0.282617435 -0.13262207 -0.04924190
14 -0.09584281 -0.24585841 -0.13440171  0.320851345  0.03455753 -0.05654963
15  0.28056080  0.43072662 -0.12774776 -0.227467617 -0.10063227 -0.03736424
16 -0.04925242  0.31813470  0.18209171 -0.284348868 -0.13109707 -0.04867567
17 -0.09775070 -0.25072377 -0.09196990  0.312471768 -0.01502035 -0.05768540
18 -0.09911792 -0.25349746 -0.03081913  0.285307177 -0.06085249 -0.05849224
19  0.13338383  0.48106687 -0.07361699 -0.253326550 -0.11207232 -0.04161188
20 -0.09634749 -0.18777395  0.31045475 -0.003896692 -0.14902196 -0.05685733
              8
1  -0.008662997
2  -0.008919841
3  -0.008031027
4  -0.008379636
5  -0.003670094
6  -0.006023264
7  -0.008578997
8  -0.005458166
9  -0.007114610
10 -0.006852248
11 -0.007908023
12  0.141382240
13 -0.007484421
14 -0.008596646
15 -0.005679102
16 -0.007398359
17 -0.008767774
18 -0.008890408
19 -0.006324712
20 -0.008641913
attr(,"S.arg")
             [,1]        [,2]         [,3]        [,4]        [,5]        [,6]
[1,]  1.228215441 -0.63881992  0.525388625  0.31060151  0.18259572  0.05788106
[2,] -0.638819916  1.02037952 -0.856680086 -0.01529768  0.01697425 -0.01913379
[3,]  0.525388625 -0.85668009  1.499173421 -0.75354784  0.37673294  0.03331664
[4,]  0.310601512 -0.01529768 -0.753547836  1.07398770 -0.74887051  0.20407042
[5,]  0.182595717  0.01697425  0.376732941 -0.74887051  1.35918228 -0.80855445
[6,]  0.057881061 -0.01913379  0.033316643  0.20407042 -0.80855445  1.04501535
[7,]  0.005613929 -0.01107481 -0.004314089 -0.01208529  0.20703963 -0.41803641
             [,7]
[1,]  0.005613929
[2,] -0.011074812
[3,] -0.004314089
[4,] -0.012085291
[5,]  0.207039631
[6,] -0.418036411
[7,]  0.207839589
attr(,"degree")
[1] 3
attr(,"knots")
 [1] -0.54555108 -0.35180864 -0.15806620  0.03567624  0.22941868  0.42316112
 [7]  0.61690356  0.81064600  1.00438845  1.19813089  1.39187333  1.58561577
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
Warning message:
In vgam.fit(x = x, y = y, w = w, mf = mf, Xm2 = Xm2, Ym2 = Ym2,  :
  convergence not obtained in 10 outer iterations

Call:
vgam(formula = mode ~ sm.ps(gcost.air, gcost.trn, gcost.bus) + 
    ns(junkx2, 4) + sm.ps(incom.air, incom.trn, incom.bus) + 
    wait, family = multinomial(parallel = FALSE ~ 1), data = TravelMode2, 
    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, crit = "coef", 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))

Pearson residuals:
                       Min       1Q   Median        3Q     Max
log(mu[,1]/mu[,4]) -1.2621 -0.56407 -0.11802  0.200377 10.7175
log(mu[,2]/mu[,4]) -1.1754 -0.49852 -0.17142  0.324806 26.3111
log(mu[,3]/mu[,4]) -1.0092 -0.29804 -0.11527 -0.030093  6.1641

Parametric coefficients: 
                Estimate Std. Error z value  Pr(>|z|)    
(Intercept):1   4.923398   1.116642  4.4091 1.038e-05 ***
(Intercept):2   3.980882   0.911226  4.3687 1.250e-05 ***
(Intercept):3   3.209328   0.922100  3.4805 0.0005006 ***
ns(junkx2, 4)1  0.020934   0.895317  0.0234 0.9813461    
ns(junkx2, 4)2 -0.258297   0.906580 -0.2849 0.7757107    
ns(junkx2, 4)3 -0.125030   1.976111 -0.0633 0.9495508    
ns(junkx2, 4)4  0.617522   0.851002  0.7256 0.4680584    
wait           -0.097017   0.010465 -9.2704 < 2.2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                                          edf Est.rank  Chi.sq p-value  
sm.ps(gcost.air, gcost.trn, gcost.bus) 1.9205        4 13.1574 0.01053 *
sm.ps(incom.air, incom.trn, incom.bus) 3.9428        8  3.3788 0.90839  
---
Signif. codes:  0***0.001**0.01*0.05.’ 0.1 ‘ ’ 1

Number of linear/additive predictors:    3 

Names of linear/additive predictors: log(mu[,1]/mu[,4]), log(mu[,2]/mu[,4]), log(mu[,3]/mu[,4]) 

Dispersion Parameter for multinomial family:   1

Residual deviance:  391.6181 on 640 degrees of freedom

Log-likelihood: -195.8091 on 640 degrees of freedom

Number of outer iterations:  10 

Number of IRLS iterations at final outer iteration:  3 

VGAM documentation built on Jan. 16, 2021, 5:21 p.m.