Tutorial 7.4b - Single factor ANOVA (Bayesian)

The observed response ($y_i$) are assumed to be drawn from a normal distribution with a given mean ($\mu$) and standard deviation ($\sigma$). The expected values ($\mu$) are themselves determined by the linear predictor ($\beta_0 + \beta X_i$). In this case, $\beta_0$ represents the mean of the first group and the set of $\beta$'s represent the differences between each other group and the first group.

MCMC sampling requires priors on all parameters. We will employ weakly informative priors. Specifying 'uninformative' priors is always a bit of a balancing act. If the priors are too vague (wide) the MCMC sampler can wander off into nonscence areas of likelihood rather than concentrate around areas of highest likelihood (desired when wanting the outcomes to be largely driven by the data). On the other hand, if the priors are too strong, they may have an influence on the parameters. In such a simple model, this balance is very forgiving - it is for more complex models that prior choice becomes more important.

For this simple model, we will go with zero-centered Gaussian (normal) priors with relatively large standard deviations (100) for both the intercept and the treatment effect and a wide half-cauchy (scale=5) for the standard deviation. $$ \begin{align} y_i &\sim{} N(\mu_i, \sigma)\\ \mu_i &= \beta_0 + \beta X_i\\[1em] \beta_0 &\sim{} N(0,100)\\ \beta &\sim{} N(0,100)\\ \sigma &\sim{} cauchy(0,5)\\ \end{align} $$

Specific formulation

For very simple models such as this example, we can write the models as: $$\begin{align} y_i&\sim{}N(\mu_i, \tau)\\ \mu_i &= \beta_0 + \beta X_i\\ \beta_0&\sim{}N(0,1.0{E-6}) \hspace{1cm}\mathsf{non-informative~prior~for~interept}\\ \beta_j&\sim{}N(0,1.0{E-6}) \hspace{1cm}\mathsf{non-informative~prior~for~partial~slopes}\\ \tau &= 1/\sigma^2\\ \sigma&\sim{}U(0,100)\\ \end{align} $$

Define the model

Note the following example as group means calculated as derived posteriors

modelString = "
  model {
  #Likelihood
  for (i in 1:n) {
  y[i]~dnorm(mean[i],tau.res)
  mean[i] <- alpha+beta[x[i]]
  }

  #Priors and derivatives
  alpha ~ dnorm(0,1.0E-6)
  beta[1] <- 0
  for (i in 2:ngroups) {
  beta[i] ~ dnorm(0, 1.0E-6) #prior
  }
  sigma.res ~ dunif(0, 100)
  tau.res <- 1 / (sigma.res * sigma.res)
  sigma.group <- sd(beta[])

  #Group mean posteriors (derivatives)
  for (i in 1:ngroups) {
  Group.means[i] <- beta[i]+alpha
  }
  }
  "

Arrange the data as a list (as required by BUGS). As input, JAGS will need to be supplied with:

• the response variable (y)
• a numeric representation of the predictor variable (x)
• the total number of observed items (n)
• the number of groups

data.list <- with(data, list(y = y, x = as.numeric(x), n = nrow(data),
    ngroups = length(levels(data$x))))

Define the MCMC chain parameters

Next we should define the behavioural parameters of the MCMC sampling chains. Include the following:

• the nodes (estimated parameters) to monitor (return samples for)
• the number of MCMC chains (3)
• the number of burnin steps (1000)
• the thinning factor (10)
• the number of MCMC iterations - determined by the number of samples to save, the rate of thinning and the number of chains

params <- c("alpha", "beta", "sigma", "Group.means")
nChains = 3
burnInSteps = 3000
thinSteps = 10
numSavedSteps = 15000  #across all chains
nIter = ceiling(burnInSteps + (numSavedSteps * thinSteps)/nChains)
nIter

[1] 53000

Fit the model

Now run the JAGS code via the R2jags interface. Note that the first time jags is run after the R2jags package is loaded, it is often necessary to run any kind of randomization function just to initiate the .Random.seed variable.

## load the R2jags package
library(R2jags)

data.r2jags <- jags(data = data.list, inits = NULL, parameters.to.save = params,
    model.file = textConnection(modelString), n.chains = nChains, n.iter = nIter,
    n.burnin = burnInSteps, n.thin = thinSteps)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 50
   Unobserved stochastic nodes: 6
   Total graph size: 137

Initializing model

print(data.r2jags)

Inference for Bugs model at "5", fit using jags,
 3 chains, each with 53000 iterations (first 3000 discarded), n.thin = 10
 n.sims = 15000 iterations saved
               mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
Group.means[1]  40.390   0.847  38.688  39.838  40.389  40.955  42.065 1.001 15000
Group.means[2]  45.742   0.840  44.110  45.174  45.739  46.295  47.417 1.001 13000
Group.means[3]  54.591   0.844  52.961  54.026  54.589  55.141  56.299 1.001 15000
Group.means[4]  40.370   0.844  38.749  39.795  40.366  40.929  42.049 1.001  7700
Group.means[5]  30.396   0.844  28.747  29.833  30.394  30.949  32.057 1.001 15000
alpha           40.390   0.847  38.688  39.838  40.389  40.955  42.065 1.001 15000
beta[1]          0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
beta[2]          5.353   1.196   2.982   4.568   5.355   6.144   7.713 1.001 14000
beta[3]         14.201   1.190  11.894  13.410  14.195  14.980  16.597 1.001 15000
beta[4]         -0.020   1.196  -2.342  -0.823  -0.029   0.772   2.364 1.001 15000
beta[5]         -9.994   1.194 -12.317 -10.789 -10.004  -9.198  -7.615 1.001 15000
deviance       237.628   3.789 232.446 234.858 236.901 239.639 247.035 1.001 15000

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)
pD = 7.2 and DIC = 244.8
DIC is an estimate of expected predictive error (lower deviance is better).

data.mcmc.list <- as.mcmc(data.r2jags)

Model matrix formulation

For very simple models such as this example, we can write the models as: $$\begin{align} y_i&\sim{}N(\mu_i, \tau)\\ \mu_i &= \beta_0 + \beta X_i\\ \beta_0&\sim{}N(0,1.0{E-6}) \hspace{1cm}\mathsf{non-informative~prior~for~interept}\\ \beta_j&\sim{}N(0,1.0{E-6}) \hspace{1cm}\mathsf{non-informative~prior~for~partial~slopes}\\ \tau &= 1/\sigma^2\\ \sigma&\sim{}U(0,100)\\ \end{align} $$

Define the model

modelString = "
  model {
  #Likelihood
  for (i in 1:n) {
  y[i]~dnorm(mean[i],tau)
  mean[i] <- inprod(beta[],X[i,])
  }
  #Priors
  for (i in 1:ngroups) {
  beta[i] ~ dnorm(0, 1.0E-6) 
  }
  sigma ~ dunif(0, 100)
  tau <- 1 / (sigma * sigma)
  }
  "

Define the data

Arrange the data as a list (as required by BUGS). As input, JAGS will need to be supplied with:

• the response variable (y)
• the predictor model matrix (X)
• the total number of observed items (n)
• the number of predictor terms (nX)

X <- model.matrix(~x, data)
data.list <- with(data, list(y = y, X = X, n = nrow(data), ngroups = ncol(X)))

Define the MCMC chain parameters

Next we should define the behavioural parameters of the MCMC sampling chains. Include the following:

• the nodes (estimated parameters) to monitor (return samples for)
• the number of MCMC chains (3)
• the number of burnin steps (1000)
• the thinning factor (10)
• the number of MCMC iterations - determined by the number of samples to save, the rate of thinning and the number of chains

params <- c("beta", "sigma")
nChains = 3
burnInSteps = 3000
thinSteps = 10
numSavedSteps = 15000  #across all chains
nIter = ceiling(burnInSteps + (numSavedSteps * thinSteps)/nChains)
nIter

[1] 53000

Fit the model

Now run the JAGS code via the R2jags interface. Note that the first time jags is run after the R2jags package is loaded, it is often necessary to run any kind of randomization function just to initiate the .Random.seed variable.

data.r2jags <- jags(data = data.list, inits = NULL, parameters.to.save = params,
    model.file = textConnection(modelString), n.chains = nChains, n.iter = nIter,
    n.burnin = burnInSteps, n.thin = thinSteps)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 50
   Unobserved stochastic nodes: 6
   Total graph size: 380

Initializing model

print(data.r2jags)

Inference for Bugs model at "5", fit using jags,
 3 chains, each with 53000 iterations (first 3000 discarded), n.thin = 10
 n.sims = 15000 iterations saved
         mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
beta[1]   40.387   0.841  38.745  39.825  40.386  40.944  42.037 1.001 15000
beta[2]    5.363   1.190   3.058   4.556   5.363   6.164   7.710 1.001 15000
beta[3]   14.204   1.183  11.866  13.419  14.204  14.998  16.518 1.001 15000
beta[4]   -0.042   1.192  -2.358  -0.837  -0.031   0.745   2.283 1.001  8400
beta[5]   -9.987   1.193 -12.348 -10.780  -9.977  -9.199  -7.645 1.001 15000
sigma      2.646   0.287   2.157   2.447   2.619   2.817   3.284 1.001 15000
deviance 237.605   3.796 232.404 234.815 236.902 239.573 247.009 1.001 15000

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)
pD = 7.2 and DIC = 244.8
DIC is an estimate of expected predictive error (lower deviance is better).

data.mcmc.list <- as.mcmc(data.r2jags)

Whilst Gibbs sampling provides an elegantly simple MCMC sampling routine, very complex hierarchical models can take enormous numbers of iterations (often prohibitory large) to converge on a stable posterior distribution.

To address this, Andrew Gelman (and other collaborators) have implemented a variation on Hamiltonian Monte Carlo (HMC: a sampler that selects subsequent samples in a way that reduces the correlation between samples, thereby speeding up convergence) called the No-U-Turn (NUTS) sampler. All of these developments are brought together into a tool called Stan ("Sampling Through Adaptive Neighborhoods").

By design (to appeal to the vast BUGS users), Stan models are defined in a manner reminiscent of BUGS. Stan first converts these models into C++ code which is then compiled to allow very rapid computation.

Consistent with the use of C++, the model must be accompanied by variable declarations for all inputs and parameters.

One important difference between Stan and JAGS is that whereas BUGS (and thus JAGS) use precision rather than variance, Stan uses variance.

Stan itself is a stand-alone command line application. However, conveniently, the authors of Stan have also developed an R interface to Stan called Rstan which can be used much like R2jags.

Model matrix formulation

• the normal distribution is defined by variance rather than precision
• rather than using a uniform prior for sigma, I am using a half-Cauchy

We now translate the likelihood model into STAN code.
$$\begin{align} y_i&\sim{}N(\mu_i, \sigma)\\ \mu_i &= \beta_0+\beta X_i\\ \beta_0&\sim{}N(0,100)\\ \beta&\sim{}N(0,100)\\ \sigma&\sim{}Cauchy(0,5)\\ \end{align} $$

Define the model

