--- title: "Which Hierarchical Model" author: Jeff Rouder date: Sept. 28th, 2018 output: html_document --- # First Problem There are often different model specifications with different shrinkage. Here is some data: {r} library('curl') urlName="http://pcl.missouri.edu/jeff/sites/pcl.missouri.edu.jeff/files/dat1_0.txt" dat=read.table(url(urlName),head=T) tapply(dat$rt,dat$cond,mean) sm=tapply(dat$rt,list(dat$sub,dat$cond),mean)  Your Job. Estimate each person-by-condition cell mean. * All Groups: Meet. Write a formal model(s). Put it on the board. Then implement it in R, stan, jags, chicken blood, whatever. {r} library('rstan') m1C <- " data { int n; int I; int J; vector[n] y; int sub[n]; int cond[n]; } parameters{ matrix[I,J] mu; real sigma2; } transformed parameters { real sigma; vector[n] obsMean; for (k in 1:n){obsMean[k]=mu[sub[k],cond[k]];} sigma = sqrt(sigma2);} model { for (kk in 1:I){ mu[kk,] ~ normal(.5,.5);} sigma2 ~ inv_gamma(.1,.1); y ~ normal(obsMean,sigma); }" n=length(dat$rt) I=length(unique(dat$sub)) J=length(unique(dat$cond)) data <- list(n=n, I=I, J=J, y=dat$rt, sub=dat$sub, cond=dat$cond) m1 <- stan_model(model_code = m1C) samples <- sampling(m1, data=data, iter=1000, chains=1) mu <- extract(samples)$mu sigma <- extract(samples)$sigma pm1=apply(mu,c(2,3),mean) person.pm1=apply(pm1,1,mean) residual1=pm1-person.pm1 matplot(t(residual1),typ='l') plot(sm,pm1) abline(0,1) t=as.matrix(read.table('truth1.txt')) err0=sqrt(mean((t-sm)^2)) err1= sqrt(mean((t-pm1)^2))  {r} m2C <- " data { int n; int I; int J; vector[n] y; int sub[n]; int cond[n]; } parameters{ matrix[I,J] mu; real sigma2; vector[I] nu; real tau2; } transformed parameters { real sigma; real tau; vector[n] obsMean; for (k in 1:n){obsMean[k]=mu[sub[k],cond[k]];} sigma = sqrt(sigma2); tau= sqrt(tau2);} model { nu~ normal(.5,.3); tau2~ inv_gamma(2,.2^2); for (kk in 1:I){ mu[kk,] ~ normal(nu[kk],tau);} sigma2 ~ inv_gamma(.1,.1); y ~ normal(obsMean,sigma); }" n=length(dat$rt) I=length(unique(dat$sub)) J=length(unique(dat$cond)) data <- list(n=n, I=I, J=J, y=dat$rt, sub=dat$sub, cond=dat$cond) m2 <- stan_model(model_code = m2C) samples <- sampling(m2, data=data, iter=1000, chains=1) mu <- extract(samples)$mu nu <- extract(samples)$nu sigma <- extract(samples)$sigma pm2=apply(mu,c(2,3),mean) matplot(t(pm2),typ='l') plot(sm,pm2) abline(0,1) person.pm2=apply(nu,2,mean) residual2=pm2-person.pm2 matplot(t(residual2),typ='l') err2= sqrt(mean((t-pm2)^2))  {r} m3C <- " data { int n; int I; int J; vector[n] y; int sub[n]; int cond[n]; } parameters{ matrix[I,J] mu; real sigma2; vector[I] nu; vector[J] beta; real muNu; real tau2; real rho2; } transformed parameters { real sigma; real tau; real rho; vector[n] obsMean; for (k in 1:n){obsMean[k]=mu[sub[k],cond[k]];} sigma = sqrt(sigma2); tau= sqrt(tau2); rho=sqrt(rho2);} model { muNu~ normal(.5,.3); tau2~inv_gamma(2,.2^2); rho2~ inv_gamma(2,.2^2); nu ~ normal(muNu,tau); beta ~normal(0,.1); for (kk in 1:I){ mu[kk,] ~ normal(nu[kk]+beta,rho);} sigma2 ~ inv_gamma(.1,.1); y ~ normal(obsMean,sigma); }" n=length(dat$rt) I=length(unique(dat$sub)) J=length(unique(dat$cond)) data <- list(n=n, I=I, J=J, y=dat$rt, sub=dat$sub, cond=dat$cond) m3 <- stan_model(model_code = m3C) samples <- sampling(m3, data=data, iter=1000, chains=1) mu <- extract(samples)$mu nu <- extract(samples)$nu sigma <- extract(samples)sigma pm3=apply(mu,c(2,3),mean) matplot(t(pm3),typ='l') plot(sm,pm3) abline(0,1) person.pm3=apply(nu,2,mean) residual3=pm3-person.pm3 matplot(t(residual3),typ='l') err3= sqrt(mean((t-pm3)^2))  # A Tree Let's Do One-Threshold Tree: One Condition, One Person: \begin{aligned} h &= d+(1-d)g\\ f &= g \end{aligned} More formally, letN_s$and$N_n$be the number of signal and noise trials, respectively, and let$Y_s$and$Y_nbe the respective number of "yes" responses in these trials. \begin{aligned} Y_s|d,g &\sim \mbox{Binomial}(d+(1-d)g,N_s)\\ Y_n|g &\sim \mbox{Binomial}(g,N_n)\\ d&\sim \mbox{Beta}(a,b)\\ g&\sim \mbox{Beta}(a,b) \end{aligned} {r} library(rstan) one <- " data { int nSig; // #number of signal trials int nNoise; // #number of noise trials int ySig; // #number of Hits int yNoise; // #number of FA } parameters { real d; real g; } transformed parameters{ real probH; real probF; probH=d+(1-d)*g; probF=g; } model { g ~ beta(.5, .5); d ~ beta(.5, .5); yNoise ~ binomial(nNoise,probF); ySig ~ binomial(nSig,probH); }" data <- list(nSig=100,nNoise=100,ySig=80,yNoise=20) m <- stan_model(model_code = one) samples <- sampling(m, data=data, iter=1000, chains=1, warmup=200) d <- (extract(samples)d) g <- (extract(samples)g)  # Many PPl Change in variables to probits \begin{aligned} Y_s|d^*,g^* &\sim \mbox{Binomial}(d^*+(1-d^*)g^*,N_s)\\ Y_n|g^* &\sim \mbox{Binomial}(g^*,N_n)\\ d^* &= \Phi(d)\\ g^* &= \Phi(g)\\ d&\sim \mbox{Normal}(a,b)\\ g&\sim \mbox{Normal}(a,b) \end{aligned} {r} one.probit <- " data { int nSig; // #number of signal trials int nNoise; // #number of noise trials int ySig; // #number of Hits int yNoise; // #number of FA } parameters { real d; real g; } transformed parameters{ real probH; real probF; real pd; real pg; pd=Phi(d); pg=Phi(g); probH=pd+(1-pd)*pg; probF=pg; } model { d ~ normal(0,1); g ~ normal(0,1); yNoise ~ binomial(nNoise,probF); ySig ~ binomial(nSig,probH); }" data <- list(nSig=100,nNoise=100,ySig=80,yNoise=20) m <- stan_model(model_code = one.probit) samples <- sampling(m, data=data, iter=1000, chains=1, warmup=200) d <- (extract(samples)d) g <- (extract(samples)g) lab=c(.01,.1,.3,.5,.7,.9,.99) plot(hist(d),axes=F,xlim=c(-3,3)) axis(1,at=qnorm(lab),lab=lab)  Now, time for ppl Change in variables to probits \begin{aligned} Y_{si}|d^*_i,g^*_i &\sim \mbox{Binomial}(d^*_i+(1-d^*_i)g^*_i,N_{si})\\ Y_{ni}|g^*_i &\sim \mbox{Binomial}(g^*_i,N_{ni})\\ d^*_i &= \Phi(d_i)\\ g^*_i &= \Phi(g_i)\\ d_i&\sim \mbox{Normal}(a,b)\\ g_i&\sim \mbox{Normal}(a,b) \end{aligned} {r} many <- " data { int I; //#number of ppl int nSig[I]; // #number of signal trials int nNoise[I]; // #number of noise trials int ySig[I]; // #number of Hits int yNoise[I]; // #number of FA } parameters { vector[I] d; vector[I] g; } transformed parameters{ vector[I] probH; vector[I] probF; vector[I] pd; vector[I] pg; for (i in 1:I){ pd[i]=Phi(d[i]); pg[i]=Phi(g[i]); probH[i]=pd[i]+(1-pd[i])*pg[i]; probF[i]=pg[i];} } model { d ~ normal(0,1); g ~ normal(0,1); yNoise ~ binomial(nNoise,probF); ySig ~ binomial(nSig,probH); }" I=20 t.s=rbeta(I,10,5) t.n=rbeta(I,5,15) nSig=nNoise=rep(100,I) ySig=rbinom(I,nSig,t.s) yNoise=rbinom(I,nNoise,t.n) data <- list(I=I,nSig=nSig,nNoise=nNoise,ySig=ySig,yNoise=yNoise) m <- stan_model(model_code = many) samples <- sampling(m, data=data, iter=1000, chains=1, warmup=200) d <- (extract(samples)d) g <- (extract(samples)$g) lab=c(.01,.1,.3,.5,.7,.9,.99) plot(hist(d),axes=F,xlim=c(-3,3)) axis(1,at=qnorm(lab),lab=lab)  Now make it hierarchical and compare. # Tree Models for Conditions: Suppose we have three conditions in the above high threshold model. We think these conditions should affect$d$and$g$. We also know the conditions were run "within," so, we think that$g$cannot depend on condition. Can we construct a high threshold model for several people and conditions? Here are the data: {r} urlName="http://pcl.missouri.edu/jeff/sites/pcl.missouri.edu.jeff/files/dat2.txt" dat=read.table(url(urlName),head=T) tapply(dat$resp,list(dat$cond,dat$stim),mean)  {r} library(rstan) code <- " data { int I; //# number of ppl int J; //# number of conditions int nSig[I,J]; // #number of signal trials int nNoise[I,J]; // #number of noise trials int ySig[I,J]; // #number of Hits int yNoise[I,J]; // #number of FA } parameters { vector[I] alpha; vector[J] beta; vector[I] g; real muAlpha; real muG; real s2Alpha; real s2G; } transformed parameters{ real sG; real sAlpha; matrix[I,J] d; matrix[I,J] probH; matrix[I,J] probF; matrix[I,J] pd; matrix[I,J] pg; for (i in 1:I){ for (j in 1:J){ d[i,j]=alpha[i]+beta[j]; pd[i,j]=Phi(d[i,j]); pg[i,j]=Phi(g[i]); probH[i,j]=pd[i,j]+(1-pd[i,j])*pg[i,j]; probF[i,j]=pg[i,j];}} sAlpha=sqrt(s2Alpha); sG=sqrt(s2G); } model { muAlpha ~ normal(0,1); muG ~ normal(0,1); s2Alpha ~ inv_gamma(2,.25); s2G ~ inv_gamma(2,.25); alpha ~ normal(muAlpha,sAlpha); beta ~ normal(0,.1); g ~ normal(muG,sG); for (i in 1:I){ for (j in 1:J){ yNoise[i,j] ~ binomial(nNoise[i,j],probF[i,j]); ySig[i,j] ~ binomial(nSig[i,j],probH[i,j]);}} }" m <- stan_model(model_code = code) I=max(dat$sub) J=max(dat$cond) nSig=as.matrix(table(list(dat$sub[dat$stim==1],dat$cond[dat$stim==1]))) nNoise=as.matrix(table(list(dat$sub[dat$stim==2],dat$cond[dat$stim==2]))) hits=tapply(dat$resp,list(dat$sub,dat$cond,dat$stim),sum) ySig=hits[,,1] yNoise=hits[,,2] data <- list(I=I,J=J,nSig=nSig,nNoise=nNoise,ySig=ySig,yNoise=yNoise) samples <- sampling(m, data=data, iter=1000, chains=1, warmup=200) d <- (extract(samples)$d) alpha <- colMeans(extract(samples)$alpha) beta <- extract(samples)$beta g <- colMeans(extract(samples)$g) plot.ecdf(pnorm(alpha)) plot.ecdf(pnorm(g)) lab=c(.01,.1,.3,.5,.7,.9,.99) plot(g,alpha,axes=F,xlim=c(-3,3),ylim=c(-3,3),pch=21,bg='green') axis(1,at=qnorm(lab),lab=lab) axis(2,at=qnorm(lab),lab=lab) boxplot(beta[,1],beta[,2],beta[,3],beta[,4])