--- title: "Going Hierarchical" author: Jeff Rouder date: Sept. 28th, 2018 output: html_document --- A bunch of people are going to provide a bunch of responses: Let $Y_{ij}$ be the $jth$ response for the $i$th person. Here is some data: {r} set.seed(789) I=30 J=20 sub=rep(1:I,each=J) replicate=rep(1:J,I) true.mean=seq(.5,.55,length=I) #true.mean[sub] dat=rnorm(I*J,true.mean[sub],.15)  First, sample means. How well do sample means in the first half predict those in the second half {r} samp.mean=tapply(dat,sub,mean) plot(true.mean,samp.mean,xlim=c(.425,.625),ylim=c(.425,.625)) abline(0,1) err.sm=sqrt(sum((samp.mean-true.mean)^2)/I)  A nonhierarchical model: \begin{aligned} Y_{ij} &\sim \mbox{Normal}(\mu_i,\sigma^2)\\ \sigma^2 &\sim \mbox{Inverse Gamma}(r,s)\\ \mu_i &\sim \mbox{Normal}(a,b) \end{aligned} stan: {r} library(rstan) nonhier <- " data { int N; // #obs int I; // #ppl int sub[N]; // which person real y[N]; } parameters { vector[I] mu; real sigma2; } transformed parameters { real sigma; vector[N] obsmean; sigma = sqrt(sigma2); for (k in 1:N){ obsmean[k]=mu[sub[k]];} } model { mu ~ normal(.5,.3); sigma2 ~ inv_gamma(2,.2^2); y ~ normal(obsmean, sigma); }" N=length(dat) data <- list(N=N,I=I,sub=sub,y=dat) m <- stan_model(model_code = nonhier) samples <- sampling(m, data=data, iter=1000, chains=1, warmup=200) mu <- (extract(samples)$mu) sigma2 <- (extract(samples)$sigma2) nonhier.pm=apply(mu,2,mean)  {r} plot(true.mean,nonhier.pm,xlim=c(.425,.625),ylim=c(.425,.625)) abline(0,1) err.nonhier=sqrt(sum((nonhier.pm-true.mean)^2)/I) plot(samp.mean,nonhier.pm) abline(0,1)  Now, a hierarchical model: \begin{aligned} Y_{ij} &\sim \mbox{Normal}(\mu_i,\sigma^2)\\ \sigma^2 &\sim \mbox{Inverse Gamma}(r_1,s_1)\\ \mu_i &\sim \mbox{Normal}(\nu,\delta)\\ \nu &\sim \mbox{Normal}(a,b)\\ \delta & \sim \mbox{Inverse Gamma} \end{aligned} stan: {r} hier <- " data { int N; // #obs int I; // #ppl int sub[N]; // which person real y[N]; } parameters { vector[I] mu; real sigma2; real nu; real delta2; } transformed parameters { real sigma; real delta; vector[N] obsmean; sigma = sqrt(sigma2); delta=sqrt(delta2); for (k in 1:N){ obsmean[k]=mu[sub[k]];} } model { nu ~ normal(.5,.3); delta2 ~ inv_gamma(2,.1^2); sigma2 ~ inv_gamma(2,.2^2); mu ~normal(nu,delta); y ~ normal(obsmean, sigma); }" N=length(dat) data <- list(N=N,I=I,sub=sub,y=dat) mh <- stan_model(model_code = hier) samples <- sampling(mh, data=data, iter=1000, chains=1, warmup=200) mu <- (extract(samples)$mu) sigma2 <- (extract(samples)$sigma2) hier.pm=apply(mu,2,mean)  {r} plot(true.mean,hier.pm,xlim=c(.425,.625),ylim=c(.425,.625)) abline(0,1) err.hier=sqrt(sum((hier.pm-true.mean)^2)/I) plot(samp.mean,hier.pm) abline(0,1)  Same Thing in R! Nonhierarchical {r} library('MCMCpack') M=1000 mu=matrix(nrow=1000,ncol=I) s2=1:M #prior settings mu.m=.5 mu.v=.3^2 s2.scale=.2^2 s2.shape=2 ###start values mu[1,]=rep(.5,I) s2[1]=.1^2 for (m in 2:M){ #mu v=1/(J/s2[m-1]+1/mu.v) c=samp.mean*J/s2[m-1]+mu.m/mu.v mu[m,]=rnorm(I,v*c,sqrt(v)) #s2 scale=sum((dat-mu[m,sub])^2)/2+s2.scale s2[m]=rinvgamma(1,shape=s2.shape+N/2,scale=scale) } nonhier.pm=apply(mu,2,mean)  Hierarchical {r} M=1000 mu=matrix(nrow=1000,ncol=I) s2=1:M nu=1:M delta=1:M #prior settings nu.m=.5 nu.v=.3^2 s2.scale=.2^2 s2.shape=2 delta.scale=.1^2 delta.shape=2 ###start values mu[1,]=rep(.5,I) s2[1]=.1^2 nu[1]=.5 delta[1]=.1^2 for (m in 2:M){ #mu v=1/(J/s2[m-1]+1/delta[m-1]) c=samp.mean*J/s2[m-1]+nu[m-1]/delta[m-1] mu[m,]=rnorm(I,v*c,sqrt(v)) #s2 scale=sum((dat-mu[m,sub])^2)/2+s2.scale s2[m]=rinvgamma(1,shape=s2.shape+N/2,scale=scale) #nu v=1/(I/delta[m-1]+1/nu.v) c=I*mean(mu[m,])/delta[m-1]+nu.m/nu.v nu[m]=rnorm(1,v*c,sqrt(v)) #delta scale=sum((mu[m,]-nu[m])^2)/2+delta.scale delta[m]=rinvgamma(1,shape=delta.shape+I/2,scale=scale) } hier.pm=apply(mu,2,mean)  Your Turn * Group 1,2, manipulate I and J (low I, high J vs. low J, high I) * Group 3,4, manipuate true s2 vs true differences between ppl * Group 5,6, alternative ways of stan w/o warnings