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

'

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