Instructions and Examples for P-spline S-plus Code: ps( )

The ps() function can be used within the gam() function, just as the commonly used smoothing splines s( ) and loess lo( ) terms. However the arguments for ps( ) differ:

“ps” <- function(x, ps.intervals=NULL, lambda=0, degree=3, order=3, ridge.adj=0.0001, ridge.inv=0.0001)

· x: regressor of interest

· ps.intervals: number of equally-spaced B-spline intervals (knots=ps.int+2*degree+1)

· lambda: non-negative regularization parameter for difference penalty

· degree: degree of B-splines

· order: order of difference penalty (natural number, 0 is ridge penalty)

· ridge.adj, ridge.inv: small positive numbers to stabilize linear dependencies among B-splines bases.

References:

Eilers, P.H.C. and Marx, B.D. (2002). Generalized linear additive smooth structures (GLASS). Journal of Computational and Graphical Statistics 11(4): 758-783.

Marx, B.D. and Eilers, P.H.C. (1998). Direct generalized additive modeling with penalized likelihood. Computational Statistics and Data Analysis 28(2): 193-209.

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.

Some comments:

The spirit of P-splines is to purposely overfit using a modest number (say 20) of equally spaced B-spline bases. Further smoothness is ensured through a difference penalty on adjacent B-spline coefficients. Optimal lambda can be found by minimizing AIC, BIC or CV. The ps() function is used directly in the gam() function. Examples follow below. The plot.ps(), ps.wam(), gam.slist, and gam.wlist are support functions for glass(). The use of smooth ps() terms can co-exist with linear (polynomial) or factor terms in a GAM. However ps() terms cannot co-exist with other s() or lo() terms in a GAM. The default of glass() is identical to ps().

Examples:

source('ps.code')

# GAM logistic regression with 2 smooths and 1 linear and 1 block term

ps1 <- gam(binary.Y ~ ps(X1, ps.int=10, lambda=1, degree=2, order=2) + X2

+ ps(X3, ps.int=20, varying.index=t.col, lambda=10)

+ factor(block), family=binomial)

names(ps1)

ps1$nl.dim

ps1$lin.dim

aic1 <- ps1$dev + 2*(sum(ps1$nl.dim)+sum(ps1$lin.dim))

aic1

plot.ps(ps1, type=’smooth’, index=1, t.domain=X1)

plot.glass(glass1, type=’smooth’, index=2, t.domain=X3)

# Univariate Smoother using P-splines and a binomial response (same as ps( ))

pgam2 <- gam(Kyphosis ~ glass(Age, ps.int=20, lambda=.1, degree=4, order=4), family=binomial)

plot.ps(pgam2)