modelString = "
  data {
  int<lower=1> n;
  int<lower=1> nX;
  vector [n] y;
  matrix [n,nX] X;
  }
  parameters {
  vector[nX] beta;
  real<lower=0> sigma;
  }
  transformed parameters {
  vector[n] mu;

  mu = X*beta;
  }
  model {
  #Likelihood
  y~normal(mu,sigma);
  
  #Priors
  beta ~ normal(0,1000);
  sigma~cauchy(0,5);
  }
  generated quantities {
  vector[n] log_lik;
  
  for (i in 1:n) {
  log_lik[i] = normal_lpdf(y[i] | mu[i], sigma); 
  }
  }
  "

Define the data

Arrange the data as a list (as required by BUGS). As input, JAGS will need to be supplied with:

• the response variable (y)
• the predictor model matrix (X)
• the total number of observed items (n)
• the number of predictor terms (nX)

Xmat <- model.matrix(~x, data)
data.list <- with(data, list(y = y, X = Xmat, nX = ncol(Xmat), n = nrow(data)))

Fit the model

Now run the JAGS code via the R2jags interface. Note that the first time jags is run after the R2jags package is loaded, it is often necessary to run any kind of randomization function just to initiate the .Random.seed variable.

## load the rstan package
library(rstan)

data.rstan <- stan(data = data.list, model_code = modelString, chains = 3,
    iter = 2000, warmup = 500, thin = 3)

In file included from /usr/local/lib/R/site-library/BH/include/boost/config.hpp:39:0,
                 from /usr/local/lib/R/site-library/BH/include/boost/math/tools/config.hpp:13,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core/var.hpp:7,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core.hpp:12,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/mat.hpp:4,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math.hpp:4,
                 from /usr/local/lib/R/site-library/StanHeaders/include/src/stan/model/model_header.hpp:4,
                 from file1c233b59b6c8.cpp:8:
/usr/local/lib/R/site-library/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
 #  define BOOST_NO_CXX11_RVALUE_REFERENCES
 ^
<command-line>:0:0: note: this is the location of the previous definition

SAMPLING FOR MODEL '3b057d3d81cbed2078ce678376a94574' NOW (CHAIN 1).

