lecture 2 code

gistfile1.txt

### R code from vignette source '02_ADALecture2.Rnw'

###################################################
### code chunk number 1: 02_ADALecture2.Rnw:66-72
###################################################
library(mvtnorm)
u0 <- u1 <- seq(from = -3, to = 3, length.out = 30)
Sigma1<-diag(2)
f <- function(u0, u1) dmvnorm(cbind(u0, u1), mean = c(0, 0),sigma = Sigma1)
z <- outer(u0, u1, FUN = f)
persp(u0, u1, z, theta = -30, phi = 30, ticktype = "detailed")


###################################################
### code chunk number 2: 02_ADALecture2.Rnw:81-85
###################################################
Sigma2<-matrix(c(1,.6,.6,1),byrow=FALSE,ncol=2)
f <- function(u0, u1) dmvnorm(cbind(u0, u1), mean = c(0, 0),sigma = Sigma2)
z <- outer(u0, u1, FUN = f)
persp(u0, u1, z, theta = -30, phi = 30, ticktype = "detailed")


###################################################
### code chunk number 3: 02_ADALecture2.Rnw:94-98
###################################################
Sigma3<-matrix(c(1,-.6,-.6,1),byrow=FALSE,ncol=2)
f <- function(u0, u1) dmvnorm(cbind(u0, u1), mean = c(0, 0),sigma = Sigma3)
z <- outer(u0, u1, FUN = f)
persp(u0, u1, z, theta = -30, phi = 30, ticktype = "detailed")


###################################################
### code chunk number 4: 02_ADALecture2.Rnw:230-231
###################################################
pnorm(600,mean=600,sd=sqrt(50))


###################################################
### code chunk number 5: 02_ADALecture2.Rnw:236-237
###################################################
qnorm(0.5,mean=600,sd=sqrt(50))


###################################################
### code chunk number 6: 02_ADALecture2.Rnw:242-243
###################################################
qnorm(0.5,mean=600,sd=sqrt(50),lower.tail=FALSE)


###################################################
### code chunk number 7: 02_ADALecture2.Rnw:255-256
###################################################
pnorm(610,mean=600,sd=sqrt(50))-pnorm(490,mean=600,sd=sqrt(50))


###################################################
### code chunk number 8: 02_ADALecture2.Rnw:265-271
###################################################
x<-rnorm(10000,mean=600,sd=sqrt(50))
## proportion of cases where 
## x is less than 500:
mean(x<590)
## theoretical value:
pnorm(590,mean=600,sd=sqrt(50))


###################################################
### code chunk number 9: 02_ADALecture2.Rnw:320-323
###################################################
plot(function(x) dnorm(x), -6, 6,
      main = "Normal density N(0,1)",ylim=c(0,.4),
              ylab="density",xlab="X")


###################################################
### code chunk number 10: 02_ADALecture2.Rnw:346-347
###################################################
pnorm(-2)


###################################################
### code chunk number 11: 02_ADALecture2.Rnw:362-366
###################################################
##P(Z < -x):
pnorm(-2)
##P(Z > x):
pnorm(2,lower.tail=FALSE)


###################################################
### code chunk number 12: 02_ADALecture2.Rnw:384-390
###################################################
## code source: http://www.r-bloggers.com/creating-shaded-areas-in-r/
coord.x <- c(-3,seq(-3,-2,0.01),-2) 
 coord.y <- c(0,dnorm(seq(-3,-2,0.01)),0) 
 curve(dnorm(x,0,1),xlim=c(-3,3),main='Standard Normal') 
 polygon(coord.x,coord.y,col='skyblue')
text(x=-2.5,y=0.1,label=expression(phi(-2)))


###################################################
### code chunk number 13: 02_ADALecture2.Rnw:410-411
###################################################
round(qnorm(0.025,lower.tail=FALSE),digits=2)


