# install.packages(lmerTest) # install.packages(MASS) nLevel0 = 30 nLevel1 = 10 ### fixed effects b0 = 20 b1 = 10 ### standard deviation random effects varRandomInterceptLevel1 = 4 varRandomSlopeLevel1 = 3 ### standard deviation error sdError = 1 ### covariation between random effects covRandom = 3 ### calculate confidence intervals? yes = 1, no = 0 confInt = 1 ################################################################# ################################################################# ################################################################# expand = .5 yexpand = 1 xexpand = 1 perfectRepresentative = 0 ################################################################# ################################################################# ################################################################# # load libraries library(lmerTest) library(MASS) ## perfectRepresentative = ifelse(perfectRepresentative==1,TRUE,FALSE) # predictor from standard normal distibution pred = rnorm(nLevel0 * nLevel1,0,1) ######### ######### ######### Sigma = cbind(c(varRandomInterceptLevel1,covRandom),c(covRandom,varRandomSlopeLevel1)) if(varRandomInterceptLevel1 != 0 & varRandomSlopeLevel1 != 0){ randomParam = mvrnorm(n = nLevel1, c(0,0), Sigma, empirical =perfectRepresentative) }else{ randomParam = cbind(rnorm(nLevel1,0,varRandomInterceptLevel1),rnorm(nLevel1,0,varRandomSlopeLevel1)) } randomInterceptLevel1 = randomParam[,1] randomInterceptLevel1 = rep(randomInterceptLevel1,each=nLevel0) RandomSlopeLevel1 = randomParam[,2] RandomSlopeLevel1 = rep(RandomSlopeLevel1,each=nLevel0) subject = rep(1:nLevel1,each=nLevel0) crit = (b0 + randomInterceptLevel1) + (b1 + RandomSlopeLevel1) * pred errorLevel1 = rnorm(nLevel0*nLevel1,0,sdError) crit = crit + errorLevel1 ## fir multilevel model modMLR = lmer(crit ~ pred + (1+pred|subject),REML=0) ## plotting colPart = rep(rainbow(nLevel1),each=nLevel0) plot(pred,crit,col=colPart,pch=16,cex=1.2,xlab="Predictor", ylab = "Criterion",ylim = c(min(crit)-1,max(crit)+1) * yexpand, xlim = c(min(pred)-1,max(pred)+1) * xexpand) fixedEffects = unlist(fixef(modMLR)) randomEffects = ranef(modMLR) randomEffects = matrix(unlist(randomEffects),ncol=2) colUnique = unique(colPart) for(loop in 1 : nLevel1){ f = function(x){fixedEffects[1] + randomEffects[loop,1] + (fixedEffects[2] + randomEffects[loop,2] ) * x} curve(f,-1000,1000,add=T,col=colUnique[loop],lwd=2) } abline(h=0,col="black",lty=3) abline(v=0,col="black",lty=3) ############ output print("Summary Model") summary(modMLR) print("True model") b0 b1 sqrt(varRandomInterceptLevel1) sqrt(varRandomSlopeLevel1) sdError print("Fixed Effects") fixedEffects print("Random Effects") randomEffects print("Variance Covariance Matrix") VarCorr(modMLR) if(confInt == 1){ print("Confidence Intervals for Fixed and Random Effects based on Likelihood Ratio Test]") print("") print("") confint(modMLR) }