### ch9.r ### ### This is a provisional, undocumented r script ### ### (which may require modifications for running well) ### that was used to produce Example 9.3 in ### ### Snijders, Tom A.B., and Bosker, Roel J. ### ### Multilevel Analysis: ### ### An Introduction to Basic and Advanced Multilevel Modeling, ### ### second edition. ### ### London etc.: Sage Publishers, 2012 ### ### Script written by Tom Snijders. ### Later on, I will go back to this script and clean & document it. ### Watch out - the script moves between nlme and lme4 ### without always detaching and attaching properly. mlbook_mm <- read.table("mlbook2_mm.dat", header=TRUE) attach(mlbook_mm) # check: id2[3950:4016] names(mlbook_mm) library(lme4) # attach functions for multiple imputation source("ml2_impute_b.R") mean(mlbook_mm,na.rm=TRUE) mismlb <- is.na(as.matrix(mlbook_mm)) apply(mismlb,2,sum) # report numbers of missings manymiss <- is.na(meetings) + is.na(currmeet) + is.na(homework) hist(sch_manymiss) hist(manymiss) thedata <- cbind(mlbook_mm$langPRET,mlbook_mm$langPOST,mlbook_mm$ses, mlbook_mm$IQ_verb,mlbook_mm$IQ_perf,mlbook_mm$Minority) mean(thedata[,3],na.rm=T) # A lot of preparations follow for the functions # defined in ml2_impute_b.R. ti <- start_mvn(thedata) # basic first imputation ti <- cbind(1,ti) # add first column of 1 nvar1 <- dim(ti)[2] # number of level-1 variables tim <- thevarmean(ti,id2) # calculate group means names1 <- c("constant","langPRET","langPOST","ses","IQ_verb","IQ_perf","Minority") names2 <- paste(names1,"mean",sep=".") names1 names(ti) dim(ti) #check ti[1:30,] iqv_ses <- ti[,4]*ti[,5] miqv_mses <- (tim$V4)*(tim$V5) iqv_mmin <- ti[,5]*(tim$V7) thedata2 <- cbind(sch_manymiss,iqv_ses,miqv_mses,iqv_mmin) names3 <- c("manymiss","iqv_ses","miqv_mses","iqv_mmin") varnames <- c(names1,names2,names3) varnames impb <- list(imp_data=cbind(ti,tim,thedata2)) # for starting; names(impb[[1]]) <- varnames # also for starting; dim(impb[[1]]) thedata1 <- cbind(1,thedata) thedata.m <- thevarmean(thedata1,id2) thedata0 <- cbind(thedata1,thedata.m,thedata2) names(thedata0) <- varnames ### in the definition of impb after tim may come additional ### level-2 variables without missings mis_var <- cbind(FALSE,is.na(thedata)) # calculate missings indicators apply(mis_var,2,sum) # report numbers of missings depvar <- c(2,3,4,5,6) # give numbers of imputed variables in data frame ti ### is.binary <- c(F,F,F,F,F) table(mis_var[,2],mis_var[,3]) table(mis_var[,5],mis_var[,6]) table(mis_var[,4],mis_var[,5]) formula1 <- (langPRET ~ langPOST + ses + IQ_verb + IQ_perf + Minority + langPRET.mean + ses.mean + IQ_verb.mean + IQ_perf.mean + Minority.mean + iqv_ses + miqv_mses + iqv_mmin + (1|id2)) formula2 <- (langPOST ~ langPRET + ses + IQ_verb + IQ_perf + Minority + langPRET.mean + ses.mean + IQ_verb.mean + IQ_perf.mean + Minority.mean + iqv_ses + miqv_mses + iqv_mmin + (1|id2)) formula3 <- (ses ~ langPRET + langPOST + IQ_verb + IQ_perf + Minority + langPRET.mean + langPOST.mean + IQ_verb.mean + IQ_perf.mean + Minority.mean + (1|id2)) formula4 <- (IQ_verb ~ langPRET + langPOST + ses + IQ_perf + Minority + langPRET.mean + langPOST.mean + ses.mean + IQ_perf.mean + Minority.mean + (1|id2)) formula5 <- (IQ_perf ~ langPRET + langPOST + ses + IQ_verb + Minority + langPRET.mean + langPOST.mean + ses.mean + IQ_verb.mean + Minority.mean + (1|id2)) formulas <- c(formula1,formula2,formula3,formula4,formula5) di <- length(depvar) #formula1 <- (langPOST ~ langPRET + ses + IQ_verb + IQ_perf + langPRET.mean + ses.mean + IQ_verb.mean + IQ_perf.mean + manymiss + (1|id2)) #formula2 <- (ses ~ langPRET + langPOST + IQ_verb + IQ_perf + langPRET.mean + langPOST.mean + IQ_verb.mean + IQ_perf.mean + manymiss + (1|id2)) #formula3 <- (IQ_verb ~ langPRET + langPOST + ses + IQ_perf + langPRET.mean + langPOST.mean + ses.mean + IQ_perf.mean + manymiss + (1|id2)) #formulas <- c(formula1,formula2,formula3) formula_moi <- (langPOST ~ ses + IQ_verb + Minority + ses.mean + IQ_verb.mean + iqv_ses + miqv_mses + iqv_mmin + (1|id2)) formula_moi_av <- (langPOST ~ ses + IQ_verb + sch_ses + sch_iqv + ses:IQ_verb + (1|id2)) formula_moi_0 <- (langPOST ~ ses + IQ_verb + ses*IQ_verb + sch_ses + sch_iqv + (1|id2)) dep.name <- "langPOST :" # Model with listwise deletion of missings: fit0_moi <- lmer(formula_moi,data=thedata0) fit0_moi # Prepare design matrices XX for use in impute.mcmc.ml2 di <- length(formulas) XX <- vector("list",di) for (i in 1:di) { fiti <- lmer(formulas[[i]],data=impb[[1]]) XX[[i]] <- fiti@X } eff.names <- c('Intercept ',attr(terms(formula_moi,keep.order=TRUE),"term.labels"),"level-1 s.d.") eff.names report.model() formula_moi formula1 <- (langPOST ~ langPRET + ses + IQ_verb + IQ_perf + langPRET.mean + ses.mean + IQ_verb.mean + IQ_perf.mean + manymiss + (1|id2)) formula2 <- (ses ~ langPRET + langPOST + IQ_verb + IQ_perf + langPRET.mean + langPOST.mean + IQ_verb.mean + IQ_perf.mean + manymiss + (1|id2)) formula3 <- (IQ_verb ~ langPRET + langPOST + ses + IQ_perf + langPRET.mean + langPOST.mean + ses.mean + IQ_perf.mean + manymiss + (1|id2)) formulas <- c(formula1,formula2,formula3) formula_moi <- (langPOST ~ ses + IQ_verb + ses.mean + IQ_verb.mean + ses:IQ_verb + (1|id2)) report.model() formula1a <- (langPRET ~ langPOST + ses + IQ_verb + IQ_perf + langPOST.mean + ses.mean + IQ_verb.mean + IQ_perf.mean + ses:IQ_verb + IQ_verb:IQ_perf + langPOST:IQ_verb + ses:IQ_perf + (1|id2)) formula2a <- (langPOST ~ langPRET + ses + IQ_verb + IQ_perf + langPRET.mean + ses.mean + IQ_verb.mean + IQ_perf.mean + ses:IQ_verb + IQ_verb:IQ_perf + langPRET:IQ_verb + ses:IQ_perf + (1|id2)) formula3a <- (ses ~ langPRET + langPOST + IQ_verb + IQ_perf + langPRET.mean + langPOST.mean + IQ_verb.mean + IQ_perf.mean + langPOST:IQ_verb + IQ_verb:IQ_perf + langPRET:IQ_verb + langPOST:IQ_perf + (1|id2)) formula4a <- (IQ_verb ~ langPRET + langPOST + ses + IQ_perf + langPRET.mean + langPOST.mean + ses.mean + IQ_perf.mean + langPOST:ses + (1|id2)) formula5a <- (IQ_perf ~ langPRET + langPOST + ses + IQ_verb + langPRET.mean + langPOST.mean + ses.mean + IQ_verb.mean + langPOST:ses + (1|id2)) formulas <- c(formula1a,formula2a,formula3a,formula4a,formula5a) di <- length(formulas) eff.names <- c('Intercept ',attr(terms(formula_moi,keep.order=TRUE),"term.labels"),"level-1 s.d.") eff.names report.model() # Prepare design matrices XX for use in impute.mcmc.ml2 di <- length(formulas) XX <- vector("list",di) for (i in 1:di) {fiti <- lmer(formulas[[i]],data=impb[[1]]) XX[[i]] <- fiti@X } formula1b <- (langPRET ~ langPOST + ses + IQ_verb + IQ_perf + langPOST.mean + ses.mean + IQ_verb.mean + IQ_perf.mean + ses:IQ_verb + IQ_verb:IQ_perf + langPOST:IQ_verb + ses:IQ_perf + manymiss + (1|id2)) formula2b <- (langPOST ~ langPRET + ses + IQ_verb + IQ_perf + langPRET.mean + ses.mean + IQ_verb.mean + IQ_perf.mean + ses:IQ_verb + IQ_verb:IQ_perf + langPRET:IQ_verb + ses:IQ_perf + manymiss + (1|id2)) formula3b <- (ses ~ langPRET + langPOST + IQ_verb + IQ_perf + langPRET.mean + langPOST.mean + IQ_verb.mean + IQ_perf.mean + langPOST:IQ_verb + IQ_verb:IQ_perf + langPRET:IQ_verb + langPOST:IQ_perf + manymiss +(1|id2)) formula4b <- (IQ_verb ~ langPRET + langPOST + ses + IQ_perf + langPRET.mean + langPOST.mean + ses.mean + IQ_perf.mean + langPOST:ses + manymiss + (1|id2)) formula5b <- (IQ_perf ~ langPRET + langPOST + ses + IQ_verb + langPRET.mean + langPOST.mean + ses.mean + IQ_verb.mean + langPOST:ses + manymiss + (1|id2)) formulas <- c(formula1b,formula2b,formula3b,formula4b,formula5b) di <- length(formulas) eff.names <- c('Intercept ',attr(terms(formula_moi,keep.order=TRUE),"term.labels"),"level-1 s.d.") eff.names report.model() # Prepare design matrices XX for use in impute.mcmc.ml2 di <- length(formulas) XX <- vector("list",di) for (i in 1:di) {fiti <- lmer(formulas[[i]],data=impb[[1]]) XX[[i]] <- fiti@X } depvar <- c(3,4,5) # give numbers of imputed variables in data frame ti ### formula1 <- (langPOST ~ ses + IQ_verb + ses:IQ_verb + ses.mean + IQ_verb.mean + (1|id2)) formula2 <- (ses ~ langPOST + IQ_verb + langPOST.mean + IQ_verb.mean + (1|id2)) formula3 <- (IQ_verb ~ langPOST + ses + langPOST.mean + ses.mean + (1|id2)) formulas <- c(formula1,formula2,formula3) # Prepare design matrices XX for use in impute.mcmc.ml2 di <- length(formulas) XX <- vector("list",di) for (i in 1:di) {fiti <- lmer(formulas[[i]],data=impb[[1]]) XX[[i]] <- fiti@X } # following are names for us in output in print.multi.imp, # refers to model of interest. ### : eff.names <- c('Intercept :','IQ_verb :','IQ_v mean :') # all same length, no trailing blanks dep.name <- "langPOST" #model1 <- lme(langPOST ~ IQ_verb + IQ_verb.mean,data=impb[[1]], # random=~1|id2,method="ML") #sm1 <- summary(model1) ### : mip1 <- multi.impute(100,50,impb) # first is inner loop, second is number imputations report.model() print.multi.imp(mip1) # results for model of interest list.multi.fit(mip1) # results for model of interest fit0 <- lmer(langPOST ~ IQ_verb + sch_iqv + (1|schoolnr), control=lmeControl(maxIter=100,msMaxIter=100)) print(fit0) print.multi.imp(mip1) length(mip1[[1]]) xyplot(ti[,3] ~ ti[,5], groups=mis_var[,3]) xyplot(mip1[[1]][[50]][,3] ~ mip1[[1]][[50]][,5], groups=mis_var[,3]) # blue = observed, red = imputed ================================================== # Now for checking something else. hist(langPRET) use1 <- (langPRET >= 30) table(use1) langPOST[!use1] <- NA thedata.use1 <- cbind(mlbook_mm$langPRET,langPOST,mlbook_mm$ses, mlbook_mm$IQ_verb,mlbook_mm$IQ_perf) # thedata <- cbind(mlbook_mm$langPRET,mlbook_mm$langPOST,mlbook_mm$ses, mlbook_mm$IQ_verb,mlbook_mm$IQ_perf) thedata2 <- cbind(sch_manymiss) names1 <- c("constant","langPRET","langPOST","ses","IQ_verb","IQ_perf") names2 <- paste(names1,"mean",sep=".") names3 <- c("manymiss") varnames <- c(names1,names2,names3) varnames thedata.use1 <- cbind(1,thedata.use1) # add first column of 1 thedata.use1.m <- thevarmean(thedata.use1,id2) # calculate group means alldata.use1 <- cbind(thedata.use1,thedata.use1.m,thedata2) names(alldata.use1) <- varnames ti.use1 <- start_mvn(thedata.use1) # basic first imputation ti.use1 <- cbind(1,ti.use1) # add first column of 1 nvar1 <- dim(ti)[2] # number of level-1 variables tim.use1 <- thevarmean(ti.use1,id2) # calculate group means impb.use1 <- list(imp_data=cbind(ti.use1,tim.use1,thedata2)) # for starting; names(impb.use1[[1]]) <- varnames mis_var <- cbind(FALSE,is.na(thedata.use1)) # calculate missings indicators apply(mis_var,2,sum) # report numbers of missings fit0_moi.use1 <- lmer(formula_moi,data=data.frame(alldata.use1)) fit0_moi.use1 fit1_moi.use1 <- lmer(formula_moi,data=data.frame(impb.use1[[1]])) fit1_moi.use1 mip.use1 <- multi.impute(20,20,impb.use1) # first is inner loop, second is number imputations print.multi.imp(mip.use1) # results for model of interest list.multi.fit(mip.use1) # results for model of interest # Prepare design matrices XX for use in impute.mcmc.ml2 di <- length(formulas) XX <- vector("list",di) for (i in 1:di) { fiti <- lmer(formulas[[i]],data=impb.use1[[1]]) XX[[i]] <- fiti@X } mip3 <- multi.mc.impute(10,10,100,50,impb) # nmc0,nmc,k,m, mip.use2 <- multi.mc.impute(10,10,100,20,impb.use1) # nmc0,nmc,k,m, mip4 <- multi.impute(100,50,impb.use1) mip4.use1 <- mip4 rm(mip4) print.multi.imp(mip1) # results for model of interest print.multi2.mc.imp(mip3) # results for model of interest mip4 <- multi.impute(100,50,impb) formu.f <- "ses + IQ_verb + ses.mean + IQ_verb.mean + ses:IQ_verb" formu.f <- formula(formu.f) # attach nlme library(nlme) # attach utilities source("nlme_utils.R") imp.nlme.1 <- multi.fit.nlme(mip1, langPOST ~ ses + IQ_verb + ses.mean + IQ_verb.mean + ses:IQ_verb , ~1|id2) print.multi.imp(imp.nlme.1) imp.nlme.4 <- multi.fit.nlme(mip4, langPOST ~ ses + IQ_verb + ses.mean + IQ_verb.mean + ses:IQ_verb , ~1|id2) print.multi.imp(imp.nlme.4) print.multi.imp(mip4) save(mip4,file="mip4.Rdata") eff.names <- c('Intercept ',attr(terms(formula_moi,keep.order=TRUE),"term.labels"),"level-1 s.d.") ===================================== # iqv_ses <- ti[4]*ti[5] # miqv_mses <- (tim$V4)*(tim$V5) # iqv_mmin <- ti[5]*(tim$V7) # thedata2 <- cbind(sch_manymiss,iqv_ses,miqv_mses,iqv_mmin) varnames extraCalculations <- function(imput){ # This can be redefined, and serves for last calculations, e.g., # calculation of interactions or dummy variables from basic imputed variables. # If extraCalculations is redefined, then also imput_b must be read again. # Example: tt <- imput$imp_data tt$iqv_ses <- (tt$IQ_verb)*(tt$ses) tt$miqv_mses <- (tt$IQ_verb.mean)*(tt$ses.mean) tt$iqv_mmin <- (tt$IQ_verb)*(tt$Minority.mean) return(list(imp_data=tt, imp_coefs=imput$imp_coefs, imp_sters=imput$imp_sters)) } fit0_moi <- lmer(formula_moi,data=thedata0) fit1_moi <- lmer(formula_moi,data=impb[[1]]) fit1_moi fit0_moi mip5 <- multi.impute(100,50,impb) print.multi.imp(mip5) save(mip5,file="mip5.Rdata") report.model() formula1 <- (langPRET ~ langPOST + ses + IQ_verb + IQ_perf + Minority + langPRET.mean + ses.mean + IQ_verb.mean + IQ_perf.mean + Minority.mean + iqv_ses + miqv_mses + iqv_mmin + (1|id2)) formula2 <- (langPOST ~ langPRET + ses + IQ_verb + IQ_perf + Minority + langPRET.mean + ses.mean + IQ_verb.mean + IQ_perf.mean + Minority.mean + iqv_ses + miqv_mses + iqv_mmin + (1|id2)) formula3 <- (ses ~ langPRET + langPOST + IQ_verb + IQ_perf + Minority + langPRET.mean + langPOST.mean + IQ_verb.mean + IQ_perf.mean + Minority.mean + (1|id2)) formula4 <- (IQ_verb ~ langPRET + langPOST + ses + IQ_perf + Minority + langPRET.mean + langPOST.mean + ses.mean + IQ_perf.mean + Minority.mean + (1|id2)) formula5 <- (IQ_perf ~ langPRET + langPOST + ses + IQ_verb + Minority + langPRET.mean + langPOST.mean + ses.mean + IQ_verb.mean + Minority.mean + (1|id2)) formulas <- c(formula1,formula2,formula3,formula4,formula5) di <- length(depvar) #formula1 <- (langPOST ~ langPRET + ses + IQ_verb + IQ_perf + langPRET.mean + ses.mean + IQ_verb.mean + IQ_perf.mean + manymiss + (1|id2)) #formula2 <- (ses ~ langPRET + langPOST + IQ_verb + IQ_perf + langPRET.mean + langPOST.mean + IQ_verb.mean + IQ_perf.mean + manymiss + (1|id2)) #formula3 <- (IQ_verb ~ langPRET + langPOST + ses + IQ_perf + langPRET.mean + langPOST.mean + ses.mean + IQ_perf.mean + manymiss + (1|id2)) #formulas <- c(formula1,formula2,formula3) formula_moi <- (langPOST ~ ses + IQ_verb + Minority + ses.mean + IQ_verb.mean + iqv_ses + miqv_mses + iqv_mmin + (1|id2)) eff.names <- c('Intercept ',attr(terms(formula_moi,keep.order=TRUE),"term.labels"),"level-1 s.d.") eff.names print.multi.imp(mip6) imp.nlme.6 <- multi.fit.nlme(mip6, langPOST ~ ses + IQ_verb + ses.mean + IQ_verb.mean + ses:IQ_verb + ses.mean:IQ_verb.mean, ~1|id2) print.multi.imp(imp.nlme.6) model6 <- lme(langPOST ~ses + IQ_verb + sch_ses + sch_iqv + ses:IQ_verb + sch_ses:sch_iqv,data=mlbook_mm, random=~1|id2,na.action=na.omit,method="ML") summary(model6) GetRandomVars1(model6) GetRandomSD1(model6) print.multi.imp(imp.nlme.6)