Homework 2 answers

I created this R markdown as a self contained document with code, output and conclusions. You don’t have to use this format. But if you want to learn to do this, feel free to contact me and I will try to explain. also, if you are an advanced R programmer experienced with R markdown language and want to share your knowledge, feel free to contribute through discussion topics.

setwd("C:/Users/marti/OneDrive/Documents/job/bayesian inference/2020/Lab 3 convergence diagnostics/")
library(coda)

load("MCMC.RData")
class(ll)
## [1] "mcmc.list"
colnames(ll[[1]])
## [1] "trt[1]"       "trt[2]"       "trt[3]"       "trt[4]"       "sigmablk"    
## [6] "sigmaepsilon"

Question 1: describe parameters of the MCMC sampler

This is done by askig for a summry, which includes a header with the required information

summary(ll)
## 
## Iterations = 101:5000
## Thinning interval = 1 
## Number of chains = 4 
## Sample size per chain = 4900 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                Mean     SD Naive SE Time-series SE
## trt[1]       81.303 14.234 0.101670          1.695
## trt[2]       82.288 14.281 0.102009          1.683
## trt[3]       86.303 14.242 0.101729          1.733
## trt[4]       83.309 14.251 0.101796          1.662
## sigmablk      9.089 21.311 0.152222          2.324
## sigmaepsilon  4.926  1.163 0.008306          0.014
## 
## 2. Quantiles for each variable:
## 
##                 2.5%    25%    50%    75%  97.5%
## trt[1]       20.5409 81.569 83.880 86.024 91.445
## trt[2]       20.5183 82.534 84.829 87.115 92.215
## trt[3]       25.7330 86.556 88.893 91.087 96.262
## trt[4]       22.3031 83.581 85.828 88.061 93.365
## sigmablk      0.5456  2.863  4.454  7.010 72.253
## sigmaepsilon  3.2222  4.104  4.728  5.533  7.712

There are 4 chains, burn in was 100, thin=1 and 4900 iterations were retained.

Question 2. Was burn in adequate?

Because we have 4 chains here, we can use gelman and rubin’s plot

gelman.plot(ll)

Clearly, burn in was not adequate. Especially for some treatment effects, the Gelman-Rubin statistic was well above 1.2.

Question 3 using similar supporting evidence as for 2.

From looking at the graphic of the Gelman and Rubin criteria, a burn in of 1000 may be apropriate.

subs<-window(ll,1001,5000)
summary(subs)
## 
## Iterations = 1001:5000
## Thinning interval = 1 
## Number of chains = 4 
## Sample size per chain = 4000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                Mean    SD Naive SE Time-series SE
## trt[1]       84.049 3.779  0.02987        0.08571
## trt[2]       85.047 3.765  0.02977        0.08868
## trt[3]       89.040 3.789  0.02996        0.09239
## trt[4]       86.053 3.826  0.03025        0.09045
## sigmablk      5.365 4.771  0.03772        0.11532
## sigmaepsilon  4.938 1.169  0.00924        0.01556
## 
## 2. Quantiles for each variable:
## 
##                 2.5%    25%    50%    75%  97.5%
## trt[1]       76.5309 81.927 84.036 86.117 91.661
## trt[2]       77.4728 82.940 85.001 87.239 92.441
## trt[3]       81.2701 86.876 89.054 91.208 96.454
## trt[4]       78.4911 83.887 86.004 88.204 93.601
## sigmablk      0.5239  2.783  4.303  6.531 16.726
## sigmaepsilon  3.2203  4.111  4.742  5.548  7.735
gelman.plot(subs)

This is clearly a good burn in.

Question 4.

effectiveSize(subs)
##       trt[1]       trt[2]       trt[3]       trt[4]     sigmablk sigmaepsilon 
##     2134.791     2088.989     1915.630     2044.037     1882.040     5742.360

This sample size is acceptable to estimate centrality measures such as means, medians, modes, and dispersion measures such as standard deviations. They are also acceptable to estimate quantiles, as long as they are not very extreme. For instance, I would feel comfortable estimating the percentile 5% or 95% of the posterior distribution with an effective sample size of 2000.

Question 5.

