# GLMM.R

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

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

. . .```