author: Nancy Reid
date: March 4
library(faraway)
data(jsp)
head(jsp)
## school class gender social raven id english math year
## 1 1 1 girl 9 23 1 72 23 0
## 2 1 1 girl 9 23 1 80 24 1
## 3 1 1 girl 9 23 1 39 23 2
## 4 1 1 boy 2 15 2 7 14 0
## 5 1 1 boy 2 15 2 17 11 1
## 6 1 1 boy 2 22 3 88 36 0
jspr = jsp[jsp$year == 2,]
========================================================= The data =========================================================
names(jspr)
## [1] "school" "class" "gender" "social" "raven" "id" "english"
## [8] "math" "year"
“Raven” is test score on entering (measure of ability?)
“Social” is measure of socio-economic status
50 schools; up to 4 classes in each school;
glin <- lm(math ~ raven * gender * social, data = jspr)
anova(glin)
## Analysis of Variance Table
##
## Response: math
## Df Sum Sq Mean Sq F value Pr(>F)
## raven 1 11480.5 11480.5 368.0625 < 2.2e-16 ***
## gender 1 44.1 44.1 1.4142 0.234668
## social 8 779.4 97.4 3.1233 0.001725 **
## raven:gender 1 0.0 0.0 0.0004 0.984718
## raven:social 8 582.6 72.8 2.3347 0.017460 *
## gender:social 8 450.1 56.3 1.8038 0.072742 .
## raven:gender:social 8 234.6 29.3 0.9400 0.482355
## Residuals 917 28602.8 31.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lm(math ~ raven + social, data = jspr))
##
## Call:
## lm(formula = math ~ raven + social, data = jspr)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.8430 -3.2426 0.7726 3.7765 14.0825
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.02481 1.37451 12.386 <2e-16 ***
## raven 0.58040 0.03256 17.826 <2e-16 ***
## social2 0.04950 1.12938 0.044 0.9651
## social3 -0.42893 1.19568 -0.359 0.7199
## social4 -1.77452 1.05993 -1.674 0.0944 .
## social5 -0.78228 1.18924 -0.658 0.5108
## social6 -2.49373 1.26094 -1.978 0.0483 *
## social7 -3.04851 1.29065 -2.362 0.0184 *
## social8 -3.11746 1.77494 -1.756 0.0793 .
## social9 -0.63278 1.12731 -0.561 0.5747
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.632 on 943 degrees of freedom
## Multiple R-squared: 0.2907, Adjusted R-squared: 0.2839
## F-statistic: 42.93 on 9 and 943 DF, p-value: < 2.2e-16
library(lme4)
## Loading required package: Matrix
## Loading required package: Rcpp
schoolmod <- lmer(math ~ raven*social*gender +
(1 | school) + (1 | school:class),
data = jspr)
anova(schoolmod)
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## raven 1 10217.9 10217.9 374.4012
## social 8 615.7 77.0 2.8203
## gender 1 21.6 21.6 0.7917
## raven:social 8 577.5 72.2 2.6449
## raven:gender 1 2.5 2.5 0.0902
## social:gender 8 275.3 34.4 1.2608
## raven:social:gender 8 187.2 23.4 0.8572
jspr$craven <- jspr$raven - mean(jspr$raven)
schoolmod2 <- lmer(math ~ craven*social + (1|school) + (1 | school:class), data = jspr)
summary(schoolmod2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: math ~ craven * social + (1 | school) + (1 | school:class)
## Data: jspr
##
## REML criterion at convergence: 5921.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9872 -0.5250 0.1468 0.6237 2.6341
##
## Random effects:
## Groups Name Variance Std.Dev.
## school:class (Intercept) 1.177 1.085
## school (Intercept) 3.148 1.774
## Residual 27.141 5.210
## Number of obs: 953, groups: school:class, 90; school, 48
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 31.91123 1.19554 26.692
## craven 0.60584 0.18854 3.213
## social2 0.02361 1.27220 0.019
## social3 -0.63067 1.30888 -0.482
## social4 -1.96700 1.19707 -1.643
## social5 -1.35847 1.30023 -1.045
## social6 -2.26870 1.37374 -1.651
## social7 -2.55185 1.40555 -1.816
## social8 -3.39495 1.80137 -1.885
## social9 -0.83128 1.25347 -0.663
## craven:social2 -0.13206 0.20579 -0.642
## craven:social3 -0.22433 0.21888 -1.025
## craven:social4 0.03582 0.19488 0.184
## craven:social5 -0.15034 0.20889 -0.720
## craven:social6 -0.03860 0.23259 -0.166
## craven:social7 0.39826 0.23177 1.718
## craven:social8 0.25600 0.26154 0.979
## craven:social9 -0.08102 0.20550 -0.394
##
## Correlation of Fixed Effects:
## (Intr) craven socil2 socil3 socil4 socil5 socil6 socil7 socil8
## craven -0.513
## social2 -0.879 0.479
## social3 -0.860 0.469 0.805
## social4 -0.942 0.510 0.882 0.862
## social5 -0.863 0.468 0.807 0.784 0.863
## social6 -0.826 0.446 0.769 0.754 0.829 0.754
## social7 -0.810 0.435 0.750 0.739 0.810 0.743 0.712
## social8 -0.627 0.348 0.587 0.576 0.627 0.573 0.546 0.539
## social9 -0.901 0.489 0.839 0.823 0.901 0.823 0.791 0.775 0.596
## craven:scl2 0.476 -0.910 -0.519 -0.435 -0.476 -0.436 -0.417 -0.404 -0.321
## craven:scl3 0.445 -0.860 -0.414 -0.417 -0.442 -0.405 -0.387 -0.375 -0.302
## craven:scl4 0.499 -0.966 -0.466 -0.456 -0.492 -0.456 -0.433 -0.424 -0.337
## craven:scl5 0.463 -0.899 -0.435 -0.424 -0.461 -0.396 -0.403 -0.393 -0.313
## craven:scl6 0.413 -0.810 -0.381 -0.377 -0.407 -0.374 -0.297 -0.344 -0.279
## craven:scl7 0.417 -0.813 -0.389 -0.381 -0.412 -0.377 -0.363 -0.295 -0.283
## craven:scl8 0.374 -0.718 -0.352 -0.346 -0.373 -0.343 -0.325 -0.324 -0.242
## craven:scl9 0.471 -0.916 -0.441 -0.432 -0.470 -0.430 -0.410 -0.399 -0.322
## socil9 crvn:2 crvn:3 crvn:4 crvn:5 crvn:6 crvn:7 crvn:8
## craven
## social2
## social3
## social4
## social5
## social6
## social7
## social8
## social9
## craven:scl2 -0.455
## craven:scl3 -0.425 0.786
## craven:scl4 -0.476 0.881 0.832
## craven:scl5 -0.441 0.821 0.775 0.871
## craven:scl6 -0.394 0.736 0.700 0.783 0.728
## craven:scl7 -0.397 0.740 0.699 0.786 0.733 0.658
## craven:scl8 -0.359 0.657 0.622 0.695 0.647 0.580 0.587
## craven:scl9 -0.428 0.834 0.788 0.886 0.826 0.743 0.743 0.656
par(mfrow=c(1,2))
qqnorm(resid(schoolmod2), main = "")
plot(fitted(schoolmod2), resid(schoolmod2),xlab = "fitted", ylab = "residuals")