Recovering correlations

gistfile1.txt

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

###################################################
### code chunk number 1: recoveringcorrelations.Rnw:93-151
###################################################
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: recoveringcorrelations.Rnw:156-165
###################################################
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: recoveringcorrelations.Rnw:170-171
###################################################
nsim<-100


###################################################
### code chunk number 4: recoveringcorrelations.Rnw:179-190
###################################################
library(lme4)
subjcorr<-rep(NA,nsim)
itemcorr<-rep(NA,nsim)

for(i in 1:nsim){
#print(i)
dat<-gendata()  
m3<-lmer(rt~x+(1+x|subj)+(1+x|item),dat)
subjcorr[i]<-attr(VarCorr(m3)$subj,"correlation")[1,2]
itemcorr[i]<-attr(VarCorr(m3)$item,"correlation")[1,2]
}


###################################################
### code chunk number 5: recoveringcorrelations.Rnw:193-200
###################################################
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: recoveringcorrelations.Rnw:205-215
###################################################
subjcorr<-rep(NA,nsim)
itemcorr<-rep(NA,nsim)

for(i in 1:nsim){
#print(i)
dat<-gendata(subjects=50,items=30)  
m3<-lmer(rt~x+(1+x|subj)+(1+x|item),dat)
subjcorr[i]<-attr(VarCorr(m3)$subj,"correlation")[1,2]
itemcorr[i]<-attr(VarCorr(m3)$item,"correlation")[1,2]
}


###################################################
### code chunk number 7: recoveringcorrelations.Rnw:218-225
###################################################
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)