---
title: "Computing and Plotting *Differences* in Effects"
author: "John Fox"
date: "May 31, 2018"
output: html_document
---

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


As before, I'll fit Georges's model to the full HSB ("High School and Beyond" study) data set, but this time I'll save the object produced by the `Effect()` function:
```{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 )
eff <- Effect(c("ses", "mean.ses", "Sector"), fitc.full,  
            xlevels=list(mean.ses=c(-1, 0, 1), 
                         ses=seq(-2, 2, by=0.1)))
class(eff)
```

This produces an object of class `"eff"`, which includes a great deal of information. We can get a sense of what's there by printing the names of the components of the object. Also see either `str(eff)` (not shown) or the RStudio Environment tab and `?Effect`.
```{r}
names(eff)
```

As it turns out, the estimated effect is in `fit`, its covariance matrix is in `vcov`, and the various combinations of the predictors for the effect are in `x`:
```{r}
head(eff$fit)
dim(eff$fit)
eff$vcov[1:5, 1:5] # first 5 rows and cols
dim(eff$vcov)
head(eff$x)
tail(eff$x)
```

There are 246 rows because there are 41 values of `ses` between -2 and +2, times 3 values of `mean.ses`, -1, 0, and +1, times 2 levels of `Sector`, `"Catholic"` and `"Public"`.

I'll save these elements in global variables to make it easier to compute the differences in effects between the sectors and the standard errors of these differences:
```{r}
fit <- as.vector(eff$fit) # reshape from 1-column matrix
V <- eff$vcov
x <- eff$x
(n <- 41*3*2)
```

Next, I'll construct an **L** matrix to take differences between corresponding Catholic fits (in the first 132 values) and Public fits (in the last 132 values):
```{r}
L <- matrix(0, nrow=n/2, ncol=n)
for (i in 1:(n/2)){
  L[i, c(i, i + (n/2))] <- c(1, -1)
}
dim(L)
L[1:2, ] # first 2 rows, note 1s in columns 1, 2 and -1s in columns 124, 125

diff <- as.vector(L %*% fit)
se <- sqrt(diag(L %*% V %*% t(L)))
head(diff)
length(diff)
head(se)
length(se)
```

The differences are computed as **d = L e** (if **e** is the vector of effects), and the standard errors of the differences are the square-root diagonal elements of **L V L'** (where **V** is the covariance matrix of the effects and **L'** is the transpose of **L**).

Then I compute confidence bounds for the differences, and put the differences and confidence bounds into a data frame that starts with the first 123 rows of `x`, which include all combinations of `ses` and `mean.ses`:
```{r}
lower <- diff - 1.96*se
upper <- diff + 1.96*se
Diffs <- x[1:(n/2), ]
Diffs$Sector <- NULL # no longer relevant (all Catholic in first 123 rows)
Diffs <- cbind(Diffs, diff, lower, upper)
head(Diffs)
tail(Diffs)
dim(Diffs)
```

Finally, I use `xyplot()` in the **lattice** package to plot the differences along with confidence bands:
```{r}
library("lattice")
Diffs$mean.ses <- as.factor(Diffs$mean.ses) # a trick to get "strip" labels with values
xyplot(diff ~ ses | mean.ses, data=Diffs, lower=Diffs$lower, upper=Diffs$upper,
       panel=function(x, y, lower, upper, subscripts){
         llines(x, y, lwd=2) # difference in effects
         panel.abline(c(0, 0), lty=2, col="black") # horizontal line at y = 0
         panel.polygon(c(x, rev(x)), c(upper[subscripts], rev(lower[subscripts])),
                       col="magenta", border = FALSE, alpha=0.2) # confidence band
       },
       strip=function(...) strip.default(..., sep=" = ", strip.names=c(TRUE, TRUE)),
       ylab="Difference (Catholic - Public)")

```

Thus, `Effect()` takes care of generating combinations of the predictors and computing effects and their covariance matrix, but to get differences and their standard errors, and to plot the results require custom computations.