sm.ps {VGAM}

Defining Penalized Spline Smooths in VGAM Formulas

Description

This 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 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

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)