Paris lecture 3 code

ParisLecture3.R

### R code from vignette source '03_ADALecture3.Rnw'

###################################################
### code chunk number 1: 03_ADALecture3.Rnw:85-92
###################################################
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)
headnoun<-subset(data,region=="headnoun")


###################################################
### code chunk number 2: 03_ADALecture3.Rnw:96-106
###################################################

## data for JAGS and Stan:
headnoun.dat <- list(subj=
                sort(as.integer(factor(headnoun$subj) )),
                     item=sort(as.integer(factor(headnoun$item))),
                    rrt = headnoun$rrt,
                    x = headnoun$x,
                    I = nrow(headnoun),
                    J =length( unique(headnoun$subj) ),
                    K =length( unique(headnoun$item)))


###################################################
### code chunk number 3: 03_ADALecture3.Rnw:299-362
###################################################
cat("
data
{
    zero[1] <- 0
    zero[2] <- 0
}
    model
{
    # Intercept and slope for each subj
    for( j in 1:J )
{
    u[j,1:2] ~ dmnorm(zero,Omega.u)
}    
    # Intercept and slope for each item
    for( k in 1:K )
{
    w[k,1:2] ~ dmnorm(zero,Omega.w)
}
    # Define model for each observational unit
    for( i in 1:I )
{
    mu[i] <- ( beta[1] + u[subj[i],1] + w[item[i],1]) +
    ( beta[2] + u[subj[i],2] +  w[item[i],2]) * ( x[i] ) 
    rrt[i] ~ dnorm( mu[i], tau.e )
}
    # Priors:   
    # Fixed intercept and slope (uninformative)
    beta[1] ~ dnorm(0.0,1.0E-5)
    beta[2] ~ dnorm(0.0,1.0E-5)
    # Residual variance
    tau.e <- pow(sigma.e,-2)
    sigma.e  ~ dunif(0,100)    
    # Define prior for the variance-covariance matrix of the random effects for subjects
    ## precision:
    Omega.u  ~ dwish( R.u, 2 )
    ## R matrix:
    R.u[1,1] <- pow(sigma.a,2)
    R.u[2,2] <- pow(sigma.b,2)
    R.u[1,2] <- rho.u*sigma.a*sigma.b
    R.u[2,1] <- R.u[1,2]
    ## Vcov matrix:
    Sigma.u <- inverse(Omega.u)
    ## priors for var int. var slopes
    sigma.a ~ dunif(0,10)
    sigma.b ~ dunif(0,10)  
    ## prior for correlation:
    rho.u ~ dunif(-1,1)
    # Between-item variation
    Omega.w  ~ dwish( R.w, 2 )
    ## R matrix:
    R.w[1,1] <- pow(sigma.c,2)
    R.w[2,2] <- pow(sigma.d,2)
    R.w[1,2] <- rho.w*sigma.c*sigma.d
    R.w[2,1] <- R.w[1,2]
    ## Vcov matrix:
    Sigma.w <- inverse(Omega.w)
    ## priors for var int. var slopes
    sigma.c ~ dunif(0,10)
    sigma.d ~ dunif(0,10)  
    ## prior for correlation:
    rho.w ~ dunif(-1,1)
}",
     file="gwmaximal.jag" )


###################################################
### code chunk number 4: 03_ADALecture3.Rnw:505-509
###################################################
track.variables<-c("beta","sigma.e",
                   "sigma.a","sigma.b",
                   "sigma.c","sigma.d",
                   "rho.u","rho.w")


###################################################
### code chunk number 5: 03_ADALecture3.Rnw:512-518
###################################################
library(rjags)
headnoun.mod <- jags.model( 
  file="gwmaximal.jag",
  data = headnoun.dat,
  n.chains = 4,
  n.adapt =2000 , quiet=T)


###################################################
### code chunk number 6: 03_ADALecture3.Rnw:526-530
###################################################
headnoun.res <- coda.samples(headnoun.mod,
                var = track.variables,
                n.iter = 10000,
                thin = 1)


###################################################
### code chunk number 7: 03_ADALecture3.Rnw:542-548
###################################################
mcmcChain<-as.matrix(headnoun.res)
#head(round(mcmcChain,digits=3))
hist(mcmcChain[,2],freq=F,main="Posterior distribution",xlab=expression(beta[1]))

x<-seq(-0.6,0.4,by=0.001)
lines(x,dnorm(x,mean=-0.0773,sd=0.1027))


