1/30/12

giving a guest lecture in baby bayes course for Brad on 2/2/12

# set up the data-here a1 represents individual level random effects

y1=matrix(0,10,4)
a1=rnorm(10,0,.5)
y1[,1]=a1+rnorm(10,4,1)
y1[,2]=a1+rnorm(10,4,1)
y1[,3]=a1+rnorm(10,5,1)
y1[,4]=a1+rnorm(10,5,1)

# data reduction 1-just use mean of replicates

y2=matrix(0,10,2)
y2[,1]=0.5*(y1[,1]+y1[,2])
y2[,2]=0.5*(y1[,3]+y1[,4])

# data reduction 2-look at differences over time for each subject

diff=y2[,2]-y2[,1]

# simple approach-use a 1 sample t-test

> t.test(diff)

        One Sample t-test

data:  dbar 
t = 2.8982, df = 9, p-value = 0.01765
alternative hypothesis: true mean is not equal to 0 
95 percent confidence interval:
 0.2494491 2.0237415 
sample estimates:
mean of x 
 1.136595 

# or use a Gibbs sampler

dbar=mean(diff)

# get an initial value for mu

rnorm(1,dbar,sqrt(var(diff)))
[1] 3.058596

N=10

postSmp1=matrix(0,1000,2)

postSmp1[1,2]=1/rgamma(1,0.5*N+1, 0.5*N*(3.1-dbar)^2+1)
postSmp1[1,1]=rnorm(1,dbar,sd=sqrt(postSmp1[1,2]/10))

for(i in 2:1000){
  postSmp1[i,2]=1/rgamma(1,0.5*N+1, 0.5*N*(postSmp1[i-1,1]-dbar)^2+1)
  postSmp1[i,1]=rnorm(1,dbar,sd=sqrt(postSmp1[i,2]/10))
}

> hist(postSmp1[,1]) 
> hist(postSmp1[,2])

> quantile(postSmp1[,1],c(.025,.975))
     2.5%     97.5% 
0.8292262 1.4194737 

# now suppose some data is missing-simulate that as iid

msInd=matrix(rbinom(40,1,.5),nrow=10)

b1=mean(c(y1[msInd[,1]==0,1],y1[msInd[,2]==0,2]))
a1=mean(c(y1[msInd[,3]==0,3],y1[msInd[,4]==0,4]))

yDat=matrix(0,N,4)
for(i in 1:N){
  for(j in 1:2){
    if(msInd[i,j]==1) yDat[i,j]=b1+rnorm(1,0,.2)
    else yDat[i,j]=y1[i,j]
  }
  for(j in 3:4){
    if(msInd[i,j]==1) yDat[i,j]=a1+rnorm(1,0,.2)
    else yDat[i,j]=y1[i,j]
  }
}

alpha=rnorm(N,0,1)
tau2=var(alpha)
sig2=var(y1[,1:2])
postSmp2=matrix(0,1000,2+N+2)

for(i in 1:1000){
  mu1=rnorm(1,mean(0.5*(yDat[,1]+yDat[,2])-alpha),sqrt(sig2/(2*N)))
  postSmp2[i,1]=mu1
  mu2=rnorm(1,mean(0.5*(yDat[,3]+yDat[,4])-alpha),sqrt(sig2/(2*N)))
  postSmp2[i,2]=mu2
  for(j in 1:N){
    mm=0.25*((yDat[j,1]-mu1)+(yDat[j,2]-mu1)+(yDat[j,3]-mu2)+(yDat[j,4]-mu2))
    mm=(mm/sig2)/(4/sig2+1/tau2)
    alpha[j]=rnorm(1,mm,sd=sqrt(1/(4/sig2+1/tau2)))
    postSmp2[i,2+j]=alpha[j]
  }
  tau2=1/rgamma(1,0.5*N+1, 0.5*sum(alpha^2)+1)
  postSmp2[i,2+N+1]=tau2
  mm=sum(yDat[,1]-mu1)^2+sum(yDat[,2]-mu1)^2+sum(yDat[,3]-mu2)^2+
    sum(yDat[,4]-mu2)^2
  sig2=1/rgamma(1,0.5*4*N+1, 0.5*mm+1)
  postSmp2[i,2+N+2]=sig2
  for(j in 1:N){
    for(k in 1:2){
      if(msInd[j,k]==1)
        yDat[j,k]=rnorm(1,mu1+alpha[j],sd=sqrt(sig2))
    }
    for(k in 3:4){
      if(msInd[j,k]==1)
        yDat[j,k]=rnorm(1,mu2+alpha[j],sd=sqrt(sig2))
    }
  }
}

then do

hist(postSmp2[,1])
hist(postSmp2[,2])
hist(postSmp2[,2]-postSmp2[,1])

quantile(postSmp2[,2]-postSmp2[,1],c(.025,.975))
      2.5%      97.5% 
0.09261415 1.83318409