Paris Lecture 4

ParisLecture4.R

### R code from vignette source '04_ADALecture4.Rnw'

###################################################
### code chunk number 1: 04_ADALecture4.Rnw:67-72
###################################################
vpost<-function(v=0.2609^2,n=1,s=0.15^2){
  return(1/((1/v)+n/s))
}
n<-seq(1,1000,by=1)
plot(n,vpost(v=2600,n=n),type="l",ylab="posterior variance",xlab="sample size")


###################################################
### code chunk number 2: 04_ADALecture4.Rnw:110-112
###################################################
x<-rnorm(1000,mean=0,sd=1)
mean(x)


###################################################
### code chunk number 3: 04_ADALecture4.Rnw:119-120
###################################################
mean(x<1.96)


###################################################
### code chunk number 4: 04_ADALecture4.Rnw:146-148
###################################################
indep.samp<-rnorm(500,mean=0,sd=1)
head(indep.samp)


###################################################
### code chunk number 5: indepsamp
###################################################
plot(1:500,indep.samp,type="l")


###################################################
### code chunk number 6: markovchainexample
###################################################
nsim<-500
x<-rep(NA,nsim)
y<-rep(NA,nsim)
x[1]<-rnorm(1) ## initialize x
for(i in 2:nsim){
## draw i-th value based on i-1-th value:  
y[i]<-rnorm(1,mean=x[i-1],sd=1)
x[i]<-y[i]
}


###################################################
### code chunk number 7: 04_ADALecture4.Rnw:183-184
###################################################
plot(1:nsim,y,type="l")


###################################################
### code chunk number 8: 04_ADALecture4.Rnw:235-238
###################################################
fn<-function(x){(1/40) * (2*x + 3)}
x<-seq(0,5,by=0.01)
plot(x,fn(x),type="l")


###################################################
### code chunk number 9: 04_ADALecture4.Rnw:301-304
###################################################
u<-runif(1000,min=0,max=1) 

z<-(1/2) * (sqrt(160*u +9)-3)


###################################################
### code chunk number 10: 04_ADALecture4.Rnw:313-317
###################################################
hist(z,freq=FALSE)
fn<-function(x){(1/40) * (2*x + 3)}
x<-seq(0,5,by=0.01)
lines(x,fn(x))


###################################################
### code chunk number 11: 04_ADALecture4.Rnw:360-366
###################################################
fn<-function(x){(1/40) * (2*x + 3)}
x<-seq(0,5,by=0.01)
plot(x,fn(x),type="l",ylim=c(0,.5))
arrows(x0=0,y0=0.4,x1=5,y1=0.4,code=0)
arrows(x0=0,y0=0,x1=0,y1=0.4,code=0)
arrows(x0=5,y0=0,x1=5,y1=0.4,code=0)


###################################################
### code chunk number 12: 04_ADALecture4.Rnw:376-377
###################################################
x<-runif(1,0,5)


###################################################
### code chunk number 13: 04_ADALecture4.Rnw:382-383
###################################################
fn(x)


###################################################
### code chunk number 14: 04_ADALecture4.Rnw:392-393
###################################################
u<-runif(1,0,1)


###################################################
### code chunk number 15: 04_ADALecture4.Rnw:398-400
###################################################
0.4*u
0.4*u < fn(x)


###################################################
### code chunk number 16: 04_ADALecture4.Rnw:411-426
###################################################
count<-0
k<-1 
accepted<-rep(NA,1000) 
rejected<-rep(NA,1000)
while(k<1001)
{
z<-runif(1,min=0,max=5) 
r<-((1/40)*(2*z+3))/(2*.2)
if(r>runif(1,min=0,max=1)) {
  accepted[k]<-z
  k<-k+1} else {
    rejected[k]<-z
  }
count<-count+1
}


###################################################
### code chunk number 17: rejectionsampling
###################################################
hist(accepted,freq=F,
     main="Example of rejection sampling")

fn<-function(x){
  (1/40)*(2*x+3)
}

x<-seq(0,5,by=0.01)

lines(x,fn(x))


###################################################
### code chunk number 18: 04_ADALecture4.Rnw:678-686
###################################################
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")
library(lme4)


