library(foreign) mlbook <- read.dta("mlbook1.dta") mlbook[1:10,] names(mlbook) attach(mlbook) schdata <- aggregate(mlbook, by = list(mlbook$schoolnr), mean) schoolsb <- schdata[,c('schoolnr', 'langpost', 'iq_verb', 'ses', 'mixedgra', 'groupsiz')] names(schoolsb) <- c("schoolnr", "meanlp", "meaniq", "meanses", "combi", "grs") schoolsz <- table(mlbook$schoolnr) mlbook <- merge(mlbook, schoolsb[, 1:4], by='schoolnr') attach(mlbook) ##total students nrow(mlbook) ##total schools nrow(schdata) library(lattice) library(nlme) samp25 <- sample(unique(schoolnr), 25) sample25 <- groupedData(langpost ~ iq_verb | schoolnr, data = mlbook[is.element(schoolnr, samp25), ]) xyplot(langpost ~ iq_verb|schoolnr, data=sample25, main="Exploring 25 random schools", panel=function(x,y){ panel.xyplot(x, y) panel.loess(x, y, span=1) panel.lmline(x, y, lty=2) } ) ##note plot(sample25) shows you the data in a similar plot summary(model0 <- lme(langpost ~ 1, random =~1|schoolnr)) (vc <- nlme::VarCorr(model0)) ## finding the se of the variance estimate ## intervals fn uses component apVar and attribute of this, Pars ## bit of delta functioning too selog = sqrt(diag(model0$apVar)) senolog = exp(attr(model0$apVar,'Pars')) * selog sevarnolog = 2 * exp(attr(model0$apVar,'Pars')) * senolog (model1 <- lme(langpost ~ iq_verb, random=~1|schoolnr)) (vc <- VarCorr(model1)) detach("mlbook") mlbook$iqdev <- mlbook$iq_verb - mean(mlbook$iq_verb) attach(mlbook) summary(model2 <- lme(langpost ~ iqdev , random =~1|schoolnr)) (vc <- VarCorr(model2)) fixef(model1) ranef(model1) summary(model2n <- lme(langpost ~ iqdev + meaniq + ses + sex, random=~1|schoolnr)) qqnorm(model2n, ~ranef(.)) ## unstandardised ebe2 <- ranef(model2n) bv2 <- exp(attr(model2n$apVar,'Pars'))[1] sebe2 <- ebe2/sqrt(bv2) qqmath(~sebe2) summary(model3 <- lme(langpost ~iqdev + meaniq, random= ~iqdev|schoolnr)) (vc <- VarCorr(model3)) detach('mlbook') mean(mlbook$meaniq) mean(mlbook$groupsiz) mlbook$meaniq <- mlbook$meaniq - mean(mlbook$meaniq) mlbook$groupsiz <- mlbook$groupsiz - mean(mlbook$groupsiz) attach(mlbook) summary(model4 <- lme(langpost ~iqdev + meaniq + groupsiz + iqdev*groupsiz, random=~iqdev|schoolnr, method='ML')) -2*logLik(model1) ##deviance logLik(model1) p(model5 <- lmer(langpost ~ iqdev + meaniq + groupsiz + iqdev*groupsiz + (1|schoolnr), REML=FALSE)) summary(model5 <- lme(langpost ~ iqdev + meaniq + groupsiz + iqdev*groupsiz, random = ~1|schoolnr, method='ML')) anova(model5,model4) summary(model1n <- lme(langpost ~iqdev + ses + meaniq + sex, random = ~ 1 |schoolnr)) summary(model2n <- lme(langpost ~ iqdev + ses + meaniq + sex, random = ~ iqdev |schoolnr)) anova(model1n, model2n) qqnorm(model2n, ~ranef(.)) qqnorm(model2n, ~resid(.)) olse2list <- lmList(langpost ~ iqdev + ses + meaniq + sex | schoolnr, data=mlbook) ## lmList seems to demand data=! compm2 <- compareFits(coef(olse2list), coef(model2n)) olse2list[1:3] plot(compm2, mark=fixef(model2n)) ##data = mlbook[is.element(schoolnr, samp25), ] resi2 <- residuals(olse2list) plot(iq_verb,langpost) plot(jitter(iq_verb),jitter(langpost)) lines(lowess(iq_verb,langpost, f=.2), lwd=3) plot(jitter(iqdev),resi2, lines(lowess(iqdev,resi2, f=.2), lwd=3))