autocorr.diag(subs)
##            trt[1]       trt[2]     trt[3]     trt[4]    sigmablk sigmaepsilon
## Lag 0  1.00000000  1.000000000 1.00000000 1.00000000  1.00000000  1.000000000
## Lag 1  0.64742459  0.641857340 0.65300024 0.64430816  0.55469505  0.268678102
## Lag 5  0.27883551  0.283163305 0.28497624 0.28065212  0.24455484  0.071434663
## Lag 10 0.12964402  0.135620370 0.13956834 0.13156361  0.14378531  0.027719486
## Lag 50 0.02370604 -0.002518619 0.01042398 0.01370414 -0.02201192 -0.007150398

The lag10 autocorrelation ranges between 0.129 and 0.144 for treatment effects and block SD. it is substantiallylower (0.028 for the error standard deviation).

Question 6.

effectiveSize(ll)
##       trt[1]       trt[2]       trt[3]       trt[4]     sigmablk sigmaepsilon 
##     76.46100     76.50132     73.61455     78.27031     84.55938   6946.10584
autocorr.diag(ll)
##           trt[1]    trt[2]    trt[3]    trt[4]  sigmablk sigmaepsilon
## Lag 0  1.0000000 1.0000000 1.0000000 1.0000000 1.0000000  1.000000000
## Lag 1  0.9715370 0.9715061 0.9717644 0.9709444 0.8692769  0.267532346
## Lag 5  0.9295465 0.9304900 0.9300842 0.9294768 0.7693539  0.073527598
## Lag 10 0.8973059 0.8981073 0.8972537 0.8961364 0.7383326  0.029597411
## Lag 50 0.6722893 0.6720466 0.6714761 0.6708612 0.5441920 -0.004973321

The effective sample size is substaintially lower before we perform the appropriate burn in. The reason for it is because there is a very high autocorrelation.

And there is a very high autocorrelation because without and appropriate burn-in the trend in a changing mean for the MCMC chain leads to an inflated autocorrelation estimate. Remember that for correctly estimating autocorrelation in a time series it is better to have a de-trended series.

Question 7.

mc_lme<-read.table("MCMC_Example_1.txt",header=T)
lme_mc_list<-as.mcmc(mc_lme)
geweke.diag(lme_mc_list)
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##   varE   varU   h2_2 
## -1.106  1.199  1.221

This is only one chain, so Gelman and Rubin diagnostic can’t be computed. But we can compute Geweke’s. The Geweke statistic well below 2.0 in absolute value indicates that the burn in was likely appropriate.

Question 8.

summary(lme_mc_list)
## 
## Iterations = 1:12000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 12000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##        Mean      SD  Naive SE Time-series SE
## varE 0.5346 0.04327 0.0003950      0.0009275
## varU 0.5464 0.09286 0.0008477      0.0033972
## h2_2 0.5027 0.05344 0.0004878      0.0019567
## 
## 2. Quantiles for each variable:
## 
##        2.5%    25%    50%    75%  97.5%
## varE 0.4555 0.5050 0.5324 0.5626 0.6250
## varU 0.3851 0.4805 0.5400 0.6031 0.7521
## h2_2 0.3981 0.4662 0.5028 0.5388 0.6065
h2_hat=0.5464/(0.5346+0.5464)
h2_hat
## [1] 0.5054579

The point estimate of heritability is very close to the posterior meand and posterior median obtained from the full MCMC chain analysis (approximately 0.503). But with the full posterior distribution we are not limited to reporting these point estimates. We can also provide measures of variation such as posterior standard deviation and credibility intervals based on extreme quantiles.

Question 8.

library(rstan)
## Loading required package: StanHeaders
## Loading required package: ggplot2
## rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file
## 
## Attaching package: 'rstan'
## The following object is masked from 'package:coda':
## 
##     traceplot
y <- as.matrix(read.table('https://raw.github.com/wiki/stan-dev/rstan/rats.txt', header = TRUE))
dim(y)
## [1] 30  5
x <- c(8, 15, 22, 29, 36)
xbar <- mean(x)
N <- nrow(y)
T <- ncol(y)

rats<-list(y=y,x=x,xbar=xbar,N=nrow(y),T=ncol(y))

I start with 4 chains (to allow G-R diagnostics), 1000 for burn in, this is totally arbitrary and 5000 total iterations (again, arbitrary), thinning of 1 allows is to better assess autocorrelation and the need for thinning.

parameters_to_save=c("mu_alpha","mu_beta","sigmasq_y","sigmasq_alpha","sigmasq_beta","sigma_y",
                     "sigma_alpha","sigma_beta","alpha0")

