### From MathWiki

## SPIDA Day 4, June 2010
## Script for GLMMs and Splines
##
library(spida)
library(p3d)
?migraines
ds <- migraines
xqplot( ds )
ds$treat <- (ds$time > 0)*1 # start of treatment is day 0
fit <- glmmPQL ( ha ~ treat, data = ds,
random = ~ 1 | id,
family = binomial)
summary(fit)
names(ds)
# marginal relationship
xyplot( ha ~ time, ds, panel = function(x, y, ...) {
panel.xyplot( x, jitter(y),...)
panel.loess( x, y, ..., family = 'gaussian')
})
# create a measure of time that ranges within -10 to 10
# to help stability for spline
ds$tdays <- ds$time / 10
# create a spline function:
# with knots at 0, 5, 10, 30 and 50 days
# and
# quadratic in each interior interval and linear in the end intervals
# allowing a discontinuity in response at time 0
# and smooth joins at other knots
# Note:
# smooth = -1 = possibly discontinuous
# 0 = continuous but possibly kinky, i.e. kinky but no junp
# 1 = smooth of first order: no kink but possibly jerky
# - curvature[i.e. acceleration] can change
# 2 = change in acceleration, i.e. jerk but no snap
# 3 = snap but no crackle
# 4 = crackle but no pop
# 5 = pop but no (make up your own word)
# 6+ = you're on your own (as far as I know)
#
sp <- function(x) gsp( x,
knots = c(0,5,10,30,50)/10, # we express the knots in days but divide by 10
# to transform to 'dekadays'
degree = c(1,2,2, 2, 2 ,1),
smooth = c(-1,1,1, 1, 1))
fit <- glmmPQL ( ha ~ sp(tdays), data = ds,
random = ~ 1 | id,
family = binomial)
summary(fit)
pred <- expand.grid( time = seq(-30,50,.5))
pred$tdays <- pred$time / 10
pred$ha.logit <- predict( fit, pred, level = 0) # predict same as predict in glms
# for binomial gives the logit
xyplot( ha.logit ~ time, pred, type = 'l',lwd=2)
pred$ha.prob <- 1/(1+exp( -pred$ha.logit)) # transform to a probability
xyplot( ha.prob ~ time, pred, type = 'l',lwd=2)
# Let's ask questions:
?sc # spline contrast function to build L matrices for splines
sc(sp,x = c(-10,30)/10, D = 0)
zapsmall( sc(sp,x = c(-10,30)/10, D = 0) )
wald(fit)
Lp <- sc(sp,x = c(-10,30)/10, D = 0)
Lm <- cbind(0,Lp)
Ldiff <- rbind( Lm, diff= Lm[2,] - Lm[1,])
Ldiff
wald(fit, Ldiff)
## Are men and women different?
. . .