dat3 <-read.table("lecture4ex3.csv", sep=",", header=TRUE) dat3$aveage <- ave(dat3$age, dat3$tableid) dat3$diffage <- dat3$age - dat3$aveage library(lme4) out3m1 <- lmer(sat ~ 1 + age + (1|tableid), data=dat3, REML=FALSE) summary(out3m1) out3m2 <- lmer(sat ~ 1 + diffage + (1|tableid), data=dat3, REML=FALSE) summary(out3m2) out3m3 <- lmer(sat ~ 1 + age + aveage + (1|tableid), data=dat3, REML=FALSE) summary(out3m3) out3m4 <- lmer(sat ~ 1 + diffage + aveage + (1|tableid), data=dat3, REML=FALSE) summary(out3m4) dat3$sage <- dat3$age / 10 dat3$avesage <- ave(dat3$sage, dat3$tableid) dat3$diffsage <- dat3$sage - dat3$avesage out3m1 <- lmer(sat ~ 1 + sage + (1|tableid), data=dat3, REML=FALSE) summary(out3m1) out3m2 <- lmer(sat ~ 1 + diffsage + (1|tableid), data=dat3, REML=FALSE) summary(out3m2) out3m3 <- lmer(sat ~ 1 + sage + avesage + (1|tableid), data=dat3, REML=FALSE) summary(out3m3) out3m4 <- lmer(sat ~ 1 + diffsage + avesage + (1|tableid), data=dat3, REML=FALSE) summary(out3m4) out3m5 <- lmer(sat ~ 1 + sage + avesage + sage:avesage + (1 + sage|tableid), data=dat3, REML=FALSE) summary(out3m5) out3m6 <- lmer(sat ~ 1 + diffsage + avesage + diffsage:avesage + (1 + diffsage|tableid), data=dat3, REML=FALSE) summary(out3m6) dat4 <-read.table("lecture6ex1.csv", sep=",", header=TRUE) dat4$aveach <- ave(dat4$ach, dat4$schoolid) dat4$diffach <- dat4$ach - dat4$aveach out4m1 <- lmer(efficacy ~ 1 + diffach + I(aveach - 50) + private + diffach:I(aveach - 50) + (1 + diffach|schoolid), data=dat4, REML=FALSE) summary(out4m1) out4m0 <- lmer(ach ~ 1 + (1|schoolid), data=dat4, REML=FALSE) summary(out4m0) plot(dat4$diffach, dat4$efficacy, type="n", xlab="Math Achievement Deviation within School", ylab="Self Efficacy") abline(a = 35.35, b = 0.01, col="red") abline(a = 35.85, b = 0.11, col="green") abline(a = 34.85, b = 0.21, col="blue") legend(-40, 50, c("Low", "Medium", "High"), col=c("red", "green", "blue"), lty=1, title="School Math Achievement") plot(dat4$ach, dat4$efficacy, type="n", xlab="School Math Achievement", ylab="Self Efficacy") abline(a = 35 - (-0.15*50), b = -0.15, col="red") abline(a = 35.6 - (-0.05*50), b = -0.05, col="green") abline(a = 36.2 - (0.05*50), b = 0.05, col="blue") legend(20, 30, c("Low", "Medium", "High"), col=c("red", "green", "blue"), lty=1, title="Relative Position within School") sumout4 <- summary(out4m1) psi <- coef(sumout4)[,1] vpsi <- vcov(sumout4) cvec1 <- c(1, 0, -5, 0, 0) cvec2 <- c(1, 0, 5, 0, 0) cvec3 <- c(1, 0, 15, 0, 0) cvec4 <- c(0, 1, 0, 0, -5) cvec5 <- c(0, 1, 0, 0, 5) cvec6 <- c(0, 1, 0, 0, 15) 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) cvec1 <- c(1, -10, 0, 0, 0) cvec2 <- c(1, 0, 0, 0, 0) cvec3 <- c(1, 10, 0, 0, 0) cvec4 <- c(0, 0, 1, 0, -10) cvec5 <- c(0, 0, 1, 0, 0) cvec6 <- c(0, 0, 1, 0, 10) 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) ####################################### dat3$avefemale <- ave(dat3$female, dat3$tableid) out3g1 <- lmer(sat ~ 1 + female + avefemale + (1|tableid), data=dat3, REML=FALSE) out3g2 <- lmer(sat ~ 1 + female + avefemale + (1 + female|tableid), data=dat3, REML=FALSE) anova(out3g1, out3g2) out3g3 <- lmer(sat ~ 1 + female + avefemale + female:avefemale + (1 + female|tableid), data=dat3, REML=FALSE) anova(out3g2, out3g3) summary(out3g3) avefemale <- seq(0, 1, 0.01) cust <- 66.5712 + ((-3.8363 - 3.646) * avefemale) + (4.6837*avefemale*avefemale) plot(avefemale, cust, type="l", xlim=c(0, 1), ylim=c(60, 70), xlab="Female Proportion in a Table", ylab="Customer Satisfaction") plot(dat3$female, dat3$sat, type="n", xlab="Female", ylab="Customer Satisfaction") abline(a = 65.84, b = -2.904, col="red") abline(a = 64.75, b = -1.5, col="green") abline(a = 63.65, b = -0.096, col="blue") legend(0.01, 50, c("20%", "50%", "80%"), col=c("red", "green", "blue"), lty=1, title="Percent Female in the Table") sumout3 <- summary(out3g3) coef(sumout3) psi <- coef(sumout3)[,"Estimate"] vpsi <- vcov(sumout3) cvec1 <- c(1, 0, 0.2, 0) cvec2 <- c(1, 0, 0.5, 0) cvec3 <- c(1, 0, 0.8, 0) cvec4 <- c(0, 1, 0, 0.2) cvec5 <- c(0, 1, 0, 0.5) cvec6 <- c(0, 1, 0, 0.8) 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) cvec1 <- c(1, 0, 0, 0) cvec2 <- c(1, 1, 0, 0) cvec3 <- c(0, 0, 1, 0) cvec4 <- c(0, 0, 1, 1) 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) ################### dat3$avefemale <- ave(dat3$female, dat3$tableid) dat3$difffemale <- dat3$female - dat3$avefemale out3f1 <- lmer(sat ~ 1 + difffemale + avefemale + (1|tableid), data=dat3, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")) summary(out3f1) out3f2 <- lmer(sat ~ 1 + difffemale + avefemale + (1 + difffemale|tableid), data=dat3, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")) anova(out3f1, out3f2) out3f3 <- lmer(sat ~ 1 + difffemale + avefemale + difffemale:avefemale + (1 + difffemale|tableid), data=dat3, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")) anova(out3f2, out3f3) out3f4 <- lmer(sat ~ 1 + difffemale + avefemale + I(avefemale^2) + (1 + difffemale|tableid), data=dat3, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")) anova(out3f2, out3f4) out3f5 <- lmer(sat ~ 1 + difffemale + avefemale + I(avefemale^2) + avefemale:difffemale + I(avefemale^2):difffemale + (1 + difffemale|tableid), data=dat3, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead")) anova(out3f4, out3f5) summary(out3f4) avefemale <- seq(0, 1, 0.01) cust <- 67.6715 - (13.1943 * avefemale) + (10.2163*avefemale*avefemale) plot(avefemale, cust, type="l", xlim=c(0, 1), ylim=c(60, 70), xlab="Female Proportion in a Table", ylab="Customer Satisfaction") ############ R Square library(lme4) dat4 <-read.table("lecture6ex1.csv", sep=",", header=TRUE) dat4$aveach <- ave(dat4$ach, dat4$schoolid) dat4$diffach <- dat4$ach - dat4$aveach dat4$aveach50 <- dat4$aveach - 50 out4m1 <- lmer(efficacy ~ 1 + diffach + aveach50 + private + diffach:aveach50 + (1 + diffach|schoolid), data=dat4, REML=FALSE) summary(out4m1) library(r2mlm) r2mlm(out4m1) #### No Group Mean Centering library(lme4) library(r2mlm) dat5 <-read.table("lecture4ex2.csv", sep=",", header=TRUE) dat5$iqc <- dat5$iq - 100 dat5$int <- dat5$eesex * dat5$ersex out5 <- lmer(score ~ 1 + eesex + iqc + ersex + int + (1 + eesex|erid), data=dat5, REML=FALSE) summary(out5) sumout5 <- summary(out5) r2mlm:::r2mlm_long_manual(data=dat5, covs=c("eesex", "iqc", "ersex", "int"), random_covs=c("eesex"), gammas=coef(sumout5)[-1, "Estimate"], clusterID="erid", Tau=as.matrix(Matrix::bdiag(VarCorr(out5))), sigma2=getME(out5, "sigma")^2) #### R Squared Changed out4m2 <- lmer(efficacy ~ 1 + private + (1|schoolid), data=dat4, REML=FALSE) summary(out4m2) anova(out4m2, out4m1) r2mlm_comp(out4m2, out4m1) ############# Standardized Regression Coefficient dat1 <-read.table("lecture4ex1.csv", sep=",", header=TRUE) dat1$lengthx <- dat1$length/10 dat1$avelengthx <- ave(dat1$lengthx, dat1$groupid) dat1$difflengthx <- dat1$lengthx - dat1$avelengthx dat1$avegold <- ave(dat1$goldcolor, dat1$groupid) dat1$diffgold <- dat1$goldcolor - dat1$avegold dat1$zconsume <- scale(dat1$consume) dat1$zlength <- scale(dat1$length) dat1$avezlength <- ave(dat1$zlength, dat1$groupid) dat1$diffzlength <- dat1$zlength - dat1$avezlength volume <- dat1[!duplicated(dat1$groupid), "volume"] mvolume <- mean(volume) sdvolume <- sd(volume) dat1$zvolume <- (dat1$volume - mvolume) / sd(volume) out1z23 <- lmer(zconsume ~ 1 + diffzlength + avezlength + zvolume + diffzlength:zvolume + (1 + diffzlength|groupid), data=dat1, REML=FALSE) summary(out1z23) out1zb23 <- lmer(consume ~ 1 + diffzlength + avezlength + zvolume + diffzlength:zvolume + (1 + diffzlength|groupid), data=dat1, REML=FALSE) summary(out1zb23) out1bz23 <- lmer(zconsume ~ 1 + difflengthx + avelengthx + volume + difflengthx:volume + (1 + difflengthx|groupid), data=dat1, REML=FALSE) summary(out1bz23) out1b23 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + difflengthx:volume + (1 + difflengthx|groupid), data=dat1, REML=FALSE) summary(out1b23) r2mlm(out1b23) out1b23_1 <- lmer(consume ~ 1 + avelengthx + volume + (1|groupid), data=dat1, REML=FALSE) out1b23_2 <- lmer(consume ~ 1 + difflengthx + volume + difflengthx:volume + (1 + difflengthx|groupid), data=dat1, REML=FALSE) out1b23_3 <- lmer(consume ~ 1 + difflengthx + avelengthx + difflengthx + (1 + difflengthx|groupid), data=dat1, REML=FALSE) r2mlm_comp(out1b23_1, out1b23) r2mlm_comp(out1b23_2, out1b23) r2mlm_comp(out1b23_3, out1b23)