MATH 6643 Summer 2012 Applications of Mixed Models/Testing fixed and random effects in R

From Wiki1

Jump to: navigation, search

'

Testing Fixed Effects

'

 library(spidadev)   # loads car, nlme, etc.
 library(p3ddev)
 library(effects)   # need new version (perhaps on R-Forge)
 head(hs)
 dim(hs)
 xqplot(hs)
 
 hs$id <- factor(hs$school)
 
 fit.unc <- lme ( mathach ~ Sex*Sector*(ses+cvar(ses,id)), hs,
              random = ~ 1 + ses|id)
 fit.cwg <- lme ( mathach ~ Sex*Sector*(ses+cvar(ses,id)), hs,
                  random = ~ 1 + dvar(ses,id)|id)
 AIC( fit.unc, fit.cwg)
 
 fit <- fit.cwg
 summary(fit)
 

'

Overall fixed model

'

 wald(fit,-1)
 
 Anova(fit)  # Type II
 
 wald(fit,':')
 
 wald(fit, 'cvar')
 
 wald(fit, ":.*:")
 
 # Let's drop 3-way
 
 fit.2 <- update( fit, . ~ (Sex+Sector+(ses+cvar(ses,id)))^2)
 summary(fit.2)  # didn't intend to include ses x cvar(ses,id) interaction , but....
 
 wald(fit.2, "Sex")
 wald(fit.2, "Sex.*:")  # interactions not sig.
 
 fit.3 <- update( fit, . ~ Sex+(Sector+ses+cvar(ses,id))^2)
 fit.3 <- update( fit, . ~ Sex+(Sector+ses+cvar(ses,id))^2,
                  control = list(msMaxIter=500, msVerbose = T)) 
 summary(fit.3)  # Sector x cvar ? if we remove it we might load this on ses x cvar. Try it!
 
 effects(fit.3) # Try again with updated version
 # 
 # 
 #
 
 
 # Refit with explicit cvar
 
 hs$ses.m <- with( hs, cvar(ses,id))
 
 fit.3a <- update( fit.3, .~Sex+(Sector+ses+ses.m)^2)
 
 pred <- expand.grid( Sector =levels(hs$Sector), Sex = levels(hs$Sex), ses.m = -1:1, ses.d = -1:1)
 pred$ses <- with(pred, ses.m + ses.d)
 pred$idn <- factor(with(pred, paste( Sector, ses.m)))
 
 pred$mathach <- predict( fit.3a, newdata = pred, level = 0)
 
 
 xyplot( mathach ~ ses | Sector, pred, groups = idn:Sex, type = 'l',
         auto.key=T)
 
 # TO COME: error bands with wald
 # TO COME: complex L matrics with Leff
 
 # 
 
 # estimate contextual gap for Ses.m = 0 in both sectors
 L <- rbind( 'Catholic gap'= c(0,0,0,0, 1,0,0,0 ),
             'Public gap' = c(0,0,0,0, 1, 0, 1,0),
             "Difference" = c(0,0,0,0, 0,0,1,0))
 L
 wald(fit.3a, L)

'

Tests of Random effects

'


 fit.3
 # can we drop dvar(ses)
 # note this is REML fit
 fit.3i <- update(fit.3, random = ~ 1|id)
 anova(fit.3, fit.3i)  # looks like p = .968
 
 sim <- simulate( fit.3i, m2 = fit.3, nsim = 1000)
 sim
 plot(sim)   # but really more like p ~ .7
 
 # consider if nominal p = .15, really p closer to .05
 # if nominal p = .06 then really p < .001
 # need more sims for small p
Personal tools