###################################################
### code chunk number 8: 03_ADALecture3.Rnw:556-559
###################################################
## posterior probability of beta_1 < 0
## given data:
(meanbeta1<-mean(mcmcChain[,2]<0))


###################################################
### code chunk number 9: 03_ADALecture3.Rnw:589-590
###################################################
gelman.diag(headnoun.res)


###################################################
### code chunk number 10: 03_ADALecture3.Rnw:634-638
###################################################
op<-par(mfrow=c(1,2),pty="s")
hist(mcmcChain[,3],main=expression(rho[u]),
     xlab="")
hist(mcmcChain[,4],main=expression(rho[w]),xlab="")


###################################################
### code chunk number 11: 03_ADALecture3.Rnw:646-654
###################################################
op<-par(mfrow=c(1,2),pty="s")

x<-seq(0,100,by=0.01)
plot(x,dgamma(x,1.5,10^(-4)),type="l",main="Prior for sd")

x<-seq(-1,1,by=0.01)
rho<-dnorm(x)
plot(x,rho,type="l",main="Prior for correlation")


###################################################
### code chunk number 12: 03_ADALecture3.Rnw:667-739
###################################################
cat("
data
{
    zero[1] <- 0
    zero[2] <- 0
}
    model
{
    # Intercept and slope for each person, including random effects
    for( j in 1:J )
{
    u[j,1:2] ~ dmnorm(zero,Omega.u)
}    
    # Random effects for each item
    for( k in 1:K )
{
    w[k,1:2] ~ dmnorm(zero,Omega.w)
}
    
    # Define model for each observational unit
    for( i in 1:I )
{
    mu[i] <- ( beta[1] + u[subj[i],1] ) +
    ( beta[2] + u[subj[i],2] +  w[item[i],2]) * ( x[i] ) + w[item[i],1]
    rrt[i] ~ dnorm( mu[i], tau.e )
}
    # Priors:   
    # Fixed intercept and slope (uninformative)
    beta[1] ~ dnorm(0.0,1.0E-5)
    beta[2] ~ dnorm(0.0,1.0E-5)
    # Residual variance
    tau.e <- pow(sigma.e,-2)
    sigma.e  ~ dunif(0,100)    
    ## By subjects:
    Omega.u  ~ dwish( R.u, 2 )
    ## R matrix:
    R.u[1,1] <- pow(sigma.a,2)
    R.u[2,2] <- pow(sigma.b,2)
    R.u[1,2] <- rho.u*sigma.a*sigma.b
    R.u[2,1] <- R.u[1,2]
    ## Vcov matrix:
    Sigma.u <- inverse(Omega.u)
    ## priors for var int. var slopes
    tau.a ~ dgamma(1.5,10^(-4))
    tau.b ~ dgamma(1.5,10^(-4))
    sigma.a <- 1/pow(tau.a,1/2)
    sigma.b <- 1/pow(tau.b,1/2)
    ## prior for correlation:
#    rho.u2 ~ dbeta(1.5,1.5)
#    rho.u <- (2 * rho.u2) - 1
     ## truncated normal:
     rho.u ~ dnorm(0,1)T(-1,1)
# Between-item variation
    Omega.w  ~ dwish( R.w, 2 )
    ## R matrix:
    R.w[1,1] <- pow(sigma.c,2)
    R.w[2,2] <- pow(sigma.d,2)
    R.w[1,2] <- rho.w*sigma.c*sigma.d
    R.w[2,1] <- R.w[1,2]
    ## Vcov matrix:
    Sigma.w <- inverse(Omega.w)
    ## priors for var int. var slopes
    tau.c ~ dgamma(1.5,10^(-4))
    tau.d ~ dgamma(1.5,10^(-4))
    sigma.c <- 1/pow(tau.c,1/2)
    sigma.d <- 1/pow(tau.d,1/2)
    ## prior for correlation:
#    rho.w2 ~ dbeta(1.5,1.5)
#    rho.w <- (2 * rho.w2) - 1
     rho.w ~ dnorm(0,1)T(-1,1)
}",
     file="gwmaximal2.jag" )


###################################################
### code chunk number 13: 03_ADALecture3.Rnw:747-752
###################################################
headnoun.mod2 <- jags.model( 
  file="gwmaximal2.jag",
  data = headnoun.dat,
  n.chains = 4,
  n.adapt =2000 , quiet=T)