###################################################
### code chunk number 14: 02_ADALecture2.Rnw:445-446
###################################################
x<-rnorm(1,mean=600,sd=sqrt(50))


###################################################
### code chunk number 15: 02_ADALecture2.Rnw:451-453
###################################################
mu<-seq(400,800,by=0.1)
plot(mu,dnorm(x,mean=mu,sd=sqrt(50)),type="l")


###################################################
### code chunk number 16: 02_ADALecture2.Rnw:462-463
###################################################
x<-rnorm(10,mean=600,sd=sqrt(50))


###################################################
### code chunk number 17: 02_ADALecture2.Rnw:468-470
###################################################
## mu = 500
dnorm(x,mean=500,sd=sqrt(50))


###################################################
### code chunk number 18: 02_ADALecture2.Rnw:492-494
###################################################
## mu = 500
sum(dnorm(x,mean=500,sd=sqrt(50),log=TRUE))


###################################################
### code chunk number 19: 02_ADALecture2.Rnw:503-509
###################################################
mu<-seq(400,800,by=0.1)
liks<-rep(NA,length(mu))
for(i in 1:length(mu)){
liks[i]<-sum(dnorm(x,mean=mu[i],sd=sqrt(50),log=TRUE))
}
plot(mu,liks,type="l")


###################################################
### code chunk number 20: 02_ADALecture2.Rnw:568-572
###################################################
dbinom(x=46,size=100,0.4)
dbinom(x=46,size=100,0.46)
dbinom(x=46,size=100,0.5)
dbinom(x=46,size=100,0.6)


###################################################
### code chunk number 21: 02_ADALecture2.Rnw:580-583
###################################################
theta<-seq(0,1,by=0.01)
plot(theta,dbinom(x=46,size=100,theta),
     xlab=expression(theta),type="l")


###################################################
### code chunk number 22: betas
###################################################
plot(function(x) 
  dbeta(x,shape1=2,shape2=2), 0,1,
      main = "Beta density",
              ylab="density",xlab="X",ylim=c(0,4))

text(.5,1.3,"a=2,b=2")

plot(function(x) 
  dbeta(x,shape1=3,shape2=3),0,1,add=T)
text(.5,1.8,"a=3,b=3")


plot(function(x) 
  dbeta(x,shape1=6,shape2=6),0,1,add=T)
text(.5,2.6,"a=6,b=6")

plot(function(x) 
  dbeta(x,shape1=10,shape2=10),0,1,add=T)
text(.5,3.5,"a=60,b=60")


###################################################
### code chunk number 23: 02_ADALecture2.Rnw:772-792
###################################################
##lik:
plot(function(x) 
  dbeta(x,shape1=46,shape2=54),0,1,
              ylab="",xlab="X",col="red",
  ylim=c(0,10))

## prior:
plot(function(x) 
  dbeta(x,shape1=2,shape2=2), 0,1,
      main = "Prior",
              ylab="density",xlab="X",add=T,lty=2)

## posterior
plot(function(x) 
  dbeta(x,shape1=48,shape2=56), 0,1,
      main = "Posterior",
              ylab="density",xlab="X",add=T)

legend(0.1,6,legend=c("post","lik","prior"),
       lty=c(1,1,2),col=c("black","red","black"))


###################################################
### code chunk number 24: 02_ADALecture2.Rnw:801-821
###################################################
##lik:
plot(function(x) 
  dbeta(x,shape1=46,shape2=54),0,1,
              ylab="",xlab="X",col="red",
  ylim=c(0,10))

## prior:
plot(function(x) 
  dbeta(x,shape1=6,shape2=6), 0,1,
      main = "Prior",
              ylab="density",xlab="X",add=T,lty=2)

## posterior
plot(function(x) 
  dbeta(x,shape1=52,shape2=60), 0,1,
      main = "Posterior",
              ylab="density",xlab="X",add=T)

legend(0.1,6,legend=c("post","lik","prior"),
       lty=c(1,1,2),col=c("black","red","black"))


