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

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

“glass” <- function(x, ps.intervals=NULL, lambda=0, x.signal=NULL, signal.index=NULL, varying.index=NULL, degree=3, order=3, ridge.adj=0.00001, ridge.inv=0.00001)

· x: regressor of interest (also see signal.index and varying index below)

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

· lambda: non-negative regularization parameter for difference penalty

· x.signal: if SIGNAL component, this is the n X p signal matrix (p >> n possible)

· signal.index: if SIGNAL component, this the axis where B-splines are constructed (specify x above as 1:nrow(x.signal))

· varying.index: if VARYING component, this is the indexing variable where B-splines are constructed (e.g. time, age)

· 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 glass() function is used directly in the gam() function. Examples follow below. The plot.glass(), glass.wam(), gam.slist, and gam.wlist are support functions for glass(). The use of smooth glass() terms can co-exist with linear (polynomial) or factor terms in a GAM. However glass() terms cannot co-exist with other s() or lo() terms in a GAM. The default of glass() is identical to ps().

Examples:

source('glass.code')

# GLASS logistic regression with: 1 smooth component, 1 linear term, 1 varying coefficient term, 1 signal term, and 1 block term

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

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

+ glass(1:length(binary.Y), ps.int=20, signal.index=t.row,x.signal=S1.nXp, lambda=1)

+ factor(block), family=binomial)

names(glass1)

glass1$nl.dim

glass1$lin.dim

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

aic1

plot.glass(glass1, type=’smooth’, index=1, t.domain=X1)

plot.glass(glass1, type=’varying’, index=1, t.domain=t.col)

plot.glass(glass1, type=’signal’, index=1, t.domain=t.row)

# 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)