###################################################
### code chunk number 19: 04_ADALecture4.Rnw:691-692
###################################################
m0<-lmer(rrt~x+(1|subj)+(1|item),headnoun)


###################################################
### code chunk number 20: 04_ADALecture4.Rnw:700-702
###################################################
library(car)
qqPlot(residuals(m0))


###################################################
### code chunk number 21: 04_ADALecture4.Rnw:711-712
###################################################
m1<-lmer(rt~x+(1|subj)+(1|item),headnoun)


###################################################
### code chunk number 22: 04_ADALecture4.Rnw:715-716
###################################################
qqPlot(residuals(m1))


###################################################
### code chunk number 23: 04_ADALecture4.Rnw:752-755
###################################################
head.means<-aggregate(rt~subj+type,
                      mean,data=headnoun)



###################################################
### code chunk number 24: 04_ADALecture4.Rnw:760-763
###################################################
t_test_result<-t.test(subset(head.means,type=="subj-ext")$rt,
       subset(head.means,type=="obj-ext")$rt,paired=T)
t_test_result$statistic;t_test_result$p.value


###################################################
### code chunk number 25: 04_ADALecture4.Rnw:797-808
###################################################
## number of simulations:
nsim<-100
## sample size:
n<-100
## predictor:
pred<-rep(c(0,1),each=n/2)
## matrices for storing results:
store<-matrix(0,nsim,4)
storenorm<-matrix(0,nsim,4)
## true effect:
beta.1<-0.5


###################################################
### code chunk number 26: 04_ADALecture4.Rnw:815-828
###################################################
for(i in 1:nsim){
errors<-rlnorm(n)
errorsnorm<-rnorm(n)
#errors<-errors-mean(errors)
#errorsnorm<-errorsnorm-mean(errorsnorm)
y<-5 + beta.1*pred + errors ## generate data:
ynorm<-5 + beta.1*pred + errorsnorm
fm<-lm(y~pred)
fmnorm<-lm(ynorm~pred)
## store coef., SE, t-value, p-value:
store[i,1:4] <- summary(fm)$coef[2,1:4]
storenorm[i,1:4] <- summary(fmnorm)$coef[2,1:4]
}


###################################################
### code chunk number 27: 04_ADALecture4.Rnw:835-838
###################################################
## ``observed'' power for raw scores:
mean(store[,4]<0.05)
mean(storenorm[,4]<0.05)


###################################################
### code chunk number 28: powerdifferential
###################################################
size<-seq(100,1000,by=10)
powerdiff<-rep(NA,length(size))
nsim<-100
for(j in 1:length(size)){
## sample size:
n<-size[j]
## predictor:
pred<-rep(c(0,1),each=n/2)
## matrices for storing results:
store<-matrix(0,nsim,4)
storenorm<-matrix(0,nsim,4)
## true effect:
beta.1<-0.5

for(i in 1:nsim){
errors<-rlnorm(n)
errorsnorm<-rnorm(n)
#errors<-errors-mean(errors)
#errorsnorm<-errorsnorm-mean(errorsnorm)
y<-5 + beta.1*pred + errors ## generate data:
ynorm<-5 + beta.1*pred + errorsnorm
fm<-lm(y~pred)
fmnorm<-lm(ynorm~pred)
## store coef., SE, t-value, p-value:
store[i,1:4] <- summary(fm)$coef[2,1:4]
storenorm[i,1:4] <- summary(fmnorm)$coef[2,1:4]
}
powerdiff[j]<-mean(storenorm[,4]<0.05)-mean(store[,4]<0.05)
}


###################################################
### code chunk number 29: 04_ADALecture4.Rnw:877-880
###################################################
plot(size,powerdiff,xlab="sample size",
     ylab="loss of power",
     main="why violating normality is bad for null results")


###################################################
### code chunk number 30: 04_ADALecture4.Rnw:889-893
###################################################
library(car)
op<-par(mfrow=c(1,2),pty="s")
qqPlot(residuals(fm))
qqPlot(residuals(fmnorm))


###################################################
### code chunk number 31: 04_ADALecture4.Rnw:903-905
###################################################
library(MASS)
boxcox(rt~type*subj,data=headnoun)