###################################################
### code chunk number 25: 02_ADALecture2.Rnw:830-850
###################################################
##lik:
plot(function(x) 
  dbeta(x,shape1=46,shape2=54),0,1,
              ylab="",xlab="X",col="red",
  ylim=c(0,10))

## prior:
plot(function(x) 
  dbeta(x,shape1=21,shape2=21), 0,1,
      main = "Prior",
              ylab="density",xlab="X",add=T,lty=2)

## posterior
plot(function(x) 
  dbeta(x,shape1=67,shape2=75), 0,1,
      main = "Posterior",
              ylab="density",xlab="X",add=T)

legend(0.1,6,legend=c("post","lik","prior"),
       lty=c(1,1,2),col=c("black","red","black"))


###################################################
### code chunk number 26: 02_ADALecture2.Rnw:892-896
###################################################
round(qbeta(0.025,shape1=100,shape2=100),digits=1)
round(qbeta(0.975,shape1=100,shape2=100),digits=1)
## ambivalent as to whether theta <0.5 or not:
round(pbeta(0.5,shape1=100,shape2=100),digits=1)


###################################################
### code chunk number 27: 02_ADALecture2.Rnw:922-923
###################################################
qbeta(0.5,shape1=589,shape2=611,lower.tail=FALSE)


###################################################
### code chunk number 28: 02_ADALecture2.Rnw:1035-1042
###################################################
beautydata<-read.table("data/beauty.txt",header=T)
## Note: beauty level is centered.
head(beautydata)

## restate the data as a list for JAGS:
data<-list(x=beautydata$beauty,
           y=beautydata$evaluation)


###################################################
### code chunk number 29: 02_ADALecture2.Rnw:1074-1091
###################################################
library(rjags)
cat("
model
   {
    ## specify model for data:
    for(i in 1:463){ 
    y[i] ~ dnorm(mu[i],tau)
    mu[i] <- beta0 + beta1 * (x[i])
    }
    # priors:
    beta0 ~ dunif(-10,10)
    beta1 ~ dunif(-10,10)
    sigma ~ dunif(0,100)
    sigma2 <- pow(sigma,2)
    tau <- 1/sigma2
   }",
     file="JAGSmodels/beautyexample1.jag" )


###################################################
### code chunk number 30: 02_ADALecture2.Rnw:1121-1138
###################################################
## specify which variables you want to examine
## the posterior distribution of:
track.variables<-c("beta0","beta1","sigma")

## define model:
beauty.mod <- jags.model( 
  file = "JAGSmodels/beautyexample1.jag",
                     data=data,
                     n.chains = 1,
                     n.adapt =2000, 
                     quiet=T)

## sample from posterior:
beauty.res <- coda.samples(beauty.mod,
                          var = track.variables,
                          n.iter = 2000,
                          thin = 1 ) 


###################################################
### code chunk number 31: 02_ADALecture2.Rnw:1145-1147
###################################################
round(summary(beauty.res)$statistics[,1:2],digits=2)
round(summary(beauty.res)$quantiles[,c(1,3,5)],digits=2)


###################################################
### code chunk number 32: 02_ADALecture2.Rnw:1155-1160
###################################################
lm_summary<-summary(lm(evaluation~beauty,
                       beautydata))

round(lm_summary$coef,digits=2)
round(lm_summary$sigma,digits=2)


###################################################
### code chunk number 33: 02_ADALecture2.Rnw:1171-1172
###################################################
densityplot(beauty.res)


###################################################
### code chunk number 34: 02_ADALecture2.Rnw:1181-1183
###################################################
op<-par(mfrow=c(1,3),pty="s")
traceplot(beauty.res)


###################################################
### code chunk number 35: 02_ADALecture2.Rnw:1192-1200
###################################################
MCMCsamp<-as.matrix(beauty.res)
op<-par(mfrow=c(1,3),pty="s")
hist(MCMCsamp[,1],main=expression(beta[0]),
     xlab="",freq=FALSE)