Gradient evaluation took 1.5e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  501 / 2000 [ 25%]  (Sampling)
Iteration:  700 / 2000 [ 35%]  (Sampling)
Iteration:  900 / 2000 [ 45%]  (Sampling)
Iteration: 1100 / 2000 [ 55%]  (Sampling)
Iteration: 1300 / 2000 [ 65%]  (Sampling)
Iteration: 1500 / 2000 [ 75%]  (Sampling)
Iteration: 1700 / 2000 [ 85%]  (Sampling)
Iteration: 1900 / 2000 [ 95%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 0.037412 seconds (Warm-up)
               0.084862 seconds (Sampling)
               0.122274 seconds (Total)


SAMPLING FOR MODEL '3b057d3d81cbed2078ce678376a94574' NOW (CHAIN 2).

Gradient evaluation took 1.1e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.11 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  501 / 2000 [ 25%]  (Sampling)
Iteration:  700 / 2000 [ 35%]  (Sampling)
Iteration:  900 / 2000 [ 45%]  (Sampling)
Iteration: 1100 / 2000 [ 55%]  (Sampling)
Iteration: 1300 / 2000 [ 65%]  (Sampling)
Iteration: 1500 / 2000 [ 75%]  (Sampling)
Iteration: 1700 / 2000 [ 85%]  (Sampling)
Iteration: 1900 / 2000 [ 95%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 0.035185 seconds (Warm-up)
               0.077264 seconds (Sampling)
               0.112449 seconds (Total)


SAMPLING FOR MODEL '3b057d3d81cbed2078ce678376a94574' NOW (CHAIN 3).

Gradient evaluation took 1.9e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.19 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  501 / 2000 [ 25%]  (Sampling)
Iteration:  700 / 2000 [ 35%]  (Sampling)
Iteration:  900 / 2000 [ 45%]  (Sampling)
Iteration: 1100 / 2000 [ 55%]  (Sampling)
Iteration: 1300 / 2000 [ 65%]  (Sampling)
Iteration: 1500 / 2000 [ 75%]  (Sampling)
Iteration: 1700 / 2000 [ 85%]  (Sampling)
Iteration: 1900 / 2000 [ 95%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 0.036548 seconds (Warm-up)
               0.061915 seconds (Sampling)
               0.098463 seconds (Total)

print(data.rstan, par = c("beta", "sigma"))

Inference for Stan model: 3b057d3d81cbed2078ce678376a94574.
3 chains, each with iter=2000; warmup=500; thin=3; 
post-warmup draws per chain=500, total post-warmup draws=1500.

         mean se_mean   sd   2.5%    25%   50%   75% 97.5% n_eff Rhat
beta[1] 40.40    0.02 0.84  38.75  39.85 40.38 40.95 42.06  1142    1
beta[2]  5.36    0.03 1.20   3.04   4.53  5.38  6.17  7.80  1302    1
beta[3] 14.21    0.03 1.19  11.88  13.37 14.20 15.01 16.68  1168    1
beta[4] -0.04    0.03 1.19  -2.35  -0.83 -0.06  0.79  2.23  1202    1
beta[5] -9.96    0.03 1.17 -12.14 -10.78 -9.95 -9.13 -7.76  1234    1
sigma    2.65    0.01 0.29   2.14   2.44  2.63  2.83  3.23  1332    1

Samples were drawn using NUTS(diag_e) at Mon Aug 28 20:56:23 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

The STAN team has put together pre-compiled modules (functions) to make specifying and applying STAN models much simpler. Each function offers a consistent interface that is Also reminiscent of major frequentist linear modelling routines in R.

Whilst it is not necessary to specify priors when using rstanarm functions (as defaults will be generated), there is no guarantee that the routines for determining these defaults will persist over time. Furthermore, it is always better to define your own priors if for no other reason that it forces you to thing about what you re doing. Consistent with the pure STAN version, we will employ the following priors:

• weakly informative Gaussian prior for the intercept $\beta_0 \sim{} N(0, 100)$
• weakly informative Gaussian prior for the treatment effect $\beta_1 \sim{} N(0, 100)$
• half-cauchy prior for the variance $\sigma \sim{} Cauchy(0, 5)$

Note, I am using the refresh=0 option so as to suppress the larger regular output in the interest of keeping output to what is necessary for this tutorial. When running outside of a tutorial context, the regular verbose output is useful as it provides a way to gauge progress.

library(rstanarm)
library(broom)
library(coda)

data.rstanarm = stan_glm(y ~ x, data = data, iter = 2000, warmup = 500,
    chains = 3, thin = 2, refresh = 0, prior_intercept = normal(0, 100),
    prior = normal(0, 100), prior_aux = cauchy(0, 2))

Gradient evaluation took 4.6e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.46 seconds.
Adjust your expectations accordingly!



 Elapsed Time: 0.121397 seconds (Warm-up)
               0.251909 seconds (Sampling)
               0.373306 seconds (Total)


Gradient evaluation took 1.2e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
Adjust your expectations accordingly!



 Elapsed Time: 0.117929 seconds (Warm-up)
               0.159803 seconds (Sampling)
               0.277732 seconds (Total)


Gradient evaluation took 1.7e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.17 seconds.
Adjust your expectations accordingly!



 Elapsed Time: 0.207204 seconds (Warm-up)
               0.153911 seconds (Sampling)
               0.361115 seconds (Total)

print(data.rstanarm)

stan_glm
 family:  gaussian [identity]
 formula: y ~ x
------

Estimates:
            Median MAD_SD
(Intercept)  40.4    0.8 
xB            5.3    1.2 
xC           14.2    1.2 
xD            0.0    1.1 
xE          -10.0    1.1 
sigma         2.6    0.3 

Sample avg. posterior predictive 
distribution of y (X = xbar):
         Median MAD_SD
mean_PPD 42.3    0.5  

------
For info on the priors used see help('prior_summary.stanreg').

tidyMCMC(data.rstanarm, conf.int = TRUE, conf.method = "HPDinterval")

         term    estimate std.error   conf.low conf.high
1 (Intercept) 40.37013202 0.8369156  38.641206 41.904030
2          xB  5.35134132 1.1832942   3.008367  7.741613
3          xC 14.23032089 1.1924590  11.961883 16.741929
4          xD  0.01681653 1.1703009  -2.498183  2.092344
5          xE -9.99518905 1.1917167 -12.239427 -7.563808
6       sigma  2.65469957 0.2833189   2.149518  3.229648

The brms package serves a similar goal to the rstanarm package - to provide a simple user interface to STAN. However, unlike the rstanarm implementation, brms simply converts the formula, data, priors and family into STAN model code and data before executing stan with those elements.

Whilst it is not necessary to specify priors when using brms functions (as defaults will be generated), there is no guarantee that the routines for determining these defaults will persist over time. Furthermore, it is always better to define your own priors if for no other reason that it forces you to thing about what you are doing. Consistent with the pure STAN version, we will employ the following priors:

• weakly informative Gaussian prior for the intercept $\beta_0 \sim{} N(0, 100)$
• weakly informative Gaussian prior for the treatment effect $\beta_1 \sim{} N(0, 100)$
• half-cauchy prior for the variance $\sigma \sim{} Cauchy(0, 5)$

Note, I am using the refresh=0. option so as to suppress the larger regular output in the interest of keeping output to what is necessary for this tutorial. When running outside of a tutorial context, the regular verbose output is useful as it provides a way to gauge progress.

library(brms)
library(broom)
library(coda)

data.brms = brm(y ~ x, data = data, iter = 2000, warmup = 500, chains = 3,
    thin = 2, refresh = 0, prior = c(prior(normal(0, 100), class = "Intercept"),
        prior(normal(0, 100), class = "b"), prior(cauchy(0, 5), class = "sigma")))

Gradient evaluation took 2.1e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.21 seconds.
Adjust your expectations accordingly!



 Elapsed Time: 0.033766 seconds (Warm-up)
               0.050702 seconds (Sampling)
               0.084468 seconds (Total)


Gradient evaluation took 8e-06 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
Adjust your expectations accordingly!



 Elapsed Time: 0.03249 seconds (Warm-up)
               0.053189 seconds (Sampling)
               0.085679 seconds (Total)


Gradient evaluation took 6e-06 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
Adjust your expectations accordingly!



 Elapsed Time: 0.034406 seconds (Warm-up)
               0.061338 seconds (Sampling)
               0.095744 seconds (Total)

print(data.brms)

 Family: gaussian(identity) 
Formula: y ~ x 
   Data: data (Number of observations: 50) 
Samples: 3 chains, each with iter = 2000; warmup = 500; thin = 2; 
         total post-warmup samples = 2250
    ICs: LOO = NA; WAIC = NA; R2 = NA
 
Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept    40.37      0.85    38.73    41.99       1830    1
xB            5.38      1.20     2.97     7.77       1834    1
xC           14.22      1.19    11.90    16.69       1816    1
xD            0.03      1.18    -2.30     2.33       1822    1
xE           -9.96      1.20   -12.30    -7.65       1575    1

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     2.63      0.29     2.13     3.24       1967    1

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

tidyMCMC(data.brms, conf.int = TRUE, conf.method = "HPDinterval")

         term   estimate std.error   conf.low conf.high
1 b_Intercept 40.3681829 0.8453842  38.733357 41.995095
2        b_xB  5.3835983 1.1963019   3.133155  7.926054
3        b_xC 14.2246055 1.1934948  12.002788 16.737062
4        b_xD  0.0322263 1.1761590  -2.398031  2.192691
5        b_xE -9.9603899 1.1972640 -12.385544 -7.745285
6       sigma  2.6319625 0.2902885   2.123581  3.222859

In addition to the regular model diagnostic checks (such as residual plots), for Bayesian analyses, it is necessary to explore the characteristics of the MCMC chains and the sampler in general. Recall that the purpose of MCMC sampling is to replicate the posterior distribution of the model likelihood and priors by drawing a known number of samples from this posterior (thereby formulating a probability distribution). This is only reliable if the MCMC samples accurately reflect the posterior.

Unfortunately, since we only know the posterior in the most trivial of circumstances, it is necessary to rely on indirect measures of how accurately the MCMC samples are likely to reflect the likelihood. I will breifly outline the most important diagnostics, however, please refer to Tutorial 4.3, Secton 3.1: Markov Chain Monte Carlo sampling for a discussion of these diagnostics.

• Traceplots for each parameter illustrate the MCMC sample values after each successive iteration along the chain. Bad chain mixing (characterized by any sort of pattern) suggests that the MCMC sampling chains may not have completely traversed all features of the posterior distribution and that more iterations are required to ensure the distribution has been accurately represented.
  
• Autocorrelation plot for each paramter illustrate the degree of correlation between MCMC samples separated by different lags. For example, a lag of 0 represents the degree of correlation between each MCMC sample and itself (obviously this will be a correlation of 1). A lag of 1 represents the degree of correlation between each MCMC sample and the next sample along the Chain and so on. In order to be able to generate unbiased estimates of parameters, the MCMC samples should be independent (uncorrelated). In the figures below, this would be violated in the top autocorrelation plot and met in the bottom autocorrelation plot.
  
• Rhat statistic for each parameter provides a measure of sampling efficiency/effectiveness. Ideally, all values should be less than 1.05. If there are values of 1.05 or greater it suggests that the sampler was not very efficient or effective. Not only does this mean that the sampler was potentiall slower than it could have been, more importantly, it could indicate that the sampler spent time sampling in a region of the likelihood that is less informative. Such a situation can arise from either a misspecified model or overly vague priors that permit sampling in otherwise nonscence parameter space.

Prior to inspecting any summaries of the parameter estimates, it is prudent to inspect a range of chain convergence diagnostics.

• Trace plots
  View trace plots

  library(MCMCpack)
  plot(data.mcmcpack)
  

  

  
  Trace plots show no evidence that the chains have not reasonably traversed the entire multidimensional parameter space.
• Raftery diagnostic
  View Raftery diagnostic

  library(MCMCpack)
  raftery.diag(data.mcmcpack)
  

  Quantile (q) = 0.025
  Accuracy (r) = +/- 0.005
  Probability (s) = 0.95 
                                                     
               Burn-in  Total Lower bound  Dependence
               (M)      (N)   (Nmin)       factor (I)
   (Intercept) 2        3865  3746         1.030     
   xB          2        3771  3746         1.010     
   xC          2        3802  3746         1.010     
   xD          2        3929  3746         1.050     
   xE          2        3865  3746         1.030     
   sigma2      2        3741  3746         0.999     
  

  The Raftery diagnostics estimate that we would require about 3900 samples to reach the specified level of confidence in convergence. As we have 10,000 samples, we can be confidence that convergence has occurred.
• Autocorrelation diagnostic
  View autocorrelations

  library(MCMCpack)
  autocorr.diag(data.mcmcpack)
  

          (Intercept)            xB           xC           xD           xE       sigma2
  Lag 0   1.000000000  1.000000e+00  1.000000000  1.000000000  1.000000000  1.000000000
  Lag 1  -0.002026451 -2.584595e-03  0.002845588  0.001221385 -0.001192370  0.113629174
  Lag 5   0.005186105  6.716212e-05 -0.005876489  0.001195213 -0.005830629 -0.002827759
  Lag 10 -0.004695655 -3.112666e-03 -0.017000405 -0.016469962 -0.005235883 -0.008952653
  Lag 50  0.001227380  1.032229e-02 -0.003024500 -0.001542249 -0.004604603 -0.014062511
  

  A lag of 1 appears to be mainly sufficient to avoid autocorrelation.

Again, prior to examining the summaries, we should have explored the convergence diagnostics.

library(coda)
data.mcmc = as.mcmc(data.r2jags)

• Trace plots

  plot(data.mcmc)
  

  

  
  Trace plots show no evidence that the chains have not reasonably traversed the entire multidimensional parameter space.

  When there are a lot of parameters, this can result in a very large number of traceplots. To focus on just certain parameters (such as $\beta$s)

  preds <- c("beta[1]", "beta[2]", "beta[3]", "beta[4]", "beta[5]")
  plot(as.mcmc(data.r2jags)[, preds])
  

  

  
• Raftery diagnostic

  raftery.diag(data.mcmc)
  

  [[1]]
  
  Quantile (q) = 0.025
  Accuracy (r) = +/- 0.005
  Probability (s) = 0.95 
                                                  
            Burn-in  Total Lower bound  Dependence
            (M)      (N)   (Nmin)       factor (I)
   beta[1]  20       36800 3746          9.82     
   beta[2]  20       39300 3746         10.50     
   beta[3]  20       36800 3746          9.82     
   beta[4]  20       38030 3746         10.20     
   beta[5]  20       36200 3746          9.66     
   deviance 20       37410 3746          9.99     
   sigma    20       36800 3746          9.82     
  
  
  [[2]]
  
  Quantile (q) = 0.025
  Accuracy (r) = +/- 0.005
  Probability (s) = 0.95 
                                                  
            Burn-in  Total Lower bound  Dependence
            (M)      (N)   (Nmin)       factor (I)
   beta[1]  20       36810 3746          9.83     
   beta[2]  20       38030 3746         10.20     
   beta[3]  20       36800 3746          9.82     
   beta[4]  20       36800 3746          9.82     
   beta[5]  20       38660 3746         10.30     
   deviance 20       38030 3746         10.20     
   sigma    20       38660 3746         10.30     
  
  
  [[3]]
  
  Quantile (q) = 0.025
  Accuracy (r) = +/- 0.005
  Probability (s) = 0.95 
                                                  
            Burn-in  Total Lower bound  Dependence
            (M)      (N)   (Nmin)       factor (I)
   beta[1]  20       35610 3746          9.51     
   beta[2]  20       39950 3746         10.70     
   beta[3]  20       37410 3746          9.99     
   beta[4]  20       36200 3746          9.66     
   beta[5]  20       38660 3746         10.30     
   deviance 20       38030 3746         10.20     
   sigma    20       38030 3746         10.20     
  

  The Raftery diagnostics for each chain estimate that we would require no more than 5000 samples to reach the specified level of confidence in convergence. As we have 16,667 samples, we can be confidence that convergence has occurred.
• Autocorrelation diagnostic

  autocorr.diag(data.mcmc)
  

               beta[1]       beta[2]      beta[3]      beta[4]       beta[5]      deviance
  Lag 0    1.000000000  1.0000000000  1.000000000  1.000000000  1.0000000000  1.0000000000
  Lag 10   0.000511858  0.0141941964  0.003049321 -0.001101961 -0.0136748526  0.0073553829
  Lag 50  -0.008312100 -0.0097313263 -0.003298050  0.003033837  0.0006659111  0.0002969917
  Lag 100 -0.006018466 -0.0008556030 -0.009141989 -0.001865391  0.0025599386 -0.0034075443
  Lag 500  0.014270778  0.0005677835  0.012646093  0.012218425 -0.0020017030  0.0002759488
                sigma
  Lag 0   1.000000000
  Lag 10  0.002197296
  Lag 50  0.010001090
  Lag 100 0.002718464
  Lag 500 0.007906380
  

  A lag of 10 appears to be sufficient to avoid autocorrelation (poor mixing).

Again, prior to examining the summaries, we should have explored the convergence diagnostics. There are numerous ways of working with STAN model fits (for exploring diagnostics and summarization).

• extract the mcmc samples and convert them into a mcmc.list to leverage the various coda routines
• use the numerous routines that come with the rstan package
• use the routines that come with the bayesplot package
• explore the diagnostics interactively via shinystan

• via coda
  • Traceplots

  library(coda)
  s = as.array(data.rstan)
  wch = grep("beta", dimnames(s)$parameters)
  s = s[, , wch]
  mcmc <- do.call(mcmc.list, plyr:::alply(s[, , -(length(s[1, 1, ]))], 2, as.mcmc))
  plot(mcmc)
  

  
  Trace plots show no evidence that the chains have not reasonably traversed the entire multidimensional parameter space.
  • Autocorrelation

  library(coda)
  s = as.array(data.rstan)
  wch = grep("beta", dimnames(s)$parameters)
  s = s[, , wch]
  mcmc <- do.call(mcmc.list, plyr:::alply(s[, , -(length(s[1, 1, ]))], 2, as.mcmc))
  autocorr.diag(mcmc)
  

               beta[1]    beta[2]      beta[3]       beta[4]
  Lag 0   1.000000e+00 1.00000000  1.000000000  1.0000000000
  Lag 1   8.980238e-02 0.03370305  0.086635647  0.0373446611
  Lag 5  -2.807365e-02 0.02098545  0.003301172  0.0004472765
  Lag 10 -7.017194e-03 0.01037456 -0.003796991 -0.0008328611
  Lag 50 -6.760447e-05 0.01494083 -0.018636469  0.0418145402
  

  Trace plots show no evidence that the chains have not reasonably traversed the entire multidimensional parameter space.
• via rstan
  • Traceplots

    stan_trace(data.rstan)
    

    
    Trace plots show no evidence that the chains have not reasonably traversed the entire multidimensional parameter space.
  • Raftery diagnostic

    raftery.diag(data.rstan)
    

    Quantile (q) = 0.025
    Accuracy (r) = +/- 0.005
    Probability (s) = 0.95 
    
    You need a sample size of at least 3746 with these values of q, r and s
    

    The Raftery diagnostics for each chain estimate that we would require no more than 5000 samples to reach the specified level of confidence in convergence. As we have 16,667 samples, we can be confidence that convergence has occurred.
  • Autocorrelation diagnostic

    stan_ac(data.rstan)
    

    
    A lag of 2 appears broadly sufficient to avoid autocorrelation (poor mixing).
  • Rhat values. These values are a measure of sampling efficiency/effectiveness. Ideally, all values should be less than 1.05. If there are values of 1.05 or greater it suggests that the sampler was not very efficient or effective. Not only does this mean that the sampler was potentiall slower than it could have been, more importantly, it could indicate that the sampler spent time sampling in a region of the likelihood that is less informative. Such a situation can arise from either a misspecified model or overly vague priors that permit sampling in otherwise nonscence parameter space.

    stan_rhat(data.rstan)
    

    
    In this instance, all rhat values are well below 1.05 (a good thing).
  • Another measure of sampling efficiency is Effective Sample Size (ess). ess indicate the number samples (or proportion of samples that the sampling algorithm deamed effective. The sampler rejects samples on the basis of certain criterion and when it does so, the previous sample value is used. Hence while the MCMC sampling chain may contain 1000 samples, if there are only 10 effective samples (1%), the estimated properties are not likely to be reliable.

    stan_ess(data.rstan)
    

    
    In this instance, most of the parameters have reasonably high effective samples and thus there is likely to be a good range of values from which to estimate paramter properties.
• via bayesplot
  • Trace plots and density plots

    library(bayesplot)
    mcmc_trace(as.matrix(data.rstan), regex_pars = "beta|sigma")
    

    

    library(bayesplot)
    mcmc_combo(as.matrix(data.rstan), regex_pars = "beta|sigma")
    

    
    Trace plots show no evidence that the chains have not reasonably traversed the entire multidimensional parameter space.
  • Density plots

    library(bayesplot)
    mcmc_dens(as.matrix(data.rstan), regex_pars = "beta|sigma")
    

    
    Density plots sugggest mean or median would be appropriate to describe the fixed posteriors and median is appropriate for the sigma posterior.
• via shinystan

  library(shinystan)
  launch_shinystan(data.rstan)
  

• It is worth exploring the influence of our priors.

Again, prior to examining the summaries, we should have explored the convergence diagnostics. There are numerous ways of working with STANARM model fits (for exploring diagnostics and summarization).

• extract the mcmc samples and convert them into a mcmc.list to leverage the various coda routines
• use the numerous routines that come with the rstan package
• use the routines that come with the bayesplot package
• explore the diagnostics interactively via shinystan

• via coda
  • Traceplots

  library(coda)
  s = as.array(data.rstanarm)
  mcmc <- do.call(mcmc.list, plyr:::alply(s[, , -(length(s[1, 1, ]))], 2, as.mcmc))
  plot(mcmc)
  

  

  
  Trace plots show no evidence that the chains have not reasonably traversed the entire multidimensional parameter space.
  • Autocorrelation

  library(coda)
  s = as.array(data.rstanarm)
  mcmc <- do.call(mcmc.list, plyr:::alply(s[, , -(length(s[1, 1, ]))], 2, as.mcmc))
  autocorr.diag(mcmc)
  

          (Intercept)          xB          xC          xD            xE
  Lag 0   1.000000000  1.00000000  1.00000000 1.000000000  1.0000000000
  Lag 1   0.027385218  0.01674547  0.01421235 0.006356097  0.0192142041
  Lag 5   0.005678036  0.02326828  0.01627096 0.020524683 -0.0148176836
  Lag 10 -0.015002246  0.01468946 -0.03007735 0.007219375  0.0074786588
  Lag 50 -0.024405944 -0.03054356  0.01621183 0.000803842 -0.0001606121
  

  Trace plots show no evidence that the chains have not reasonably traversed the entire multidimensional parameter space.
• via rstan
  • Traceplots

    stan_trace(data.rstanarm)
    

    
    Trace plots show no evidence that the chains have not reasonably traversed the entire multidimensional parameter space.
  • Raftery diagnostic

    raftery.diag(data.rstanarm)
    

    Quantile (q) = 0.025
    Accuracy (r) = +/- 0.005
    Probability (s) = 0.95 
    
    You need a sample size of at least 3746 with these values of q, r and s
    

    The Raftery diagnostics for each chain estimate that we would require no more than 5000 samples to reach the specified level of confidence in convergence. As we have 16,667 samples, we can be confidence that convergence has occurred.
  • Autocorrelation diagnostic

    stan_ac(data.rstanarm)
    

    
    A lag of 2 appears broadly sufficient to avoid autocorrelation (poor mixing).
  • Rhat values. These values are a measure of sampling efficiency/effectiveness. Ideally, all values should be less than 1.05. If there are values of 1.05 or greater it suggests that the sampler was not very efficient or effective. Not only does this mean that the sampler was potentiall slower than it could have been, more importantly, it could indicate that the sampler spent time sampling in a region of the likelihood that is less informative. Such a situation can arise from either a misspecified model or overly vague priors that permit sampling in otherwise nonscence parameter space.

    stan_rhat(data.rstanarm)
    

    
    In this instance, all rhat values are well below 1.05 (a good thing).
  • Another measure of sampling efficiency is Effective Sample Size (ess). ess indicate the number samples (or proportion of samples that the sampling algorithm deamed effective. The sampler rejects samples on the basis of certain criterion and when it does so, the previous sample value is used. Hence while the MCMC sampling chain may contain 1000 samples, if there are only 10 effective samples (1%), the estimated properties are not likely to be reliable.

    stan_ess(data.rstanarm)
    

    
    In this instance, most of the parameters have reasonably high effective samples and thus there is likely to be a good range of values from which to estimate paramter properties.
• via bayesplot
  • Trace plots and density plots

    mcmc_trace(as.array(data.rstanarm), regex_pars = "Intercept|x|sigma")
    

    

    mcmc_combo(as.array(data.rstanarm))
    

    
    Trace plots show no evidence that the chains have not reasonably traversed the entire multidimensional parameter space.
  • Density plots

    mcmc_dens(as.array(data.rstanarm))
    

    
    Density plots sugggest mean or median would be appropriate to describe the fixed posteriors and median is appropriate for the sigma posterior.
• via rstanarm
  The rstanarm package provides additional posterior checks.
  • Posterior vs Prior - this compares the posterior estimate for each parameter against the associated prior. If the spread of the priors is small relative to the posterior, then it is likely that the priors are too influential. On the other hand, overly wide priors can lead to computational issues.

    library(rstanarm)
    posterior_vs_prior(data.rstanarm, color_by = "vs", group_by = TRUE,
        facet_args = list(scales = "free_y"))
    

    Gradient evaluation took 3.9e-05 seconds
    1000 transitions using 10 leapfrog steps per transition would take 0.39 seconds.
    Adjust your expectations accordingly!
    
    
    
     Elapsed Time: 0.468186 seconds (Warm-up)
                   0.072428 seconds (Sampling)
                   0.540614 seconds (Total)
    
    
    Gradient evaluation took 1.1e-05 seconds
    1000 transitions using 10 leapfrog steps per transition would take 0.11 seconds.
    Adjust your expectations accordingly!
    
    
    
     Elapsed Time: 0.267889 seconds (Warm-up)
                   0.070828 seconds (Sampling)
                   0.338717 seconds (Total)
    

    
• via shinystan

  						  library(shinystan) 
  						  launch_shinystan(data.rstanarm))      
  

Again, prior to examining the summaries, we should have explored the convergence diagnostics. Rather than dublicate this for both additive and multiplicative models, we will only explore the multiplicative model. There are numerous ways of working with STAN model fits (for exploring diagnostics and summarization).

• extract the mcmc samples and convert them into a mcmc.list to leverage the various coda routines
• use the numerous routines that come with the rstan package
• use the routines that come with the bayesplot package
• explore the diagnostics interactively via shinystan

• via coda
  • Traceplots

  library(coda)
  mcmc = as.mcmc(data.brms)
  plot(mcmc)
  

  

  
  Trace plots show no evidence that the chains have not reasonably traversed the entire multidimensional parameter space.
  • Autocorrelation

  library(coda)
  mcmc = as.mcmc(data.brms)
  autocorr.diag(mcmc)
  

  Error in ts(x, start = start(x), end = end(x), deltat = thin(x)): invalid time series parameters specified
  

  Trace plots show no evidence that the chains have not reasonably traversed the entire multidimensional parameter space.
• via rstan
  • Traceplots

    stan_trace(data.brms$fit)
    

    
    Trace plots show no evidence that the chains have not reasonably traversed the entire multidimensional parameter space.
  • Raftery diagnostic

    raftery.diag(data.brms)
    

    Quantile (q) = 0.025
    Accuracy (r) = +/- 0.005
    Probability (s) = 0.95 
    
    You need a sample size of at least 3746 with these values of q, r and s
    

    The Raftery diagnostics for each chain estimate that we would require no more than 5000 samples to reach the specified level of confidence in convergence. As we have 16,667 samples, we can be confidence that convergence has occurred.
  • Autocorrelation diagnostic

    stan_ac(data.brms$fit)
    

    
    A lag of 2 appears broadly sufficient to avoid autocorrelation (poor mixing).
  • Rhat values. These values are a measure of sampling efficiency/effectiveness. Ideally, all values should be less than 1.05. If there are values of 1.05 or greater it suggests that the sampler was not very efficient or effective. Not only does this mean that the sampler was potentiall slower than it could have been, more importantly, it could indicate that the sampler spent time sampling in a region of the likelihood that is less informative. Such a situation can arise from either a misspecified model or overly vague priors that permit sampling in otherwise nonscence parameter space.

    stan_rhat(data.brms$fit)
    

    
    In this instance, all rhat values are well below 1.05 (a good thing).
  • Another measure of sampling efficiency is Effective Sample Size (ess). ess indicate the number samples (or proportion of samples that the sampling algorithm deamed effective. The sampler rejects samples on the basis of certain criterion and when it does so, the previous sample value is used. Hence while the MCMC sampling chain may contain 1000 samples, if there are only 10 effective samples (1%), the estimated properties are not likely to be reliable.

    stan_ess(data.brms$fit)
    

    
    In this instance, most of the parameters have reasonably high effective samples and thus there is likely to be a good range of values from which to estimate paramter properties.

Model validation involves exploring the model diagnostics and fit to ensure that the model is broadly appropriate for the data. As such, exploration of the residuals should be routine.

For more complex models (those that contain multiple effects, it is also advisable to plot the residuals against each of the individual predictors. For sampling designs that involve sample collection over space or time, it is also a good idea to explore whether there are any temporal or spatial patterns in the residuals.

There are numerous situations (e.g. when applying specific variance-covariance structures to a model) where raw residuals do not reflect the interior workings of the model. Typically, this is because they do not take into account the variance-covariance matrix or assume a very simple variance-covariance matrix. Since the purpose of exploring residuals is to evaluate the model, for these cases, it is arguably better to draw conclusions based on standardized (or studentized) residuals.

Unfortunately the definitions of standardized and studentized residuals appears to vary and the two terms get used interchangeably. I will adopt the following definitions:

The mark of a good model is being able to predict well. In an ideal world, we would have sufficiently large sample size as to permit us to hold a fraction (such as 25%) back thereby allowing us to train the model on 75% of the data and then see how well the model can predict the withheld 25%. Unfortunately, such a luxury is still rare in ecology.

The next best option is to see how well the model can predict the observed data. Models tend to struggle most with the extremes of trends and have particular issues when the extremes approach logical boundaries (such as zero for count data and standard deviations). We can use the fitted model to generate random predicted observations and then explore some properties of these compared to the actual observed data.

Residuals are not computed directly within MCMCpack. However, we can calculate them manually form the posteriors.

mcmc = as.data.frame(data.mcmcpack)
# generate a model matrix
newdata = data
Xmat = model.matrix(~x, newdata)
## get median parameter estimates
coefs = apply(mcmc[, 1:5], 2, median)
fit = as.vector(coefs %*% t(Xmat))
resid = data$y - fit
ggplot() + geom_point(data = NULL, aes(y = resid, x = fit))

Residuals against predictors

mcmc = as.data.frame(data.mcmcpack)
# generate a model matrix
newdata = newdata
Xmat = model.matrix(~x, newdata)
## get median parameter estimates
coefs = apply(mcmc[, 1:5], 2, median)
fit = as.vector(coefs %*% t(Xmat))
resid = data$y - fit
newdata = newdata %>% cbind(fit, resid)
ggplot(newdata) + geom_point(aes(y = resid, x = x))

And now for studentized residuals

mcmc = as.data.frame(data.mcmcpack)
# generate a model matrix
newdata = data
Xmat = model.matrix(~x, newdata)
## get median parameter estimates
coefs = apply(mcmc[, 1:5], 2, median)
fit = as.vector(coefs %*% t(Xmat))
resid = data$y - fit
sresid = resid/sd(resid)
ggplot() + geom_point(data = NULL, aes(y = sresid, x = fit))

Conclusions: for this simple model, the studentized residuals yield the same pattern as the raw residuals (or the Pearson residuals for that matter).

Lets see how well data simulated from the model reflects the raw data

mcmc = as.matrix(data.mcmcpack)
# generate a model matrix
Xmat = model.matrix(~x, data)
## get median parameter estimates
coefs = mcmc[, 1:5]
fit = coefs %*% t(Xmat)
## draw samples from this model
yRep = sapply(1:nrow(mcmc), function(i) rnorm(nrow(data), fit[i,
    ], sqrt(mcmc[i, "sigma2"])))
newdata = data.frame(x = data$x, yRep) %>% gather(key = Sample,
    value = Value, -x)
ggplot(newdata) + geom_violin(aes(y = Value, x = x, fill = "Model"),
    alpha = 0.5) + geom_violin(data = data, aes(y = y, x = x,
    fill = "Obs"), alpha = 0.5) + geom_point(data = data, aes(y = y,
    x = x), position = position_jitter(width = 0.1, height = 0),
    color = "black")

Conclusions the predicted trends do encapsulate the actual data, suggesting that the model is a reasonable representation of the underlying processes. Note, these are prediction intervals rather than confidence intervals as we are seeking intervals within which we can predict individual observations rather than means.

We can also explore the posteriors of each parameter.

library(bayesplot)
mcmc_intervals(as.matrix(data.mcmcpack), regex_pars = "Intercept|x|sigma")

mcmc_areas(as.matrix(data.mcmcpack), regex_pars = "Intercept|x|sigma")

Residuals are not computed directly within JAGS. However, we can calculate them manually form the posteriors.

mcmc = data.r2jags$BUGSoutput$sims.matrix %>% as.data.frame %>%
    dplyr:::select(contains("beta"), sigma) %>% as.matrix
# generate a model matrix
newdata = data
Xmat = model.matrix(~x, newdata)
## get median parameter estimates
coefs = apply(mcmc[, 1:5], 2, median)
fit = as.vector(coefs %*% t(Xmat))
resid = data$y - fit
ggplot() + geom_point(data = NULL, aes(y = resid, x = fit))

Residuals against predictors

mcmc = data.r2jags$BUGSoutput$sims.matrix %>% as.data.frame %>%
    dplyr:::select(contains("beta"), sigma) %>% as.matrix
# generate a model matrix
newdata = newdata
Xmat = model.matrix(~x, newdata)
## get median parameter estimates
coefs = apply(mcmc[, 1:5], 2, median)
fit = as.vector(coefs %*% t(Xmat))
resid = data$y - fit
newdata = newdata %>% cbind(fit, resid)
ggplot(newdata) + geom_point(aes(y = resid, x = x))

And now for studentized residuals

mcmc = data.r2jags$BUGSoutput$sims.matrix %>% as.data.frame %>%
    dplyr:::select(contains("beta"), sigma) %>% as.matrix
# generate a model matrix
newdata = data
Xmat = model.matrix(~x, newdata)
## get median parameter estimates
coefs = apply(mcmc[, 1:5], 2, median)
fit = as.vector(coefs %*% t(Xmat))
resid = data$y - fit
sresid = resid/sd(resid)
ggplot() + geom_point(data = NULL, aes(y = sresid, x = fit))

Conclusions: for this simple model, the studentized residuals yield the same pattern as the raw residuals (or the Pearson residuals for that matter).

Lets see how well data simulated from the model reflects the raw data

mcmc = data.r2jags$BUGSoutput$sims.matrix %>% as.data.frame %>%
    dplyr:::select(contains("beta"), sigma) %>% as.matrix
# generate a model matrix
Xmat = model.matrix(~x, data)
## get median parameter estimates
coefs = mcmc[, 1:5]
fit = coefs %*% t(Xmat)
## draw samples from this model
yRep = sapply(1:nrow(mcmc), function(i) rnorm(nrow(data), fit[i,
    ], mcmc[i, "sigma"]))
newdata = data.frame(x = data$x, yRep) %>% gather(key = Sample,
    value = Value, -x)
ggplot(newdata) + geom_violin(aes(y = Value, x = x, fill = "Model"),
    alpha = 0.5) + geom_violin(data = data, aes(y = y, x = x,
    fill = "Obs"), alpha = 0.5) + geom_point(data = data, aes(y = y,
    x = x), position = position_jitter(width = 0.1, height = 0),
    color = "black")

Conclusions the predicted trends do encapsulate the actual data, suggesting that the model is a reasonable representation of the underlying processes. Note, these are prediction intervals rather than confidence intervals as we are seeking intervals within which we can predict individual observations rather than means.

We can also explore the posteriors of each parameter.

library(bayesplot)
mcmc_intervals(data.r2jags$BUGSoutput$sims.matrix, regex_pars = "beta|sigma")

mcmc_areas(data.r2jags$BUGSoutput$sims.matrix, regex_pars = "beta|sigma")

Residuals are not computed directly within RSTAN. However, we can calculate them manually form the posteriors.

mcmc = as.data.frame(data.rstan) %>% dplyr:::select(contains("beta"),
    sigma) %>% as.matrix
# generate a model matrix
newdata = data
Xmat = model.matrix(~x, newdata)
## get median parameter estimates
coefs = apply(mcmc[, 1:5], 2, median)
fit = as.vector(coefs %*% t(Xmat))
resid = data$y - fit
ggplot() + geom_point(data = NULL, aes(y = resid, x = fit))

Residuals against predictors

mcmc = as.data.frame(data.rstan) %>% dplyr:::select(contains("beta"),
    sigma) %>% as.matrix
# generate a model matrix
newdata = newdata
Xmat = model.matrix(~x, newdata)
## get median parameter estimates
coefs = apply(mcmc[, 1:5], 2, median)
fit = as.vector(coefs %*% t(Xmat))
resid = data$y - fit
newdata = newdata %>% cbind(fit, resid)
ggplot(newdata) + geom_point(aes(y = resid, x = x))

And now for studentized residuals

mcmc = as.data.frame(data.rstan) %>% dplyr:::select(contains("beta"),
    sigma) %>% as.matrix
# generate a model matrix
newdata = data
Xmat = model.matrix(~x, newdata)
## get median parameter estimates
coefs = apply(mcmc[, 1:5], 2, median)
fit = as.vector(coefs %*% t(Xmat))
resid = data$y - fit
sresid = resid/sd(resid)
ggplot() + geom_point(data = NULL, aes(y = sresid, x = fit))

Conclusions: for this simple model, the studentized residuals yield the same pattern as the raw residuals (or the Pearson residuals for that matter).

Lets see how well data simulated from the model reflects the raw data

mcmc = as.data.frame(data.rstan) %>% dplyr:::select(contains("beta"),
    sigma) %>% as.matrix
# generate a model matrix
Xmat = model.matrix(~x, data)
## get median parameter estimates
coefs = mcmc[, 1:5]
fit = coefs %*% t(Xmat)
## draw samples from this model
yRep = sapply(1:nrow(mcmc), function(i) rnorm(nrow(data), fit[i,
    ], mcmc[i, "sigma"]))
newdata = data.frame(x = data$x, yRep) %>% gather(key = Sample,
    value = Value, -x)
ggplot(newdata) + geom_violin(aes(y = Value, x = x, fill = "Model"),
    alpha = 0.5) + geom_violin(data = data, aes(y = y, x = x,
    fill = "Obs"), alpha = 0.5) + geom_point(data = data, aes(y = y,
    x = x), position = position_jitter(width = 0.1, height = 0),
    color = "black")

Conclusions the predicted trends do encapsulate the actual data, suggesting that the model is a reasonable representation of the underlying processes. Note, these are prediction intervals rather than confidence intervals as we are seeking intervals within which we can predict individual observations rather than means.

We can also explore the posteriors of each parameter.

library(bayesplot)
mcmc_intervals(as.matrix(data.rstan), regex_pars = "beta|sigma")

mcmc_areas(as.matrix(data.rstan), regex_pars = "beta|sigma")

Residuals are not computed directly within RSTANARM. However, we can calculate them manually form the posteriors.

resid = resid(data.rstanarm)
fit = fitted(data.rstanarm)
ggplot() + geom_point(data = NULL, aes(y = resid, x = fit))

Residuals against predictors

resid = resid(data.rstanarm)
dat = data %>% mutate(resid = resid)
ggplot(dat) + geom_point(aes(y = resid, x = x))

And now for studentized residuals

resid = resid(data.rstanarm)
sresid = resid/sd(resid)
fit = fitted(data.rstanarm)
ggplot() + geom_point(data = NULL, aes(y = sresid, x = fit))

Conclusions: for this simple model, the studentized residuals yield the same pattern as the raw residuals (or the Pearson residuals for that matter).

Lets see how well data simulated from the model reflects the raw data

y_pred = posterior_predict(data.rstanarm)
newdata = data %>% cbind(t(y_pred)) %>% gather(key = "Rep", value = "Value",
    -y:-x)
ggplot(newdata) + geom_violin(aes(y = Value, x = x, fill = "Model"),
    alpha = 0.5) + geom_violin(data = data, aes(y = y, x = x,
    fill = "Obs"), alpha = 0.5) + geom_point(data = data, aes(y = y,
    x = x), position = position_jitter(width = 0.1, height = 0),
    color = "black")

Conclusions the predicted trends do encapsulate the actual data, suggesting that the model is a reasonable representation of the underlying processes. Note, these are prediction intervals rather than confidence intervals as we are seeking intervals within which we can predict individual observations rather than means.

We can also explore the posteriors of each parameter.

library(bayesplot)
mcmc_intervals(as.matrix(data.rstanarm), regex_pars = "Intercept|x|sigma")

mcmc_areas(as.matrix(data.rstanarm), regex_pars = "Intercept|x|sigma")

Residuals are not computed directly within BRMS. However, we can calculate them manually form the posteriors.

resid = resid(data.brms)[, "Estimate"]
fit = fitted(data.brms)[, "Estimate"]
ggplot() + geom_point(data = NULL, aes(y = resid, x = fit))

Residuals against predictors

resid = resid(data.brms)[, "Estimate"]
dat = data %>% mutate(resid = resid)
ggplot(dat) + geom_point(aes(y = resid, x = x))

And now for studentized residuals

resid = resid(data.brms)[, "Estimate"]
sresid = resid/sd(resid)
fit = fitted(data.brms)[, "Estimate"]
ggplot() + geom_point(data = NULL, aes(y = sresid, x = fit))

Conclusions: for this simple model, the studentized residuals yield the same pattern as the raw residuals (or the Pearson residuals for that matter).

Lets see how well data simulated from the model reflects the raw data

y_pred = posterior_predict(data.brms)
newdata = data %>% cbind(t(y_pred)) %>% gather(key = "Rep", value = "Value",
    -y:-x)
ggplot(newdata) + geom_violin(aes(y = Value, x = x, fill = "Model"),
    alpha = 0.5) + geom_violin(data = data, aes(y = y, x = x,
    fill = "Obs"), alpha = 0.5) + geom_point(data = data, aes(y = y,
    x = x), position = position_jitter(width = 0.1, height = 0),
    color = "black")

Conclusions the predicted trends do encapsulate the actual data, suggesting that the model is a reasonable representation of the underlying processes. Note, these are prediction intervals rather than confidence intervals as we are seeking intervals within which we can predict individual observations rather than means.

We can also explore the posteriors of each parameter.

library(bayesplot)
mcmc_intervals(as.matrix(data.brms), regex_pars = "b_|sigma")

mcmc_areas(as.matrix(data.brms), regex_pars = "b_|sigma")

Although all parameters in a Bayesian analysis are considered random and are considered a distribution, rarely would it be useful to present tables of all the samples from each distribution. On the other hand, plots of the posterior distributions are do have some use. Nevertheless, most workers prefer to present simple statistical summaries of the posteriors. Popular choices include the median (or mean) and 95% credibility intervals.

library(coda)
mcmcpvalue <- function(samp) {
    ## elementary version that creates an empirical p-value for the
    ## hypothesis that the columns of samp have mean zero versus a general
    ## multivariate distribution with elliptical contours.

    ## differences from the mean standardized by the observed
    ## variance-covariance factor

    ## Note, I put in the bit for single terms
    if (length(dim(samp)) == 0) {
        std <- backsolve(chol(var(samp)), cbind(0, t(samp)) - mean(samp),
            transpose = TRUE)
        sqdist <- colSums(std * std)
        sum(sqdist[-1] > sqdist[1])/length(samp)
    } else {
        std <- backsolve(chol(var(samp)), cbind(0, t(samp)) - colMeans(samp),
            transpose = TRUE)
        sqdist <- colSums(std * std)
        sum(sqdist[-1] > sqdist[1])/nrow(samp)
    }

}

Matrix model (MCMCpack)

summary(data.mcmcpack)

Iterations = 1001:11000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 10000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                 Mean     SD Naive SE Time-series SE
(Intercept)  40.39328 0.8282 0.008282       0.008282
xB            5.36170 1.1669 0.011669       0.011669
xC           14.21155 1.1879 0.011879       0.011879
xD           -0.02146 1.1689 0.011689       0.011689
xE          -10.00163 1.1704 0.011704       0.011704
sigma2        6.92503 1.5370 0.015370       0.017229

2. Quantiles for each variable:

               2.5%      25%       50%     75%  97.5%
(Intercept)  38.764  39.8386  40.39788 40.9314 42.036
xB            3.047   4.5933   5.36743  6.1326  7.620
xC           11.859  13.4340  14.21123 15.0097 16.478
xD           -2.306  -0.7992  -0.02423  0.7716  2.251
xE          -12.315 -10.7635 -10.00570 -9.2575 -7.657
sigma2        4.534   5.8257   6.70891  7.7744 10.511

# OR
library(broom)
tidyMCMC(data.mcmcpack, conf.int = TRUE, conf.method = "HPDinterval")

         term     estimate std.error   conf.low conf.high
1 (Intercept)  40.39327626 0.8282101  38.810155 42.070563
2          xB   5.36169548 1.1669217   3.142349  7.693206
3          xC  14.21155435 1.1879228  11.871023 16.485575
4          xD  -0.02146077 1.1688914  -2.258333  2.294295
5          xE -10.00163460 1.1704474 -12.336540 -7.688554
6      sigma2   6.92502993 1.5370186   4.213234  9.917678

• the mean of the first group (A) is 40.3932763
• the mean of the second group (B) is 5.3616955 units greater than (A)
• the mean of the third group (C) is 14.2115544 units greater than (A)
• the mean of the forth group (D) is -0.0214608 units greater (i.e. less) than (A)
• the mean of the fifth group (E) is -10.0016346 units greater (i.e. less) than (A)

While workers attempt to become comfortable with a new statistical framework, it is only natural that they like to evaluate and comprehend new structures and output alongside more familiar concepts. One way to facilitate this is via Bayesian p-values that are somewhat analogous to the frequentist p-values for investigating the hypothesis that a parameter is equal to zero.

## since values are less than zero
mcmcpvalue(data.mcmcpack[, 2])  # effect of (B-A)

[1] 0

mcmcpvalue(data.mcmcpack[, 3])  # effect of (C-A)

[1] 0

mcmcpvalue(data.mcmcpack[, 4])  # effect of (D-A)

[1] 0.9836

mcmcpvalue(data.mcmcpack[, 5])  # effect of (E-A)

[1] 0

mcmcpvalue(data.mcmcpack[, 2:5])  # effect of (all groups)

[1] 0

There is evidence that the reponse differs between the groups. There is very little evidence that the response of group D differs from that of group A

Matrix model (JAGS)

print(data.r2jags)

Inference for Bugs model at "5", fit using jags,
 3 chains, each with 53000 iterations (first 3000 discarded), n.thin = 10
 n.sims = 15000 iterations saved
         mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
beta[1]   40.387   0.841  38.745  39.825  40.386  40.944  42.037 1.001 15000
beta[2]    5.363   1.190   3.058   4.556   5.363   6.164   7.710 1.001 15000
beta[3]   14.204   1.183  11.866  13.419  14.204  14.998  16.518 1.001 15000
beta[4]   -0.042   1.192  -2.358  -0.837  -0.031   0.745   2.283 1.001  8400
beta[5]   -9.987   1.193 -12.348 -10.780  -9.977  -9.199  -7.645 1.001 15000
sigma      2.646   0.287   2.157   2.447   2.619   2.817   3.284 1.001 15000
deviance 237.605   3.796 232.404 234.815 236.902 239.573 247.009 1.001 15000

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)
pD = 7.2 and DIC = 244.8
DIC is an estimate of expected predictive error (lower deviance is better).

# OR
library(broom)
tidyMCMC(as.mcmc(data.r2jags), conf.int = TRUE, conf.method = "HPDinterval")

      term     estimate std.error   conf.low  conf.high
1  beta[1]  40.38727065 0.8407378  38.690751  41.969038
2  beta[2]   5.36294219 1.1902714   3.079022   7.723797
3  beta[3]  14.20375455 1.1833589  11.799076  16.428188
4  beta[4]  -0.04249084 1.1921215  -2.355980   2.286026
5  beta[5]  -9.98747340 1.1933066 -12.371762  -7.677127
6 deviance 237.60478820 3.7960918 231.801132 245.265711
7    sigma   2.64600258 0.2866295   2.130840   3.230809

• the mean of the first group (A) is 40.3872707
• the mean of the second group (B) is 5.3629422 units greater than (A)
• the mean of the third group (C) is 14.2037546 units greater than (A)
• the mean of the forth group (D) is -0.0424908 units greater (i.e. less) than (A)
• the mean of the fifth group (E) is -9.9874734 units greater (i.e. less) than (A)

While workers attempt to become comfortable with a new statistical framework, it is only natural that they like to evaluate and comprehend new structures and output alongside more familiar concepts. One way to facilitate this is via Bayesian p-values that are somewhat analogous to the frequentist p-values for investigating the hypothesis that a parameter is equal to zero.

## since values are less than zero
mcmcpvalue(data.r2jags$BUGSoutput$sims.matrix[, "beta[2]"])  # effect of (B-A)

[1] 0

mcmcpvalue(data.r2jags$BUGSoutput$sims.matrix[, "beta[3]"])  # effect of (C-A)

[1] 0

mcmcpvalue(data.r2jags$BUGSoutput$sims.matrix[, "beta[4]"])  # effect of (D-A)

[1] 0.9718667

mcmcpvalue(data.r2jags$BUGSoutput$sims.matrix[, "beta[5]"])  # effect of (E-A)

[1] 0

mcmcpvalue(data.r2jags$BUGSoutput$sims.matrix[, 2:5])  # effect of (all groups)

[1] 0

There is evidence that the reponse differs between the groups. There is very little evidence that the response of group D differs from that of group A

Matrix model (RSTAN)

print(data.rstan, pars = c("beta", "sigma"))

Inference for Stan model: 3b057d3d81cbed2078ce678376a94574.
3 chains, each with iter=2000; warmup=500; thin=3; 
post-warmup draws per chain=500, total post-warmup draws=1500.

         mean se_mean   sd   2.5%    25%   50%   75% 97.5% n_eff Rhat
beta[1] 40.40    0.02 0.84  38.75  39.85 40.38 40.95 42.06  1142    1
beta[2]  5.36    0.03 1.20   3.04   4.53  5.38  6.17  7.80  1302    1
beta[3] 14.21    0.03 1.19  11.88  13.37 14.20 15.01 16.68  1168    1
beta[4] -0.04    0.03 1.19  -2.35  -0.83 -0.06  0.79  2.23  1202    1
beta[5] -9.96    0.03 1.17 -12.14 -10.78 -9.95 -9.13 -7.76  1234    1
sigma    2.65    0.01 0.29   2.14   2.44  2.63  2.83  3.23  1332    1

Samples were drawn using NUTS(diag_e) at Mon Aug 28 20:56:23 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

# OR
library(broom)
tidyMCMC(data.rstan, conf.int = TRUE, conf.method = "HPDinterval", pars = c("beta", "sigma"))

     term    estimate std.error   conf.low conf.high
1 beta[1] 40.39999711 0.8365319  38.789811 42.074467
2 beta[2]  5.36368555 1.2032887   3.040275  7.814264
3 beta[3] 14.21257316 1.1862070  11.751943 16.445171
4 beta[4] -0.04359039 1.1940964  -2.328043  2.232564
5 beta[5] -9.96494917 1.1732056 -12.145567 -7.749614
6   sigma  2.64635314 0.2850148   2.112832  3.180584

• the mean of the first group (A) is 40.3999971
• the mean of the second group (B) is 5.3636856 units greater than (A)
• the mean of the third group (C) is 14.2125732 units greater than (A)
• the mean of the forth group (D) is -0.0435904 units greater (i.e. less) than (A)
• the mean of the fifth group (E) is -9.9649492 units greater (i.e. less) than (A)

While workers attempt to become comfortable with a new statistical framework, it is only natural that they like to evaluate and comprehend new structures and output alongside more familiar concepts. One way to facilitate this is via Bayesian p-values that are somewhat analogous to the frequentist p-values for investigating the hypothesis that a parameter is equal to zero.

## since values are less than zero
mcmcpvalue(as.matrix(data.rstan)[, "beta[2]"])  # effect of (B-A)

[1] 0

mcmcpvalue(as.matrix(data.rstan)[, "beta[3]"])  # effect of (C-A)

[1] 0

mcmcpvalue(as.matrix(data.rstan)[, "beta[4]"])  # effect of (D-A)

[1] 0.9766667

mcmcpvalue(as.matrix(data.rstan)[, "beta[5]"])  # effect of (E-A)

[1] 0

mcmcpvalue(as.matrix(data.rstan)[, 2:5])  # effect of (all groups)

[1] 0

There is evidence that the reponse differs between the groups. There is very little evidence that the response of group D differs from that of group A

library(loo)
(full = loo(extract_log_lik(data.rstan)))

Computed from 1500 by 50 log-likelihood matrix

         Estimate   SE
elpd_loo   -122.3  5.9
p_loo         6.0  1.7
looic       244.7 11.9

Pareto k diagnostic values:
                         Count  Pct 
(-Inf, 0.5]   (good)     49    98.0%
 (0.5, 0.7]   (ok)        0     0.0%
   (0.7, 1]   (bad)       1     2.0%
   (1, Inf)   (very bad)  0     0.0%
See help('pareto-k-diagnostic') for details.

# now fit a model without main factor
modelString = "
  data {
  int<lower=1> n;
  int<lower=1> nX;
  vector [n] y;
  matrix [n,nX] X;
  }
  parameters {
  vector[nX] beta;
  real<lower=0> sigma;
  }
  transformed parameters {
  vector[n] mu;

  mu = X*beta;
  }
  model {
  #Likelihood
  y~normal(mu,sigma);
  
  #Priors
  beta ~ normal(0,1000);
  sigma~cauchy(0,5);
  }
  generated quantities {
  vector[n] log_lik;
  
  for (i in 1:n) {
  log_lik[i] = normal_lpdf(y[i] | mu[i], sigma); 
  }
  }
  "

Xmat <- model.matrix(~1, data)
data.list <- with(data, list(y = y, X = Xmat, n = nrow(data), nX = ncol(Xmat)))
data.rstan.red <- stan(data = data.list, model_code = modelString, chains = 3,
    iter = 2000, warmup = 500, thin = 3)

SAMPLING FOR MODEL '3b057d3d81cbed2078ce678376a94574' NOW (CHAIN 1).

Gradient evaluation took 2.4e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.24 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  501 / 2000 [ 25%]  (Sampling)
Iteration:  700 / 2000 [ 35%]  (Sampling)
Iteration:  900 / 2000 [ 45%]  (Sampling)
Iteration: 1100 / 2000 [ 55%]  (Sampling)
Iteration: 1300 / 2000 [ 65%]  (Sampling)
Iteration: 1500 / 2000 [ 75%]  (Sampling)
Iteration: 1700 / 2000 [ 85%]  (Sampling)
Iteration: 1900 / 2000 [ 95%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 0.017015 seconds (Warm-up)
               0.038261 seconds (Sampling)
               0.055276 seconds (Total)


SAMPLING FOR MODEL '3b057d3d81cbed2078ce678376a94574' NOW (CHAIN 2).

Gradient evaluation took 7e-06 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  501 / 2000 [ 25%]  (Sampling)
Iteration:  700 / 2000 [ 35%]  (Sampling)
Iteration:  900 / 2000 [ 45%]  (Sampling)
Iteration: 1100 / 2000 [ 55%]  (Sampling)
Iteration: 1300 / 2000 [ 65%]  (Sampling)
Iteration: 1500 / 2000 [ 75%]  (Sampling)
Iteration: 1700 / 2000 [ 85%]  (Sampling)
Iteration: 1900 / 2000 [ 95%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 0.013452 seconds (Warm-up)
               0.026282 seconds (Sampling)
               0.039734 seconds (Total)


SAMPLING FOR MODEL '3b057d3d81cbed2078ce678376a94574' NOW (CHAIN 3).

Gradient evaluation took 7e-06 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  501 / 2000 [ 25%]  (Sampling)
Iteration:  700 / 2000 [ 35%]  (Sampling)
Iteration:  900 / 2000 [ 45%]  (Sampling)
Iteration: 1100 / 2000 [ 55%]  (Sampling)
Iteration: 1300 / 2000 [ 65%]  (Sampling)
Iteration: 1500 / 2000 [ 75%]  (Sampling)
Iteration: 1700 / 2000 [ 85%]  (Sampling)
Iteration: 1900 / 2000 [ 95%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 0.013892 seconds (Warm-up)
               0.026678 seconds (Sampling)
               0.04057 seconds (Total)

(reduced = loo(extract_log_lik(data.rstan.red)))

Computed from 1500 by 50 log-likelihood matrix

         Estimate  SE
elpd_loo   -178.4 4.1
p_loo         1.6 0.2
looic       356.9 8.2

All Pareto k estimates are good (k < 0.5)
See help('pareto-k-diagnostic') for details.

par(mfrow = 1:2, mar = c(5, 3.8, 1, 0) + 0.1, las = 3)
plot(full, label_points = TRUE)
plot(reduced, label_points = TRUE)

compare_models(full, reduced)

Error in discrete == discrete[1]: comparison of these types is not implemented

Matrix model (RSTANARM)

summary(data.rstanarm)

Model Info:

 function:  stan_glm
 family:    gaussian [identity]
 formula:   y ~ x
 algorithm: sampling
 priors:    see help('prior_summary')
 sample:    2250 (posterior sample size)
 num obs:   50

Estimates:
                mean   sd     2.5%   25%    50%    75%    97.5%
(Intercept)     40.4    0.8   38.7   39.8   40.4   40.9   42.0 
xB               5.4    1.2    3.0    4.6    5.3    6.1    7.7 
xC              14.2    1.2   11.8   13.4   14.2   15.0   16.6 
xD               0.0    1.2   -2.4   -0.7    0.0    0.8    2.3 
xE             -10.0    1.2  -12.3  -10.8  -10.0   -9.2   -7.6 
sigma            2.7    0.3    2.2    2.5    2.6    2.8    3.3 
mean_PPD        42.3    0.5   41.2   41.9   42.3   42.6   43.3 
log-posterior -133.2    1.9 -137.7 -134.2 -132.8 -131.8 -130.6 

Diagnostics:
              mcse Rhat n_eff
(Intercept)   0.0  1.0  1983 
xB            0.0  1.0  2003 
xC            0.0  1.0  2188 
xD            0.0  1.0  2077 
xE            0.0  1.0  2070 
sigma         0.0  1.0  1517 
mean_PPD      0.0  1.0  1777 
log-posterior 0.1  1.0  1088 

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

# OR
library(broom)
tidyMCMC(data.rstanarm$stanfit, conf.int = TRUE, conf.method = "HPDinterval")

           term      estimate std.error    conf.low   conf.high
1   (Intercept)   40.37013202 0.8369156   38.641206   41.904030
2            xB    5.35134132 1.1832942    3.008367    7.741613
3            xC   14.23032089 1.1924590   11.961883   16.741929
4            xD    0.01681653 1.1703009   -2.498183    2.092344
5            xE   -9.99518905 1.1917167  -12.239427   -7.563808
6         sigma    2.65469957 0.2833189    2.149518    3.229648
7      mean_PPD   42.27470974 0.5407395   41.269346   43.338946
8 log-posterior -133.18123981 1.9005441 -136.841953 -130.277664

• the mean of the first group (A) is 40.370132
• the mean of the second group (B) is 5.3513413 units greater than (A)
• the mean of the third group (C) is 14.2303209 units greater than (A)
• the mean of the forth group (D) is 0.0168165 units greater (i.e. less) than (A)
• the mean of the fifth group (E) is -9.995189 units greater (i.e. less) than (A)

While workers attempt to become comfortable with a new statistical framework, it is only natural that they like to evaluate and comprehend new structures and output alongside more familiar concepts. One way to facilitate this is via Bayesian p-values that are somewhat analogous to the frequentist p-values for investigating the hypothesis that a parameter is equal to zero.

## since values are less than zero
mcmcpvalue(as.matrix(data.rstanarm)[, "xB"])  # effect of (B-A)

[1] 0

mcmcpvalue(as.matrix(data.rstanarm)[, "xC"])  # effect of (C-A)

[1] 0

mcmcpvalue(as.matrix(data.rstanarm)[, "xD"])  # effect of (D-A)

[1] 0.9848889

mcmcpvalue(as.matrix(data.rstanarm)[, "xE"])  # effect of (E-A)

[1] 0

mcmcpvalue(as.matrix(data.rstanarm)[, 2:5])  # effect of (all groups)

[1] 0

There is evidence that the reponse differs between the groups. There is very little evidence that the response of group D differs from that of group A

library(loo)
(full = loo(data.rstanarm))

Computed from 2250 by 50 log-likelihood matrix

         Estimate   SE
elpd_loo   -122.2  5.8
p_loo         5.8  1.5
looic       244.3 11.5

Pareto k diagnostic values:
                         Count  Pct 
(-Inf, 0.5]   (good)     49    98.0%
 (0.5, 0.7]   (ok)        1     2.0%
   (0.7, 1]   (bad)       0     0.0%
   (1, Inf)   (very bad)  0     0.0%

All Pareto k estimates are ok (k < 0.7)
See help('pareto-k-diagnostic') for details.

data.rstanarm.red = update(data.rstanarm, . ~ 1)

Gradient evaluation took 2.4e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.24 seconds.
Adjust your expectations accordingly!



 Elapsed Time: 0.021685 seconds (Warm-up)
               0.044483 seconds (Sampling)
               0.066168 seconds (Total)


Gradient evaluation took 1e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
Adjust your expectations accordingly!



 Elapsed Time: 0.023106 seconds (Warm-up)
               0.085593 seconds (Sampling)
               0.108699 seconds (Total)


Gradient evaluation took 9e-06 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
Adjust your expectations accordingly!



 Elapsed Time: 0.02003 seconds (Warm-up)
               0.041798 seconds (Sampling)
               0.061828 seconds (Total)

(reduced = loo(data.rstanarm.red))

Computed from 2250 by 50 log-likelihood matrix

         Estimate  SE
elpd_loo   -178.4 4.0
p_loo         1.5 0.2
looic       356.7 8.0

All Pareto k estimates are good (k < 0.5)
See help('pareto-k-diagnostic') for details.

par(mfrow = 1:2, mar = c(5, 3.8, 1, 0) + 0.1, las = 3)
plot(full, label_points = TRUE)
plot(reduced, label_points = TRUE)

compare_models(full, reduced)

elpd_diff        se 
    -56.2       7.4 

Matrix model (BRMS)

summary(data.brms)

 Family: gaussian(identity) 
Formula: y ~ x 
   Data: data (Number of observations: 50) 
Samples: 3 chains, each with iter = 2000; warmup = 500; thin = 2; 
         total post-warmup samples = 2250
    ICs: LOO = NA; WAIC = NA; R2 = NA
 
Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept    40.37      0.85    38.73    41.99       1830    1
xB            5.38      1.20     2.97     7.77       1834    1
xC           14.22      1.19    11.90    16.69       1816    1
xD            0.03      1.18    -2.30     2.33       1822    1
xE           -9.96      1.20   -12.30    -7.65       1575    1

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     2.63      0.29     2.13     3.24       1967    1

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

# OR
library(broom)
tidyMCMC(data.brms$fit, conf.int = TRUE, conf.method = "HPDinterval")

         term   estimate std.error   conf.low conf.high
1 b_Intercept 40.3681829 0.8453842  38.733357 41.995095
2        b_xB  5.3835983 1.1963019   3.133155  7.926054
3        b_xC 14.2246055 1.1934948  12.002788 16.737062
4        b_xD  0.0322263 1.1761590  -2.398031  2.192691
5        b_xE -9.9603899 1.1972640 -12.385544 -7.745285
6       sigma  2.6319625 0.2902885   2.123581  3.222859

• the mean of the first group (A) is 40.3681829
• the mean of the second group (B) is 5.3835983 units greater than (A)
• the mean of the third group (C) is 14.2246055 units greater than (A)
• the mean of the forth group (D) is 0.0322263 units greater (i.e. less) than (A)
• the mean of the fifth group (E) is -9.9603899 units greater (i.e. less) than (A)

While workers attempt to become comfortable with a new statistical framework, it is only natural that they like to evaluate and comprehend new structures and output alongside more familiar concepts. One way to facilitate this is via Bayesian p-values that are somewhat analogous to the frequentist p-values for investigating the hypothesis that a parameter is equal to zero.

## since values are less than zero
mcmcpvalue(as.matrix(data.brms)[, "b_xB"])  # effect of (B-A)

[1] 0

mcmcpvalue(as.matrix(data.brms)[, "b_xC"])  # effect of (C-A)

[1] 0

mcmcpvalue(as.matrix(data.brms)[, "b_xD"])  # effect of (D-A)

[1] 0.9782222

mcmcpvalue(as.matrix(data.brms)[, "b_xE"])  # effect of (E-A)

[1] 0

mcmcpvalue(as.matrix(data.brms)[, 2:5])  # effect of (all groups)

[1] 0

There is evidence that the reponse differs between the groups. There is very little evidence that the response of group D differs from that of group A

library(loo)
(full = loo(data.brms))

  LOOIC    SE
 244.51 11.83

data.brms.red = update(data.brms, . ~ 1)

SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 1).

Gradient evaluation took 9e-06 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1001 / 2000 [ 50%]  (Sampling)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 0.015071 seconds (Warm-up)
               0.010419 seconds (Sampling)
               0.02549 seconds (Total)


SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).

Gradient evaluation took 2e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.2 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1001 / 2000 [ 50%]  (Sampling)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 0.014483 seconds (Warm-up)
               0.017408 seconds (Sampling)
               0.031891 seconds (Total)


SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).

Gradient evaluation took 3e-06 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1001 / 2000 [ 50%]  (Sampling)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 0.014374 seconds (Warm-up)
               0.012957 seconds (Sampling)
               0.027331 seconds (Total)

(reduced = loo(data.brms.red))

  LOOIC   SE
 356.74 8.09

par(mfrow = 1:2, mar = c(5, 3.8, 1, 0) + 0.1, las = 3)
plot(full, label_points = TRUE)
plot(reduced, label_points = TRUE)

compare_models(full, reduced)

Error in discrete == discrete[1]: comparison of these types is not implemented

A nice graphic is often a great accompaniment to a statistical analysis. Although there are no fixed assumptions associated with graphing (in contrast to statistical analyses), we often want the graphical summaries to reflect the associated statistical analyses. After all, the sample is just one perspective on the population(s). What we are more interested in is being able to estimate and depict likely population parameters/trends.

Thus, whilst we could easily provide a plot displaying the raw data along with simple measures of location and spread, arguably, we should use estimates that reflect the fitted model. In this case, it would be appropriate to plot the credibility interval associated with each group.

Matrix model (MCMCpack)

mcmc = data.mcmcpack
## Calculate the fitted values
newdata = rbind(data.frame(x = levels(data$x)))
Xmat = model.matrix(~x, newdata)
coefs = mcmc[, 1:5]
fit = coefs %*% t(Xmat)
newdata = newdata %>% cbind(tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval"))

ggplot(newdata, aes(y = estimate, x = x)) + geom_linerange(aes(ymin = conf.low,
    ymax = conf.high)) + geom_point() + scale_y_continuous("Y") + scale_x_discrete("X") +
    theme_classic()

As this is simple single factor ANOVA, we can simple add the raw data to this figure. For more complex designs with additional predictors, it is necessary to plot partial residuals.

## Calculate partial residuals fitted values
fdata = rdata = data
fMat = rMat = model.matrix(~x, fdata)
fit = as.vector(apply(coefs, 2, median) %*% t(fMat))
resid = as.vector(data$y - apply(coefs, 2, median) %*% t(rMat))
rdata = rdata %>% mutate(partial.resid = resid + fit)

ggplot(newdata, aes(y = estimate, x = as.numeric(x) - 0.1)) + geom_blank(aes(x = x)) +
    geom_point(data = rdata, aes(y = partial.resid, x = as.numeric(x) +
        0.1), color = "gray") + geom_linerange(aes(ymin = conf.low, ymax = conf.high)) +
    geom_point() + scale_y_continuous("Y") + scale_x_discrete("") + theme_classic()

Matrix model (JAGS)

mcmc = data.r2jags$BUGSoutput$sims.matrix
## Calculate the fitted values
newdata = rbind(data.frame(x = levels(data$x)))
Xmat = model.matrix(~x, newdata)
coefs = mcmc[, c("beta[1]", "beta[2]", "beta[3]", "beta[4]", "beta[5]")]
fit = coefs %*% t(Xmat)
newdata = newdata %>% cbind(tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval"))

ggplot(newdata, aes(y = estimate, x = x)) + geom_linerange(aes(ymin = conf.low,
    ymax = conf.high)) + geom_point() + scale_y_continuous("Y") + scale_x_discrete("X") +
    theme_classic()

As this is simple single factor ANOVA, we can simple add the raw data to this figure. For more complex designs with additional predictors, it is necessary to plot partial residuals.

## Calculate partial residuals fitted values
fdata = rdata = data
fMat = rMat = model.matrix(~x, fdata)
fit = as.vector(apply(coefs, 2, median) %*% t(fMat))
resid = as.vector(data$y - apply(coefs, 2, median) %*% t(rMat))
rdata = rdata %>% mutate(partial.resid = resid + fit)

ggplot(newdata, aes(y = estimate, x = as.numeric(x) - 0.1)) + geom_blank(aes(x = x)) +
    geom_point(data = rdata, aes(y = partial.resid, x = as.numeric(x) +
        0.1), color = "gray") + geom_linerange(aes(ymin = conf.low, ymax = conf.high)) +
    geom_point() + scale_y_continuous("Y") + scale_x_discrete("") + theme_classic()

Matrix model (RSTAN)

mcmc = as.matrix(data.rstan)
## Calculate the fitted values
newdata = rbind(data.frame(x = levels(data$x)))
Xmat = model.matrix(~x, newdata)
coefs = mcmc[, c("beta[1]", "beta[2]", "beta[3]", "beta[4]", "beta[5]")]
fit = coefs %*% t(Xmat)
newdata = newdata %>% cbind(tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval"))

ggplot(newdata, aes(y = estimate, x = x)) + geom_linerange(aes(ymin = conf.low,
    ymax = conf.high)) + geom_point() + scale_y_continuous("Y") + scale_x_discrete("X") +
    theme_classic()

As this is simple single factor ANOVA, we can simple add the raw data to this figure. For more complex designs with additional predictors, it is necessary to plot partial residuals.

## Calculate partial residuals fitted values
fdata = rdata = data
fMat = rMat = model.matrix(~x, fdata)
fit = as.vector(apply(coefs, 2, median) %*% t(fMat))
resid = as.vector(data$y - apply(coefs, 2, median) %*% t(rMat))
rdata = rdata %>% mutate(partial.resid = resid + fit)

ggplot(newdata, aes(y = estimate, x = as.numeric(x) - 0.1)) + geom_blank(aes(x = x)) +
    geom_point(data = rdata, aes(y = partial.resid, x = as.numeric(x) +
        0.1), color = "gray") + geom_linerange(aes(ymin = conf.low, ymax = conf.high)) +
    geom_point() + scale_y_continuous("Y") + scale_x_discrete("") + theme_classic()

Matrix model (RSTANARM)

## Calculate the fitted values
newdata = rbind(data.frame(x = levels(data$x)))
fit = posterior_linpred(data.rstanarm, newdata = newdata)
newdata = newdata %>% cbind(tidyMCMC(as.mcmc(fit), conf.int = TRUE,
    conf.method = "HPDinterval"))

ggplot(newdata, aes(y = estimate, x = x)) + geom_linerange(aes(ymin = conf.low,
    ymax = conf.high)) + geom_point() + scale_y_continuous("Y") +
    scale_x_discrete("X") + theme_classic()

As this is simple single factor ANOVA, we can simple add the raw data to this figure. For more complex designs with additional predictors, it is necessary to plot partial residuals.

## Calculate partial residuals
rdata = data
pp = posterior_linpred(data.rstanarm, newdata = rdata)
fit = as.vector(apply(pp, 2, median))
resid = resid(data.rstanarm)
rdata = rdata %>% mutate(partial.resid = resid + fit)

ggplot(newdata, aes(y = estimate, x = as.numeric(x) - 0.1)) +
    geom_blank(aes(x = x)) + geom_point(data = rdata, aes(y = partial.resid,
    x = as.numeric(x) + 0.1), color = "gray") + geom_linerange(aes(ymin = conf.low,
    ymax = conf.high)) + geom_point() + scale_y_continuous("Y") +
    scale_x_discrete("") + theme_classic()

Matrix model (BRMS)

Although we could calculated the fitted values via matrix multiplication of the coefficients and the model matrix (as for MCMCpack, RJAGS and RSTAN), for more complex models, it is more convenient to use the marginal_effects function that comes with brms.

plot(marginal_effects(data.brms), points = TRUE)

# OR
eff = plot(marginal_effects(data.brms), points = TRUE, plot = FALSE)
eff

$x

## Calculate the fitted values
newdata = rbind(data.frame(x = levels(data$x)))
fit = fitted(data.brms, newdata = newdata, summary = FALSE)
newdata = newdata %>% cbind(tidyMCMC(as.mcmc(fit), conf.int = TRUE,
    conf.method = "HPDinterval"))

ggplot(newdata, aes(y = estimate, x = x)) + geom_linerange(aes(ymin = conf.low,
    ymax = conf.high)) + geom_point() + scale_y_continuous("Y") +
    scale_x_discrete("X") + theme_classic()

As this is simple single factor ANOVA, we can simple add the raw data to this figure. For more complex designs with additional predictors, it is necessary to plot partial residuals.

## Calculate partial residuals
rdata = data
fit = fitted(data.brms, summary = TRUE)[, "Estimate"]
resid = resid(data.brms)[, "Estimate"]
rdata = rdata %>% mutate(partial.resid = resid + fit)

ggplot(newdata, aes(y = estimate, x = as.numeric(x) - 0.1)) +
    geom_blank(aes(x = x)) + geom_point(data = rdata, aes(y = partial.resid,
    x = as.numeric(x) + 0.1), color = "gray") + geom_linerange(aes(ymin = conf.low,
    ymax = conf.high)) + geom_point() + scale_y_continuous("Y") +
    scale_x_discrete("") + theme_classic()