rats_fit <- stan(file='rats.stan',
                 data=rats,chains = 4,warmup = 1000,iter = 5000,thin = 1,
                 pars=parameters_to_save)
## 
## SAMPLING FOR MODEL 'rats' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 1: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 1: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1001 / 5000 [ 20%]  (Sampling)
## Chain 1: Iteration: 1500 / 5000 [ 30%]  (Sampling)
## Chain 1: Iteration: 2000 / 5000 [ 40%]  (Sampling)
## Chain 1: Iteration: 2500 / 5000 [ 50%]  (Sampling)
## Chain 1: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 1: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 1: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 1: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 1: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 2.299 seconds (Warm-up)
## Chain 1:                2.988 seconds (Sampling)
## Chain 1:                5.287 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'rats' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 2: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 2: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1001 / 5000 [ 20%]  (Sampling)
## Chain 2: Iteration: 1500 / 5000 [ 30%]  (Sampling)
## Chain 2: Iteration: 2000 / 5000 [ 40%]  (Sampling)
## Chain 2: Iteration: 2500 / 5000 [ 50%]  (Sampling)
## Chain 2: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 2: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 2: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 2: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 2: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 1.502 seconds (Warm-up)
## Chain 2:                2.582 seconds (Sampling)
## Chain 2:                4.084 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'rats' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 3: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 3: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1001 / 5000 [ 20%]  (Sampling)
## Chain 3: Iteration: 1500 / 5000 [ 30%]  (Sampling)
## Chain 3: Iteration: 2000 / 5000 [ 40%]  (Sampling)
## Chain 3: Iteration: 2500 / 5000 [ 50%]  (Sampling)
## Chain 3: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 3: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 3: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 3: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 3: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 1.526 seconds (Warm-up)
## Chain 3:                2.442 seconds (Sampling)
## Chain 3:                3.968 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'rats' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 4: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 4: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1001 / 5000 [ 20%]  (Sampling)
## Chain 4: Iteration: 1500 / 5000 [ 30%]  (Sampling)
## Chain 4: Iteration: 2000 / 5000 [ 40%]  (Sampling)
## Chain 4: Iteration: 2500 / 5000 [ 50%]  (Sampling)
## Chain 4: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 4: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 4: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 4: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 4: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 2.283 seconds (Warm-up)
## Chain 4:                4.9 seconds (Sampling)
## Chain 4:                7.183 seconds (Total)
## Chain 4:
rats_mc<-As.mcmc.list(rats_fit)


gelman.plot(rats_mc)

This looks good in terms of burn-in. Effective sample size and autocorrelations help determinine if we need extra iterations and eventually thinning.

effectiveSize(rats_mc)
##      mu_alpha       mu_beta     sigmasq_y sigmasq_alpha  sigmasq_beta 
##     22085.714     19104.837      9230.954     18313.685     13909.326 
##       sigma_y   sigma_alpha    sigma_beta        alpha0          lp__ 
##      9227.027     19410.329     14110.448     20626.831      4965.835
autocorr.diag(rats_mc)
##            mu_alpha      mu_beta    sigmasq_y sigmasq_alpha sigmasq_beta
## Lag 0   1.000000000  1.000000000  1.000000000  1.0000000000 1.0000000000
## Lag 1  -0.168927369 -0.090476203  0.120038018 -0.1141653250 0.0006042779
## Lag 5  -0.005091094  0.005484355  0.015045245  0.0006088409 0.0059950915
## Lag 10 -0.006618910  0.005443036 -0.005909710 -0.0078721652 0.0070178692
## Lag 50  0.007332042  0.008024469 -0.006985847 -0.0025196550 0.0038073682
##             sigma_y   sigma_alpha  sigma_beta       alpha0         lp__
## Lag 0   1.000000000  1.0000000000 1.000000000  1.000000000  1.000000000
## Lag 1   0.120292393 -0.1389826126 0.007521068 -0.137356138  0.521692533
## Lag 5   0.015656734  0.0001612969 0.007149923  0.005838470  0.035500263
## Lag 10 -0.006011440 -0.0073117178 0.006565209  0.005420293 -0.005061834
## Lag 50 -0.007021498 -0.0048099370 0.003772125  0.004656116  0.001870252

These ESS are appropriate to estimate location and dispersion measures.