Lab Day 1.R

From MathWiki

#
#  SPIDA Mixed Models Day 1
#  Example of 'wald' tests with regular expressions
#  and visualizing model with xyplot and expand.grid
#
#

# definition of reorder function to make it work
  reorder.factor <- function (x, v, FUN = mean, ...) {
  factor(x, levels(x)[order(tapply(v, x, FUN, ...))])}


library(spida)
library(p3d)

dd <- hs1

fit <- lme( mathach ~ ses*Sector*Sex , dd, random = ~ 1 + ses | school)

# doesn't converge
# can try various things: increase maxIter, etc.
# But generally easiest to try 'nicer' random model:
# e.g. drop random slope or try centering x variable, globally or within groups
# 
dd$ses.cwg <- dd$ses - capply( dd$ses, dd$school, mean)

fit <- lme( mathach ~ ses*Sector*Sex, dd, random = ~ 1 + ses.cwg |school)
      # Eureka!

summary(fit)

# big model! could we drop Sector x Sex interaction?

# Could generate L matrix

Lm <- rbind( c( 0,0,0,0,0,0,1,0), c(0,0,0,0,0,0,0,1))
Lm

wald( fit, Lm)   # answer: ok to drop

# easier to use regular expression to match lines with Sector ... Sex

wald( fit, "Sector.*Sex")  # does the same thing

# Primer on regular expressions:
#  .    matches any one character
#  *    makes the previous character match 0 or more times
#  ?    makes the previous character match 0 or 1 time
#  +    makes the previous character match 1 or more times
#  ^    matches the start of a string
#  $    matches the end of a string
#  [LIST] matches any character in the LIST    
#  \\   used to escape the meaning of the next character if it has a special meaning
#       e.g. \\. is a literal period
#  abc|def   matches abc or def
#  [a-z]  matches any character in the range a-z, e.g. [A-Z] matches a capital letter 
#
#  Examples:
#  ":.*:"  matches all 3-way or higher-order interactions
#   

wald( fit, ":.*:")
wald( fit, "Sex")   # sex matters
wald( fit, ":Sex")   # can get rid of Sex interaction according to this

# refit without Sex interactions

fit2 <- lme( mathach ~ ses*Sector + Sex, dd, random = ~ 1 + ses.cwg |school)
summary(fit2)

# Visualize fit

quantile(dd$ses, c(.1,.9))
pred <- expand.grid( ses = seq(-1,1,.1), Sex = levels( dd$Sex), Sector = levels(dd$Sector))
pred$ma1 <- predict( fit2, pred, level = 0)

xyplot( ma1 ~ ses | Sector, pred, group = Sex, type = 'l',
    auto.key = T)
xyplot( ma1 ~ ses | Sector, pred, group = Sex, type = 'l',
    auto.key = list(lines=T,points=F))