hist(MCMCsamp[,2],main=expression(beta[1]),
     xlab="",freq=FALSE)
hist(MCMCsamp[,3],main=expression(sigma),
     xlab="",freq=FALSE)


###################################################
### code chunk number 36: 02_ADALecture2.Rnw:1214-1216
###################################################
data<-list(x=c(8,15,22,29,36),
           y=c(177,236,285,350,376))


###################################################
### code chunk number 37: 02_ADALecture2.Rnw:1220-1222
###################################################
lm_summary_rats<-summary(fm<-lm(y~x,data))
round(lm_summary_rats$coef,digits=3)


###################################################
### code chunk number 38: 02_ADALecture2.Rnw:1231-1247
###################################################
cat("
model
   {
    ## specify model for data:
    for(i in 1:5){ 
    y[i] ~ dnorm(mu[i],tau)
    mu[i] <- beta0 + beta1 * (x[i]-mean(x[]))
    }
    # priors:
    beta0 ~ dunif(-500,500)
    beta1 ~ dunif(-500,500)
    tau <- 1/sigma2
    sigma2 <-pow(sigma,2)
    sigma ~ dunif(0,200)
   }",
     file="JAGSmodels/ratsexample2.jag" )


###################################################
### code chunk number 39: 02_ADALecture2.Rnw:1259-1268
###################################################
lexdec<-read.table("data/lexdec.txt",header=TRUE)
data<-lexdec[,c(1,2,3,4,5)]

contrasts(data$NativeLanguage)<-contr.sum(2)
contrasts(data$Sex)<-contr.sum(2)

lm_summary_lexdec<-summary(fm<-lm(RT~scale(Trial,scale=F)+
  NativeLanguage+Sex,data))
round(lm_summary_lexdec$coef,digits=2)


###################################################
### code chunk number 40: 02_ADALecture2.Rnw:1277-1281
###################################################
dat<-list(y=data$RT,
          Trial=(data$Trial-mean(data$Trial)),
          Lang=data$NativeLanguage,
          Sex=data$Sex)


###################################################
### code chunk number 41: 02_ADALecture2.Rnw:1289-1309
###################################################
cat("
model
   {
    ## specify model for data:
    for(i in 1:1659){ 
    y[i] ~ dnorm(mu[i],tau)
    mu[i] <- beta0 + 
             beta1 * Trial[i]+
             beta2 * Lang[i] +  beta3 * Sex[i] 
    }
    # priors:
    beta0 ~ dunif(-10,10)
    beta1 ~ dunif(-5,5)
    beta2 ~ dunif(-5,5)
    beta3 ~ dunif(-5,5)
    tau <- 1/sigma2
    sigma2 <-pow(sigma,2)
    sigma ~ dunif(0,200)
   }",
     file="JAGSmodels/multregexample1.jag" )


###################################################
### code chunk number 42: 02_ADALecture2.Rnw:1317-1331
###################################################
track.variables<-c("beta0","beta1",
                   "beta2","beta3","sigma")
library(rjags)

lexdec.mod <- jags.model( 
  file = "JAGSmodels/multregexample1.jag",
                     data=dat,
                     n.chains = 1,
                     n.adapt =2000, 
                      quiet=T)

lexdec.res <- coda.samples( lexdec.mod,
                                 var = track.variables,
                              n.iter = 3000)


###################################################
### code chunk number 43: 02_ADALecture2.Rnw:1338-1342
###################################################
round(summary(lexdec.res)$statistics[,1:2],
      digits=2)
round(summary(lexdec.res)$quantiles[,c(1,3,5)],
      digits=2)


###################################################
### code chunk number 44: 02_ADALecture2.Rnw:1351-1356
###################################################
lm_summary_lexdec<-summary(
  fm<-lm(RT~scale(Trial,scale=F)
         +NativeLanguage+Sex,
               lexdec))
