library(lme4) dat1 <- read.table("lecture10ex1.csv", sep=",", header=TRUE) psych::describe(dat1) m10 <- lmer(math ~ 1 + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE) summary(m10) m11 <- lmer(iq ~ 1 + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")) summary(m11) dat1$iqc <- (dat1$iq - 100)/15 dat1$aveschooliqc <- ave(dat1$iqc, dat1$schoolid) dat1$diffschooliqc <- dat1$iqc - dat1$aveschooliqc dat1school <- dat1[!duplicated(dat1$schoolid),] dat1school$avecountryiqc <- ave(dat1school$aveschooliqc, dat1school$countryid) dat1school$diffcountryiqc <- dat1school$aveschooliqc - dat1school$avecountryiqc dat1school$aveprivate <- ave(dat1school$private, dat1school$countryid) dat1school$diffprivate <- dat1school$private - dat1school$aveprivate matchit <- match(dat1$schoolid, dat1school$schoolid) dat1$avecountryiqc <- dat1school[matchit, "avecountryiqc"] dat1$diffcountryiqc <- dat1school[matchit, "diffcountryiqc"] dat1$aveprivate <- dat1school[matchit, "aveprivate"] dat1$diffprivate <- dat1school[matchit, "diffprivate"] m1xa1 <- lmer(math ~ 1 + diffschooliqc + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")) summary(m1xa1) m1xa2 <- lmer(math ~ 1 + iqc + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")) summary(m1xa2) m1xb1 <- lmer(math ~ 1 + diffschooliqc + diffcountryiqc + avecountryiqc + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")) summary(m1xb1) m1xb2 <- lmer(math ~ 1 + diffschooliqc + aveschooliqc + avecountryiqc + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")) summary(m1xb2) m1xb3 <- lmer(math ~ 1 + iqc + aveschooliqc + avecountryiqc + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")) summary(m1xb3) m1wa1 <- lmer(math ~ 1 + diffprivate + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")) summary(m1wa1) m1wa2 <- lmer(math ~ 1 + private + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")) summary(m1wa2) m1wb1 <- lmer(math ~ 1 + diffprivate + aveprivate + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")) summary(m1wb1) m1wb2 <- lmer(math ~ 1 + private + aveprivate + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")) summary(m1wb2) dat1country <- dat1[!duplicated(dat1$countryid),] apply(dat1country, 2, mean) m12 <- lmer(math ~ 1 + diffschooliqc + diffcountryiqc + I(avecountryiqc+0.691) + diffprivate + I(aveprivate-0.506) + I(quality-53.33) + I(opportunity-47.98) + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")) summary(m12) m13 <- lmer(math ~ 1 + diffschooliqc + diffcountryiqc + I(avecountryiqc+0.691) + diffprivate + I(aveprivate-0.506) + I(quality-53.33) + I(opportunity-47.98) + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")) m14 <- lmer(math ~ 1 + diffschooliqc + diffcountryiqc + I(avecountryiqc+0.691) + diffprivate + I(aveprivate-0.506) + I(quality-53.33) + I(opportunity-47.98) + (1 + diffschooliqc|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead", calc.derivs=FALSE)) anova(m13, m14) m15 <- lmer(math ~ 1 + diffschooliqc + diffcountryiqc + I(avecountryiqc+0.691) + diffprivate + I(aveprivate-0.506) + I(quality-53.33) + I(opportunity-47.98) + (1|schoolid) + (1 + diffschooliqc|countryid), data=dat1, REML=FALSE) anova(m13, m15) m16 <- lmer(math ~ 1 + diffschooliqc + diffcountryiqc + I(avecountryiqc+0.691) + diffprivate + I(aveprivate-0.506) + I(quality-53.33) + I(opportunity-47.98) + (1|schoolid) + (1 + diffschooliqc + diffcountryiqc|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead", calc.derivs=FALSE)) anova(m15, m16) m17 <- lmer(math ~ 1 + diffschooliqc + diffcountryiqc + I(avecountryiqc+0.691) + diffprivate + I(aveprivate-0.506) + I(quality-53.33) + I(opportunity-47.98) + (1|schoolid) + (1 + diffschooliqc + diffprivate|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")) anova(m15, m17) summary(m17) ranefm17 <- ranef(m17) names(ranefm17) ranefm17country <- ranefm17$countryid plot(ranefm17country) m18 <- lmer(math ~ 1 + diffschooliqc + diffcountryiqc + I(avecountryiqc+0.691) + diffprivate + I(aveprivate-0.506) + I(quality-53.33) + I(opportunity-47.98) + diffschooliqc*diffcountryiqc + diffschooliqc*I(avecountryiqc+0.691) + diffschooliqc*diffprivate + diffschooliqc*I(aveprivate-0.506) + diffschooliqc*I(quality-53.33) + diffschooliqc*I(opportunity-47.98) + diffprivate*I(avecountryiqc+0.691) + diffprivate*I(aveprivate-0.506) + diffprivate*I(quality-53.33) + diffprivate*I(opportunity-47.98) + (1|schoolid) + (1 + diffschooliqc + diffprivate|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")) summary(m18) mean(dat1$diffschooliqc) sd(dat1$diffschooliqc) mean(dat1country$quality) sd(dat1country$quality) mean(dat1country$opportunity) sd(dat1country$opportunity) sumoutm18 <- summary(m18) coef(sumoutm18) psi <- coef(sumoutm18)[,1] vpsi <- vcov(sumoutm18) # diffschooliqc -> math at each level of quality cvec1 <- c(1, -0.6, rep(0, 16)) cvec2 <- c(1, 0, rep(0, 16)) cvec3 <- c(1, 0.6, rep(0, 16)) cvec4 <- c(rep(0, 6), 1, rep(0, 5), -0.6, rep(0, 5)) cvec5 <- c(rep(0, 6), 1, rep(0, 5), 0, rep(0, 5)) cvec6 <- c(rep(0, 6), 1, rep(0, 5), 0.6, rep(0, 5)) cmat <- cbind(cvec1, cvec2, cvec3, cvec4, cvec5, cvec6) est <- t(cmat) %*% psi se <- sqrt(diag(t(cmat) %*% vpsi %*% cmat)) zval <- est/se pval <- round(pnorm(-abs(zval)) * 2, 6) cbind(est, se, zval, pval) # quality -> math at each level of diffschooliqc cvec1 <- c(1, rep(0, 5), -13.33, rep(0, 11)) cvec2 <- c(1, rep(0, 5), -0.33, rep(0, 11)) cvec3 <- c(1, rep(0, 5), 12.67, rep(0, 11)) cvec4 <- c(0, 1, rep(0, 10), -13.33, rep(0, 5)) cvec5 <- c(0, 1, rep(0, 10), -0.33, rep(0, 5)) cvec6 <- c(0, 1, rep(0, 10), 12.67, rep(0, 5)) cmat <- cbind(cvec1, cvec2, cvec3, cvec4, cvec5, cvec6) est <- t(cmat) %*% psi se <- sqrt(diag(t(cmat) %*% vpsi %*% cmat)) zval <- est/se pval <- round(pnorm(-abs(zval)) * 2, 6) cbind(est, se, zval, pval) # opportunity -> math at each level of private cvec1 <- c(1, rep(0, 3), 0.5, rep(0, 13)) cvec2 <- c(1, rep(0, 3), -0.5, rep(0, 13)) cvec3 <- c(rep(0, 7), 1, rep(0, 9), 0.5) cvec4 <- c(rep(0, 7), 1, rep(0, 9), -0.5) cmat <- cbind(cvec1, cvec2, cvec3, cvec4) est <- t(cmat) %*% psi se <- sqrt(diag(t(cmat) %*% vpsi %*% cmat)) zval <- est/se pval <- round(pnorm(-abs(zval)) * 2, 6) cbind(est, se, zval, pval) # private -> math at each level of opportunity cvec1 <- c(1, rep(0, 6), -7.98, rep(0, 10)) cvec2 <- c(1, rep(0, 6), 0.02, rep(0, 10)) cvec3 <- c(1, rep(0, 6), 8.02, rep(0, 10)) cvec4 <- c(rep(0, 4), 1, rep(0, 12), -7.98) cvec5 <- c(rep(0, 4), 1, rep(0, 12), 0.02) cvec6 <- c(rep(0, 4), 1, rep(0, 12), 8.02) cmat <- cbind(cvec1, cvec2, cvec3, cvec4, cvec5, cvec6) est <- t(cmat) %*% psi se <- sqrt(diag(t(cmat) %*% vpsi %*% cmat)) zval <- est/se pval <- round(pnorm(-abs(zval)) * 2, 6) cbind(est, se, zval, pval) ############################################################## dat2 <- read.table("lecture10ex2.csv", sep=",", header=TRUE) psych::describe(dat2) dat2$yearc <- dat2$year - 2017 m21 <- lmer(test ~ 1 + yearc + (1|timeid), data=dat2, REML=FALSE) m22 <- lmer(test ~ 1 + yearc + (1|timeid) + (1|schoolid), data=dat2, REML=FALSE) anova(m21, m22) summary(m22) m23 <- lmer(test ~ 1 + yearc + (1|timeid) + (1 + yearc|schoolid), data=dat2, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")) anova(m22, m23) summary(m23) plot(ranef(m23)$schoolid, xlab="School Averaged Test Score in 2017", ylab="School Change in Test Score per Year") m24 <- lmer(test ~ 1 + yearc*public + (1|timeid) + (1 + yearc|schoolid), data=dat2, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")) summary(m24) sumoutm24 <- summary(m24) coef(sumoutm24) psi <- coef(sumoutm24)[,1] vpsi <- vcov(sumoutm24) cvec1 <- c(1, 0, 1, 0) cvec2 <- c(1, 0, 0, 0) cvec3 <- c(0, 1, 0, 1) cvec4 <- c(0, 1, 0, 0) cmat <- cbind(cvec1, cvec2, cvec3, cvec4) est <- t(cmat) %*% psi se <- sqrt(diag(t(cmat) %*% vpsi %*% cmat)) zval <- est/se pval <- round(pnorm(-abs(zval)) * 2, 6) cbind(est, se, zval, pval) cvec1 <- c(1, 0, 0, 0) cvec2 <- c(1, 1, 0, 0) cvec3 <- c(1, 2, 0, 0) cvec4 <- c(0, 0, 1, 0) cvec5 <- c(0, 0, 1, 1) cvec6 <- c(0, 0, 1, 2) cmat <- cbind(cvec1, cvec2, cvec3, cvec4, cvec5, cvec6) est <- t(cmat) %*% psi se <- sqrt(diag(t(cmat) %*% vpsi %*% cmat)) zval <- est/se pval <- round(pnorm(-abs(zval)) * 2, 6) cbind(est, se, zval, pval) ###### R-squared library(r2mlm) # Clustered-mean-centered m18 <- lmer(math ~ 1 + diffschooliqc + diffcountryiqc + I(avecountryiqc+0.691) + diffprivate + I(aveprivate-0.506) + I(quality-53.33) + I(opportunity-47.98) + diffschooliqc*diffcountryiqc + diffschooliqc*I(avecountryiqc+0.691) + diffschooliqc*diffprivate + diffschooliqc*I(aveprivate-0.506) + diffschooliqc*I(quality-53.33) + diffschooliqc*I(opportunity-47.98) + diffprivate*I(avecountryiqc+0.691) + diffprivate*I(aveprivate-0.506) + diffprivate*I(quality-53.33) + diffprivate*I(opportunity-47.98) + (1|schoolid) + (1 + diffschooliqc + diffprivate|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")) dat1_1 <- model.matrix(m18) sumoutm18 <- summary(m18) lv1_covs <- c("diffschooliqc", "diffschooliqc:diffcountryiqc", "diffschooliqc:I(avecountryiqc + 0.691)", "diffschooliqc:diffprivate", "diffschooliqc:I(aveprivate - 0.506)", "diffschooliqc:I(quality - 53.33)", "diffschooliqc:I(opportunity - 47.98)") lv2_covs <- c("diffschooliqc", "diffprivate", "I(avecountryiqc + 0.691):diffprivate", "diffprivate:I(aveprivate - 0.506)", "diffprivate:I(quality - 53.33)", "diffprivate:I(opportunity - 47.98)") lv3_covs <- c("I(avecountryiqc + 0.691)", "I(aveprivate - 0.506)", "I(quality - 53.33)", "I(opportunity - 47.98)") ranefvar18 <- as.matrix(Matrix::bdiag(VarCorr(m18))) r2mlm3_manual(data=dat1_1, l1_covs=lv1_covs, l2_covs=lv2_covs, l3_covs=lv3_covs, random_covs12 = NULL, random_covs13 = "diffschooliqc", random_covs23 = "diffprivate", gamma_1 = coef(sumoutm18)[lv1_covs, "Estimate"], gamma_2 = coef(sumoutm18)[lv2_covs, "Estimate"], gamma_3 = coef(sumoutm18)[lv3_covs, "Estimate"], Tau12 = ranefvar18[1, 1,drop=FALSE], Tau13 = ranefvar18[2:3, 2:3, drop=FALSE], Tau23 = ranefvar18[4, 4, drop=FALSE], sigma2=getME(m18, "sigma")^2, clustermeancentered = TRUE) # Non-clustered-mean-centered m24 <- lmer(test ~ 1 + yearc*public + (1|timeid) + (1 + yearc|schoolid), data=dat2, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")) dat2_1 <- model.matrix(m24) varnames <- colnames(dat2_1) dat2_1 <- data.frame(dat2_1_1, timeid = dat2$timeid,schoolid = dat2$schoolid) colnames(dat2_1)[1:length(varnames)] <- varnames sumoutm24 <- summary(m24) lv2_covs <- c("yearc", "yearc:public") lv3_covs <- c("public") ranefvar24 <- as.matrix(Matrix::bdiag(VarCorr(m24))) r2mlm3_manual(data=dat2_1, l1_covs=NULL, l2_covs=lv2_covs, l3_covs=lv3_covs, random_covs12 = NULL, random_covs13 = NULL, random_covs23 = "yearc", gamma_1 = NULL, gamma_2 = coef(sumoutm24)[lv2_covs, "Estimate"], gamma_3 = coef(sumoutm24)[lv3_covs, "Estimate"], Tau2_noncmc = ranefvar24[1, 1,drop=FALSE], Tau3_noncmc = ranefvar24[2:3, 2:3, drop=FALSE], sigma2=getME(m24, "sigma")^2, l2clusterID_noncmc = "timeid", l3clusterID_noncmc = "schoolid", clustermeancentered = FALSE)