---
title: "Some Effect Plots for the HSB Data"
author: "John Fox"
date: "May 30, 2018"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


Fitting Georges's model to the full HSB ("High School and Beyond" study) data set:
```{r}
library("spida2")  
library("effects")
library("car")
library("nlme")

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)
Anova(fitc.full)
```

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:
```{r fig.width=8}
plot(Effect(c("ses", "mean.ses", "Sector"), fitc.full))
```

A customized version of the graph:
```{r}
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:
```{r fig.height=4}
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)))
```


