Maximal models in linear mixed models

recoveringcorrelationsV2.R

### R code from vignette source 'recoveringcorrelationsV2.Rnw'

###################################################
### code chunk number 1: recoveringcorrelationsV2.Rnw:98-156
###################################################
new.df <- function(cond1.rt=487, effect.size=123, 
                   sdev=544,
                   sdev.int.subj=160, sdev.slp.subj=195, 
                   rho.u=0.6,
                   nsubj=37,
                   sdev.int.items=154, sdev.slp.items=142,
                   rho.w=0.6,
                   nitems=15) {
  library(MASS)
  
  ncond <- 2
  
  subj <- rep(1:nsubj, each=nitems*ncond)
  item <- rep(1:nitems, nsubj, each=ncond)
  
  cond <- rep(0:1, nsubj*nitems)
  err  <- rnorm(nsubj*nitems*ncond, 0, sdev)
  d    <- data.frame(subj=subj, item=item, 
                     cond=cond+1, err=err)
  
  Sigma.u<-matrix(c(sdev.int.subj^2,
                    rho.u*sdev.int.subj*sdev.slp.subj,
                    rho.u*sdev.int.subj*sdev.slp.subj,
                    sdev.slp.subj^2),nrow=2)
  
  Sigma.w<-matrix(c(sdev.int.items^2,
                    rho.u*sdev.int.items*sdev.slp.items,
                    rho.u*sdev.int.items*sdev.slp.items,
                    sdev.slp.items^2),nrow=2)
  
  # Adding random intercepts and slopes for subjects:
  ## first col. has adjustment for intercept, 
  ## secdon col. has adjustment for slope
  subj.rand.effs<-mvrnorm(n=nsubj,rep(0,ncond),Sigma.u)
  
  item.rand.effs<-mvrnorm(n=nitems,rep(0,ncond),Sigma.w)
  
  re.int.subj <- subj.rand.effs[,1]
  d$re.int.subj <- rep(re.int.subj, each=nitems*ncond)
  re.slp.subj   <- subj.rand.effs[,2]
  
  d$re.slp.subj <- rep(re.slp.subj, 
                       each=nitems*ncond) * (cond - 0.5)
  
  re.int.item <- item.rand.effs[,1]
  d$re.int.item <- rep(re.int.item, nsubj, each=ncond)
  re.slp.item <- item.rand.effs[,2]
  d$re.slp.item <- rep(re.slp.item, nsubj, 
                       each=ncond) * (cond - 0.5)
  
  d$rt <- (cond1.rt + cond*effect.size
           + d$re.int.subj + d$re.slp.subj
           + d$re.int.item + d$re.slp.item
           + d$err)
  
  return(list(d,cor(re.int.subj,re.slp.subj),
              cor(re.int.item,re.slp.item)))
}


###################################################
### code chunk number 2: recoveringcorrelationsV2.Rnw:161-170
###################################################
gendata<-function(subjects=37,items=15){
  dat<-new.df(nsubj=subjects,nitems=items,
              rho.u=0.6,rho.w=0.6)
  dat <- dat[[1]]
  dat<-dat[,c(1,2,3,9)]
  dat$x<-ifelse(dat$cond==1,-0.5,0.5)
  
return(dat)
}


###################################################
### code chunk number 3: recoveringcorrelationsV2.Rnw:175-176
###################################################
nsim<-100


###################################################
### code chunk number 4: recoveringcorrelationsV2.Rnw:184-204
###################################################
library(lme4)
nsim<-100
subjcorr<-rep(NA,nsim)
itemcorr<-rep(NA,nsim)

fixef_max<-matrix(rep(NA,2*nsim),ncol=2)
fixef_min<-matrix(rep(NA,2*nsim),ncol=2)

for(i in 1:nsim){
#print(i)
dat<-gendata()  
m_max<-lmer(rt~x+(1+x|subj)+(1+x|item),dat)
m_min<-lmer(rt~x+(1|subj)+(1|item),dat)

fixef_max[i,]<-summary(m_max)$coef[2,1:2]
fixef_min[i,]<-summary(m_min)$coef[2,1:2]

subjcorr[i]<-attr(VarCorr(m_max)$subj,"correlation")[1,2]
itemcorr[i]<-attr(VarCorr(m_max)$item,"correlation")[1,2]
}