###################################################
### code chunk number 14: 03_ADALecture3.Rnw:760-764
###################################################
headnoun.res2 <- coda.samples(headnoun.mod2,
                var = track.variables,
                n.iter = 10000,
                thin = 20)


###################################################
### code chunk number 15: 03_ADALecture3.Rnw:774-775
###################################################
MCMCchain<-as.matrix(headnoun.res2)


###################################################
### code chunk number 16: 03_ADALecture3.Rnw:778-782
###################################################
par( mfrow=c(1,3) )
hist(MCMCchain[,2],xlab=expression(beta[1]),main="",freq=FALSE)
hist(MCMCchain[,3],xlab=expression(rho[u]),main="",freq=FALSE)
hist(MCMCchain[,4],xlab=expression(rho[w]),main="",freq=FALSE)


###################################################
### code chunk number 17: 03_ADALecture3.Rnw:791-792
###################################################
hist(MCMCchain[,2],xlab=expression(beta[1]),freq=FALSE,main="")


###################################################
### code chunk number 18: 03_ADALecture3.Rnw:795-796
###################################################
MCMCchain<-as.matrix(headnoun.res2)


###################################################
### code chunk number 19: 03_ADALecture3.Rnw:799-800
###################################################
mean(MCMCchain[,2]<0)


###################################################
### code chunk number 20: 03_ADALecture3.Rnw:838-917
###################################################
cat("
data
{
    zero[1] <- 0
    zero[2] <- 0
}
    model
{
    # Intercept and slope for each person, including random effects
    for( j in 1:J )
{
    u[j,1:2] ~ dmnorm(zero,Omega.u)
}    
    # Random effects for each item
    for( k in 1:K )
{
    w[k,1:2] ~ dmnorm(zero,Omega.w)
}
    
    # Define model for each observational unit
    for( i in 1:I )
{
    mu[i] <- ( beta[1] + u[subj[i],1] ) +
    ( beta[2] + u[subj[i],2] +  w[item[i],2]) * ( x[i] ) + w[item[i],1]
    rrt[i] ~ dnorm( mu[i], tau.e )
}
    # Priors:   
    beta[1] ~ dnorm(0.0,1.0E-5)
    beta[2] ~ dnorm(0.10,pow(0.10,-2))
    # Residual variance
    tau.e <- pow(sigma.e,-2)
    sigma.e  ~ dunif(0,100)    
    ## By subjects:
    Omega.u  ~ dwish( R.u, 2 )
    ## R matrix:
    R.u[1,1] <- pow(sigma.a,2)
    R.u[2,2] <- pow(sigma.b,2)
    R.u[1,2] <- rho.u*sigma.a*sigma.b
    R.u[2,1] <- R.u[1,2]
    ## Vcov matrix:
    Sigma.u <- inverse(Omega.u)
    ## priors for var int. var slopes
    tau.a ~ dgamma(1.5,10^(-4))
    tau.b ~ dgamma(1.5,10^(-4))
    sigma.a <- 1/pow(tau.a,1/2)
    sigma.b <- 1/pow(tau.b,1/2)
    ## prior for correlation:
    rho.u ~ dnorm(0,1)T(-1,1)
# Between-item variation
    Omega.w  ~ dwish( R.w, 2 )
    ## R matrix:
    R.w[1,1] <- pow(sigma.c,2)
    R.w[2,2] <- pow(sigma.d,2)
    R.w[1,2] <- rho.w*sigma.c*sigma.d
    R.w[2,1] <- R.w[1,2]
    ## Vcov matrix:
    Sigma.w <- inverse(Omega.w)
    ## priors for var int. var slopes
    tau.c ~ dgamma(1.5,10^(-4))
    tau.d ~ dgamma(1.5,10^(-4))
    sigma.c <- 1/pow(tau.c,1/2)
    sigma.d <- 1/pow(tau.d,1/2)
    ## prior for correlation:
     rho.w ~ dnorm(0,1)T(-1,1)
}",
     file="gwmaximal3.jag" )

headnoun.mod3 <- jags.model( 
  file="gwmaximal3.jag",
  data = headnoun.dat,
  n.chains = 4,
  n.adapt =2000 , quiet=T)

headnoun.res3 <- coda.samples(headnoun.mod3,
                var = track.variables,
                n.iter = 10000,
                thin = 20)

MCMCchain<-as.matrix(headnoun.res3)


###################################################
### code chunk number 21: 03_ADALecture3.Rnw:922-923
###################################################
mean(MCMCchain[,2]<0)