GLMM.R

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?

. . .