round(lm_summary_lexdec$coef,digits=2)[,1:2]


###################################################
### code chunk number 45: 02_ADALecture2.Rnw:1376-1378
###################################################
beetledata<-read.table("data/beetle.txt",header=T)
head(beetledata)


###################################################
### code chunk number 46: 02_ADALecture2.Rnw:1387-1390
###################################################
dat<-list(x=beetledata$dose-mean(beetledata$dose),
          n=beetledata$number,
          y=beetledata$killed)


###################################################
### code chunk number 47: 02_ADALecture2.Rnw:1397-1409
###################################################
cat("
model
   {
for(i in 1:8){
    y[i] ~ dbin(p[i],n[i])
    logit(p[i]) <- beta0 + beta1 * x[i]
}
    # priors:
    beta0 ~ dunif(0,pow(100,2))
    beta1 ~ dunif(0,pow(100,2))
   }",
     file="JAGSmodels/glmexample1.jag" )


###################################################
### code chunk number 48: 02_ADALecture2.Rnw:1418-1434
###################################################
track.variables<-c("beta0","beta1")
## new:
inits <- list (list(beta0=0,
                    beta1=0))

glm.mod <- jags.model( 
  file = "JAGSmodels/glmexample1.jag",
                     data=dat,
                    ## new:
                    inits=inits,
                     n.chains = 1,
                      n.adapt =2000, quiet=T)

glm.res <- coda.samples( glm.mod,
                          var = track.variables,
                              n.iter = 2000) 


###################################################
### code chunk number 49: 02_ADALecture2.Rnw:1441-1445
###################################################
round(summary(glm.res)$statistics[,1:2],
      digits=2)
round(summary(glm.res)$quantiles[,c(1,3,5)],
      digits=2)


###################################################
### code chunk number 50: 02_ADALecture2.Rnw:1454-1458
###################################################
round(coef(glm(killed/number~scale(dose,scale=F),
         weights=number,
           family=binomial(),beetledata)),
      digits=2)


###################################################
### code chunk number 51: 02_ADALecture2.Rnw:1464-1465
###################################################
plot(glm.res)


###################################################
### code chunk number 52: 02_ADALecture2.Rnw:1494-1496
###################################################
data<-list(x=c(8,15,22,29,36),
           y=c(177,236,285,350,376))


###################################################
### code chunk number 53: 02_ADALecture2.Rnw:1503-1522
###################################################
cat("
model
   {
    ## specify model for data:
    for(i in 1:5){ 
    y[i] ~ dnorm(mu[i],tau)
    mu[i] <- beta0 + beta1 * (x[i]-mean(x[]))
    }
    ## prediction
    mu45 <- beta0+beta1 * (45-mean(x[]))
    y45 ~ dnorm(mu45,tau)
    # priors:
    beta0 ~ dunif(-500,500)
    beta1 ~ dunif(-500,500)
    tau <- 1/sigma2
    sigma2 <-pow(sigma,2)
    sigma ~ dunif(0,200)
   }",
     file="JAGSmodels/ratsexample2pred.jag" )


###################################################
### code chunk number 54: 02_ADALecture2.Rnw:1529-1541
###################################################
track.variables<-c("beta0","beta1","sigma","y45")

rats.mod <- jags.model( 
  file = "JAGSmodels/ratsexample2pred.jag",
                     data=data,
                     n.chains = 1,
                      n.adapt =2000, quiet=T)

rats.res <- coda.samples( rats.mod,
                          var = track.variables,
                              n.iter = 2000,
                                thin = 1) 


###################################################
### code chunk number 55: 02_ADALecture2.Rnw:1548-1550
###################################################
round(summary(rats.res)$statistics[,1:2],
       digits=2)


###################################################
### code chunk number 56: 02_ADALecture2.Rnw:1558-1559
###################################################
densityplot(rats.res)