Paris Lecture 1 code

gistfile1.txt

### R code from vignette source '01_ADAIntroLecture.Rnw'

###################################################
### code chunk number 1: 01_ADAIntroLecture.Rnw:134-141
###################################################
## load data:
data<-read.table("~/Git/Statistics-lecture-notes-Potsdam/AdvancedDataAnalysis/data/gibsonwu2012data.txt",header=T)
## take reciprocal rt to normalize residuals:
data$rrt<- -1000/data$rt
## define predictor x:
data$x <- ifelse(
  data$type%in%c("subj-ext"),-0.5,0.5)


###################################################
### code chunk number 2: 01_ADAIntroLecture.Rnw:144-151
###################################################
headnoun<-subset(data,region=="headnoun")
## no. of subjects:
#length(unique(headnoun$subj))
## no. of items:
#length(unique(headnoun$item))
## no. of rows in data frame:
#dim(headnoun)[1]


###################################################
### code chunk number 3: 01_ADAIntroLecture.Rnw:177-178
###################################################
head(data[,c(1,2,3,4,7,10,11)])


###################################################
### code chunk number 4: 01_ADAIntroLecture.Rnw:186-187
###################################################
library(lme4)


###################################################
### code chunk number 5: 01_ADAIntroLecture.Rnw:192-194
###################################################
m1 <- lmer(rrt~x+(1+x|subj)+(1+x|item),
          subset(data,region=="headnoun"))


###################################################
### code chunk number 6: 01_ADAIntroLecture.Rnw:214-215
###################################################
summary(m1)


###################################################
### code chunk number 7: 01_ADAIntroLecture.Rnw:225-229
###################################################
m1<- lmer(rrt~x+(1+x|subj)+(1+x|item),
          headnoun)
m1a<- lmer(rrt~x+(1|subj)+(1|item),
          headnoun)


###################################################
### code chunk number 8: 01_ADAIntroLecture.Rnw:238-239
###################################################
anova(m1,m1a)


###################################################
### code chunk number 9: 01_ADAIntroLecture.Rnw:257-315
###################################################
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 10: 01_ADAIntroLecture.Rnw:325-334
###################################################
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 11: 01_ADAIntroLecture.Rnw:344-345
###################################################
nsim<-100


###################################################
### code chunk number 12: 01_ADAIntroLecture.Rnw:359-369
###################################################
library(lme4)
subjcorr<-rep(NA,nsim)
itemcorr<-rep(NA,nsim)

for(i in 1:nsim){
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 13: 01_ADAIntroLecture.Rnw:379-386
###################################################
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 14: 01_ADAIntroLecture.Rnw:396-406
###################################################
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 15: 01_ADAIntroLecture.Rnw:415-422
###################################################
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 16: 01_ADAIntroLecture.Rnw:565-566
###################################################
library(lme4)


###################################################
### code chunk number 17: 01_ADAIntroLecture.Rnw:571-573
###################################################
m1 <- lmer(rrt~x+(1+x|subj)+(1+x|item),
          headnoun)


###################################################
### code chunk number 18: 01_ADAIntroLecture.Rnw:581-582
###################################################
summary(m1)


###################################################
### code chunk number 19: 01_ADAIntroLecture.Rnw:795-856
###################################################
se <- function(x)
      {
        y <- x[!is.na(x)] # remove the missing values, if any
        sqrt(var(as.vector(y))/length(y))
}


ci <- function (scores){
m <- mean(scores,na.rm=TRUE)
stderr <- se(scores)
len <- length(scores)
upper <- m + qt(.975, df=len-1) * stderr 
lower <- m + qt(.025, df=len-1) * stderr 
return(data.frame(lower=lower,upper=upper))
}

## sample repeatedly:
lower <- rep(NA,100)
upper <- rep(NA,100)

for(i in 1:100){ 
  sample <- rnorm(100,mean=60,sd=4)
  lower[i] <- ci(sample)$lower
  upper[i] <- ci(sample)$upper
}
  
cis <- cbind(lower,upper)

store <- rep(NA,100)

pop.mean<-60
pop.sd<-4

for(i in 1:100){ 
  sample <- rnorm(100,mean=pop.mean,sd=pop.sd)
  lower[i] <- ci(sample)$lower
  upper[i] <- ci(sample)$upper
  if(lower[i]<pop.mean & upper[i]>pop.mean){
    store[i] <- TRUE} else {
      store[i] <- FALSE}
}

## need this for the plot below:
cis <- cbind(lower,upper)


## convert store to factor:
store<-factor(store)

main.title<-"95% CIs in 100 repeated samples"

line.width<-ifelse(store==FALSE,2,1)
cis<-cbind(cis,line.width)
x<-0:100
y<-seq(55,65,by=1/10)
plot(x,y,type="n",xlab="i-th repeated sample",ylab="Scores",main=main.title)
abline(60,0,lwd=2)
x0<-x
x1<-x
arrows(x0,y0=cis[,1],
       x1,y1=cis[,2],length=0,lwd=cis[,3])