# R Code 5: Hierarchical Basics

library(MCMCpack)
I = 100 # Number of People
J=10 # Number of Obs per person
N=I*J

#########################

set.seed(123)

#Set true values
#Try a lot of things here, you experiment
t.mu = seq(80,120,length=I)
#t.mu=rnorm(100,100,10)
#t.mu=80+rexp(100,.2)
#t.mu=c(rnorm(199,100,10),200)

t.s2=(20)^2

#################################
#make data
#Lets use a vector rather than matrix format.
sub=rep(1:I,each=J)
replicate=rep(1:J,I)

y=rnorm(N,t.mu[sub],sqrt(t.s2))

#quick check w/ sample means
sm=tapply(y,sub,mean)
plot(t.mu,sm)
RMSE.sm=sqrt(mean((sm-t.mu)^2))

##########################################
#Non Hierarchical Model
theta=0
delta=10000
a=.1
b=.1

M=1000

mu=matrix(nrow=M,ncol=I)
s2=1:M

#reasonable starting values
mu[1,]=rep(100,I)
s2[1]=100

for (m in 2:M)
{
#samle all mu's at once
var.prime=1/((J/s2[m-1])+(1/delta))
mean.prime=((J*sm/s2[m-1])+(theta/delta))*var.prime
mu[m,]=rnorm(I,mean.prime,sqrt(var.prime))
#sample s2
s2.scale=b+sum((y-mu[m,sub])^2)/2
s2[m]=rinvgamma(1,shape=a+N/2,scale=s2.scale)
}

keep=101:1000
plot(mu[keep,1],typ='l')
plot(s2[keep],typ='l')

nh=apply(mu[keep,],2,mean)
RMSE.nh=sqrt(mean((nh-t.mu)^2))

###########################
#hier
#these are now parameters, or more precisely, hyperparameters
theta=1:M
delta=1:M

#priors are needed for theta
theta.m=0
theta.v=10000

#priors needed for delta
delta.a=.1
delta.b=.1

#starting for theta,delta
#use old values for mu[1,] and s2[1]
theta[1]=50
delta[1]=100

for (m in 2:M)
{
#samle all mu's at once
var.prime=1/((J/s2[m-1])+(1/delta[m-1]))
mean.prime=((J*sm/s2[m-1])+(theta[m-1]/delta[m-1]))*var.prime
mu[m,]=rnorm(I,mean.prime,sqrt(var.prime))
#sample s2
s2.scale=b+sum((y-mu[m,sub])^2)/2
s2[m]=rinvgamma(1,shape=a+N/2,scale=s2.scale)
#hyperparameters
#sample theta
var.prime=1/((I/delta[m-1])+(1/theta.v))
mean.prime=((I*mean(mu[m,])/delta[m-1])+(theta.m/theta.v))*var.prime
theta[m]=rnorm(1,mean.prime,sqrt(var.prime))
#sample delta
delta.scale=delta.b+sum((mu[m,]-theta[m])^2)/2
delta[m]=rinvgamma(1,shape=delta.a+I/2,scale=delta.scale)
}

keep=101:1000
plot(mu[keep,1],typ='l')
plot(s2[keep],typ='l')

h=apply(mu[keep,],2,mean)
RMSE.h=sqrt(mean((h-t.mu)^2))

AttachmentSize
hier.pdf116.29 KB