###################################################
### code chunk number 5: recoveringcorrelationsV2.Rnw:209-216
###################################################
op<-par(mfrow=c(1,2),pty="s")
hist(subjcorr,freq=FALSE,xlab=expression(hat(rho)[u]),
     main="Distribution of subj. corr.")
abline(v=0.6,lwd=3)
hist(itemcorr,freq=FALSE,xlab=expression(hat(rho)[w]),
     main="Distribution of item corr.")
abline(v=0.6,lwd=3)


###################################################
### code chunk number 6: recoveringcorrelationsV2.Rnw:221-233
###################################################
lower_max<-fixef_max[,1]-2*fixef_max[,2]
upper_max<-fixef_max[,1]+2*fixef_max[,2]

lower_min<-fixef_min[,1]-2*fixef_min[,2]
upper_min<-fixef_min[,1]+2*fixef_min[,2]

plot(1:100,fixef_max[,1],ylim=c(-100,500),xlab="repeated samples",ylab="rt")
arrows(1:100,lower_max,1:100,upper_max,length=0)
abline(h=0)

points(1:100+0.5,fixef_min[,1],ylim=c(-100,500),bg="red",pch=21)
arrows(1:100+0.5,lower_min,1:100+0.5,upper_min,length=0,col="red")


###################################################
### code chunk number 7: recoveringcorrelationsV2.Rnw:236-243
###################################################
## incorrect failure to reject null
## under maximal model:
mean(lower_max<=0)

## incorrect failure to reject under 
## minimal model:
mean(lower_min<=0)


###################################################
### code chunk number 8: recoveringcorrelationsV2.Rnw:251-269
###################################################
nsim<-100
subjcorr<-rep(NA,nsim)
itemcorr<-rep(NA,nsim)

fixef_max<-matrix(rep(NA,2*nsim),ncol=2)
fixef_min<-matrix(rep(NA,2*nsim),ncol=2)

for(i in 1:nsim){
dat<-gendata(subjects=100,items=100) 
m_max<-lmer(rt~x+(1+x|subj)+(1+x|item),dat)
m_min<-lmer(rt~x+(1|subj)+(1|item),dat)

fixef_max[i,]<-summary(m_max)$coef[2,1:2]
fixef_min[i,]<-summary(m_min)$coef[2,1:2]

subjcorr[i]<-attr(VarCorr(m_max)$subj,"correlation")[1,2]
itemcorr[i]<-attr(VarCorr(m_max)$item,"correlation")[1,2]
}


###################################################
### code chunk number 9: recoveringcorrelationsV2.Rnw:274-281
###################################################
op<-par(mfrow=c(1,2),pty="s")
hist(subjcorr,freq=FALSE,xlab=expression(hat(rho)[u]),
     main="Distribution of subj. corr.")
abline(v=0.6,lwd=3)
hist(itemcorr,freq=FALSE,xlab=expression(hat(rho)[w]),
     main="Distribution of item corr.")
abline(v=0.6,lwd=3)


###################################################
### code chunk number 10: recoveringcorrelationsV2.Rnw:287-299
###################################################
lower_max<-fixef_max[,1]-2*fixef_max[,2]
upper_max<-fixef_max[,1]+2*fixef_max[,2]

lower_min<-fixef_min[,1]-2*fixef_min[,2]
upper_min<-fixef_min[,1]+2*fixef_min[,2]

plot(1:100,fixef_max[,1],ylim=c(-100,500),xlab="repeated samples",ylab="rt")
arrows(1:100,lower_max,1:100,upper_max,length=0)
abline(h=0)

points(1:100+0.5,fixef_min[,1],ylim=c(-100,500),bg="red",pch=21)
arrows(1:100+0.5,lower_min,1:100+0.5,upper_min,length=0,col="red")


###################################################
### code chunk number 11: recoveringcorrelationsV2.Rnw:302-309
###################################################
## incorrect failure to reject null
## under maximal model:
mean(lower_max<=0)

## incorrect failure to reject under 
## minimal model:
mean(lower_min<=0)

The function new.df doesn't work. I got the following error message:

Error in mvrnorm(n = nsubj, rep(0, ncond), Sigma.u) (from #24) : 
  incompatible arguments

I'm using R version 3.0.2 (2013-09-25) from Ubuntu 14.04.