Computing and Plotting Differences in Effects
John Fox
May 31, 2018
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:
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':
##
## getDatahsfull$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)## [1] "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.
names(eff)## [1] "term" "formula"
## [3] "response" "variables"
## [5] "fit" "x"
## [7] "x.all" "model.matrix"
## [9] "data" "discrepancy"
## [11] "offset" "residuals"
## [13] "partial.residuals.range" "x.var"
## [15] "vcov" "se"
## [17] "lower" "upper"
## [19] "confidence.level" "transformation"
## [21] "family"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:
head(eff$fit)## [,1]
## 1 7.323443
## 2 7.467184
## 3 7.610926
## 4 7.754667
## 5 7.898409
## 6 8.042150dim(eff$fit)## [1] 246 1eff$vcov[1:5, 1:5] # first 5 rows and cols## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.1298149 0.1269729 0.1241310 0.1212890 0.1184471
## [2,] 0.1269729 0.1244152 0.1218574 0.1192997 0.1167419
## [3,] 0.1241310 0.1218574 0.1195839 0.1173103 0.1150368
## [4,] 0.1212890 0.1192997 0.1173103 0.1153210 0.1133316
## [5,] 0.1184471 0.1167419 0.1150368 0.1133316 0.1116265dim(eff$vcov)## [1] 246 246head(eff$x)## ses mean.ses Sector
## 1 -2.0 -1 Catholic
## 2 -1.9 -1 Catholic
## 3 -1.8 -1 Catholic
## 4 -1.7 -1 Catholic
## 5 -1.6 -1 Catholic
## 6 -1.5 -1 Catholictail(eff$x)## ses mean.ses Sector
## 241 1.5 1 Public
## 242 1.6 1 Public
## 243 1.7 1 Public
## 244 1.8 1 Public
## 245 1.9 1 Public
## 246 2.0 1 PublicThere 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:
fit <- as.vector(eff$fit) # reshape from 1-column matrix
V <- eff$vcov
x <- eff$x
(n <- 41*3*2)## [1] 246Next, 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):
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)## [1] 123 246L[1:2, ] # first 2 rows, note 1s in columns 1, 2 and -1s in columns 124, 125## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 1 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 1 0 0 0 0 0 0 0 0 0 0 0
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## [1,] 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0
## [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35]
## [1,] 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0
## [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46]
## [1,] 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0
## [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56] [,57]
## [1,] 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0
## [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66] [,67] [,68]
## [1,] 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0
## [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77] [,78] [,79]
## [1,] 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0
## [,80] [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88] [,89] [,90]
## [1,] 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0
## [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99] [,100] [,101]
## [1,] 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0
## [,102] [,103] [,104] [,105] [,106] [,107] [,108] [,109] [,110] [,111]
## [1,] 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0
## [,112] [,113] [,114] [,115] [,116] [,117] [,118] [,119] [,120] [,121]
## [1,] 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0
## [,122] [,123] [,124] [,125] [,126] [,127] [,128] [,129] [,130] [,131]
## [1,] 0 0 -1 0 0 0 0 0 0 0
## [2,] 0 0 0 -1 0 0 0 0 0 0
## [,132] [,133] [,134] [,135] [,136] [,137] [,138] [,139] [,140] [,141]
## [1,] 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0
## [,142] [,143] [,144] [,145] [,146] [,147] [,148] [,149] [,150] [,151]
## [1,] 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0
## [,152] [,153] [,154] [,155] [,156] [,157] [,158] [,159] [,160] [,161]
## [1,] 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0
## [,162] [,163] [,164] [,165] [,166] [,167] [,168] [,169] [,170] [,171]
## [1,] 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0
## [,172] [,173] [,174] [,175] [,176] [,177] [,178] [,179] [,180] [,181]
## [1,] 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0
## [,182] [,183] [,184] [,185] [,186] [,187] [,188] [,189] [,190] [,191]
## [1,] 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0
## [,192] [,193] [,194] [,195] [,196] [,197] [,198] [,199] [,200] [,201]
## [1,] 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0
## [,202] [,203] [,204] [,205] [,206] [,207] [,208] [,209] [,210] [,211]
## [1,] 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0
## [,212] [,213] [,214] [,215] [,216] [,217] [,218] [,219] [,220] [,221]
## [1,] 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0
## [,222] [,223] [,224] [,225] [,226] [,227] [,228] [,229] [,230] [,231]
## [1,] 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0
## [,232] [,233] [,234] [,235] [,236] [,237] [,238] [,239] [,240] [,241]
## [1,] 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0
## [,242] [,243] [,244] [,245] [,246]
## [1,] 0 0 0 0 0
## [2,] 0 0 0 0 0diff <- as.vector(L %*% fit)
se <- sqrt(diag(L %*% V %*% t(L)))
head(diff)## [1] 3.723754 3.588928 3.454101 3.319275 3.184448 3.049622length(diff)## [1] 123head(se)## [1] 0.4686511 0.4582823 0.4488019 0.4402673 0.4327345 0.4262567length(se)## [1] 123The 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:
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)## ses mean.ses diff lower upper
## 1 -2.0 -1 3.723754 2.805198 4.642310
## 2 -1.9 -1 3.588928 2.690694 4.487161
## 3 -1.8 -1 3.454101 2.574450 4.333753
## 4 -1.7 -1 3.319275 2.456351 4.182199
## 5 -1.6 -1 3.184448 2.336289 4.032608
## 6 -1.5 -1 3.049622 2.214159 3.885085tail(Diffs)## ses mean.ses diff lower upper
## 118 1.5 1 -0.5222436 -1.367804 0.32331711
## 119 1.6 1 -0.6570701 -1.515178 0.20103786
## 120 1.7 1 -0.7918965 -1.664600 0.08080738
## 121 1.8 1 -0.9267229 -1.815971 -0.03747523
## 122 1.9 1 -1.0615493 -1.969182 -0.15391646
## 123 2.0 1 -1.1963757 -2.124126 -0.26862578dim(Diffs)## [1] 123 5Finally, I use xyplot() in the lattice package to plot the differences along with confidence bands:
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.