Fitting Georges’s model to the full HSB (“High School and Beyond” study) data set:

library("spida2")  
library("effects")
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
library("car")
library("nlme")
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:spida2':
## 
##     getData
hsfull$mean.ses <- with(hsfull, cvar(ses, school))

fitc.full <- lme( mathach ~ (ses + mean.ses)*Sector, data=hsfull,
             random = ~ 1 + ses | school )
S(fitc.full)
## Linear mixed model fit by REML,  Data: hsfull
## 
## Fixed Effects:
##  Formula: mathach ~ (ses + mean.ses) * Sector 
## 
##                       Estimate Std.Error   df t value Pr(>|t|)    
## (Intercept)            13.4519    0.2256 7023  59.629  < 2e-16 ***
## ses                     1.4374    0.1669 7023   8.614  < 2e-16 ***
## mean.ses                3.2536    0.5569  156   5.842 2.91e-08 ***
## SectorPublic           -1.2637    0.3052  156  -4.140 5.66e-05 ***
## ses:SectorPublic        1.3483    0.2228 7023   6.052 1.50e-09 ***
## mean.ses:SectorPublic  -0.2365    0.7636  156  -0.310    0.757    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Random effects:
##  Formula: ~1 + ses | school
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept) 1.5358681 (Intr)
## ses         0.2648212 0.435 
## Residual    6.0657219       
## 
## 
## Number of Observations: 7185
## Number of Groups: 160 
## 
##    logLik        df       AIC       BIC 
## -23256.33        10  46532.67  46601.46
Anova(fitc.full)
## Analysis of Deviance Table (Type II tests)
## 
## Response: mathach
##                    Chisq Df Pr(>Chisq)    
## ses             393.8001  1  < 2.2e-16 ***
## mean.ses         67.3858  1  2.233e-16 ***
## Sector           20.6885  1  5.404e-06 ***
## ses:Sector       36.6277  1  1.430e-09 ***
## mean.ses:Sector   0.0959  1     0.7568    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Normally one can have functions calls on the right-hand side of a model that’s graphed by the effects package, but the term cvar(ses, school) would confuse the Effect() function – it would treat it as a function of the two “predictors” ses and school rather than as a single predictor. Consequently, I created the new variable mean.ses for school-mean SES in the data set.

The default effect plot for the three predictors:

plot(Effect(c("ses", "mean.ses", "Sector"), fitc.full))

A customized version of the graph:

plot(Effect(c("ses", "mean.ses", "Sector"), fitc.full,  
            xlevels=list(mean.ses=c(-1, 0, 1),
                         ses=seq(-2, 2, by=0.1))),
     lattice=list(layout=c(3, 2)),
     axes=list(x=list(rug=FALSE)))

Putting the fitted lines for Catholic and Public schools in the same panel:

plot(Effect(c("ses", "mean.ses", "Sector"), fitc.full,  
            xlevels=list(mean.ses=c(-1, 0, 1), 
                         ses=seq(-2, 2, by=0.1))),
     lines=list(multiline=TRUE),
     lattice=list(layout=c(3, 1)),
     confint=list(style="bands"),
     axes=list(x=list(rug=FALSE)))