Tutorial 7.5b - Analysis of Covariance (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 ($\mathbf{X}\boldsymbol{\beta}$). In this case, $\boldsymbol{\beta}$ represents the vector of $\beta$'s - the intercept associated with the first group, the (effects) differences between this intercept and the intercepts for each other group as well as the slope associated with the continuous covariate. $\mathbf{X}$ is the model matrix.

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_{ij} &\sim{} N(\mu_{ij}, \sigma)\\ \mu_{ij} &= \beta_0 + \mathbf{X}\boldsymbol{\beta}\\[1em] \beta_0 &\sim{} N(0,100)\\ \beta &\sim{} N(0,10)\\ \sigma &\sim{} cauchy(0,5)\\ \end{align} $$ Note, exploratory data analysis suggests that while the intercept (intercept of Group A) and categorical predictor effects (differences between intercepts of each of the Group and Group A's intercept) could be drawn from a similar distribution (with mean in the 10's and variances in the 100's), the slope (effect associated with Group A linear relationship) is likely to be an order of magnitude less. We might therefore be tempted to provide different priors for the intercept, categorical effects and slope effect. For a simple model such as this, it is unlikely to be necessary. However, for more complex models, where prior specification becomes more critical, separate priors would probably be necessary.

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(~A + B, 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.

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: 30
   Unobserved stochastic nodes: 5
   Total graph size: 232

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]   50.201   1.803  46.617  49.021  50.207  51.404  53.741 1.001  8400
beta[2]  -17.141   1.879 -20.832 -18.393 -17.149 -15.910 -13.402 1.001  8400
beta[3]  -22.917   1.882 -26.640 -24.161 -22.922 -21.682 -19.181 1.001  5900
beta[4]   -0.414   0.065  -0.543  -0.457  -0.414  -0.372  -0.285 1.001 15000
sigma      4.175   0.619   3.178   3.735   4.108   4.531   5.587 1.001  7200
deviance 169.232   3.548 164.517 166.636 168.502 171.075 178.013 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 = 6.3 and DIC = 175.5
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_{ij}&\sim{}N(\mu_{ij}, \sigma)\\ \mu_{ij} &= \mathbf{X}\boldsymbol{\beta}\\ \beta_0&\sim{}N(0,100)\\ \beta&\sim{}N(0,10)\\ \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,100);
  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(~A + B, 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 file4be61ee4e612.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 '3d2414c9dcf4b5e12be870eadd2c894a' 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.062651 seconds (Warm-up)
               0.061975 seconds (Sampling)
               0.124626 seconds (Total)


SAMPLING FOR MODEL '3d2414c9dcf4b5e12be870eadd2c894a' 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.053706 seconds (Warm-up)
               0.057327 seconds (Sampling)
               0.111033 seconds (Total)


SAMPLING FOR MODEL '3d2414c9dcf4b5e12be870eadd2c894a' 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.054613 seconds (Warm-up)
               0.064711 seconds (Sampling)
               0.119324 seconds (Total)

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

Inference for Stan model: 3d2414c9dcf4b5e12be870eadd2c894a.
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]  50.22    0.06 1.79  46.60  49.03  50.24  51.40  53.83  1027    1
beta[2] -17.12    0.06 1.92 -20.88 -18.31 -17.14 -15.96 -13.19  1210    1
beta[3] -22.88    0.05 1.93 -26.78 -24.11 -22.86 -21.58 -19.24  1407    1
beta[4]  -0.42    0.00 0.06  -0.54  -0.46  -0.42  -0.37  -0.29  1282    1
sigma     4.11    0.02 0.58   3.12   3.70   4.05   4.45   5.37  1353    1

Samples were drawn using NUTS(diag_e) at Fri Nov  3 15:05:11 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 ~ A + B, 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 5.5e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.55 seconds.
Adjust your expectations accordingly!



 Elapsed Time: 0.149605 seconds (Warm-up)
               0.181861 seconds (Sampling)
               0.331466 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.120154 seconds (Warm-up)
               0.198657 seconds (Sampling)
               0.318811 seconds (Total)


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



 Elapsed Time: 0.182133 seconds (Warm-up)
               0.174642 seconds (Sampling)
               0.356775 seconds (Total)

print(data.rstanarm)

stan_glm
 family:  gaussian [identity]
 formula: Y ~ A + B
------

Estimates:
            Median MAD_SD
(Intercept)  50.2    1.8 
AGroup B    -17.2    1.9 
AGroup C    -22.9    1.9 
B            -0.4    0.1 
sigma         4.1    0.6 

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

------
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)  50.2286902 1.9365801  46.4428738  54.1263008
2    AGroup B -17.2114440 1.9386491 -21.0066484 -13.4338089
3    AGroup C -22.9037572 1.9679361 -26.6757211 -19.0189639
4           B  -0.4138748 0.0681432  -0.5451839  -0.2852376
5       sigma   4.1990238 0.6317351   3.0581470   5.4328442

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 ~ A + B, 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 1.8e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.18 seconds.
Adjust your expectations accordingly!



 Elapsed Time: 0.036403 seconds (Warm-up)
               0.043672 seconds (Sampling)
               0.080075 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.033241 seconds (Warm-up)
               0.041044 seconds (Sampling)
               0.074285 seconds (Total)


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



 Elapsed Time: 0.036154 seconds (Warm-up)
               0.041079 seconds (Sampling)
               0.077233 seconds (Total)

print(data.brms)

 Family: gaussian(identity) 
Formula: Y ~ A + B 
   Data: data (Number of observations: 30) 
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    50.20      1.76    46.79    53.75       2177    1
AGroupB     -17.17      1.83   -20.81   -13.68       2158    1
AGroupC     -22.89      1.85   -26.43   -19.25       2127    1
B            -0.41      0.06    -0.54    -0.29       2129    1

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     4.11      0.59     3.14     5.41       1854    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  50.1970335 1.75540856  46.6017008  53.5440980
2   b_AGroupB -17.1728328 1.83135058 -20.6803108 -13.5582289
3   b_AGroupC -22.8915861 1.85449337 -26.3846325 -19.2178870
4         b_B  -0.4122714 0.06370931  -0.5416812  -0.2961189
5       sigma   4.1081576 0.58644550   3.1075609   5.3375706

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        3834  3746         1.020     
   AGroup B    2        3771  3746         1.010     
   AGroup C    2        3962  3746         1.060     
   B           2        3741  3746         0.999     
   sigma2      2        3771  3746         1.010     
  

  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)      AGroup B     AGroup C            B       sigma2
  Lag 0   1.000000000  1.000000e+00  1.000000000  1.000000000  1.000000000
  Lag 1  -0.003346420 -4.064355e-03 -0.002480710  0.002480390  0.132659644
  Lag 5   0.016500991  1.299476e-02 -0.003939996  0.001939800 -0.018030548
  Lag 10 -0.005233676 -6.125242e-06 -0.013564685 -0.006005238  0.006407515
  Lag 50  0.006330549 -9.343650e-03 -0.011855223 -0.002673821  0.007502186
  

  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]")
  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       36200 3746          9.66     
   beta[2]  20       38660 3746         10.30     
   beta[3]  20       38030 3746         10.20     
   beta[4]  20       36200 3746          9.66     
   deviance 20       37410 3746          9.99     
   sigma    20       38030 3746         10.20     
  
  
  [[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       39300 3746         10.50     
   beta[2]  20       38660 3746         10.30     
   beta[3]  20       39950 3746         10.70     
   beta[4]  20       37410 3746          9.99     
   deviance 20       39300 3746         10.50     
   sigma    20       36800 3746          9.82     
  
  
  [[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       38030 3746         10.20     
   beta[2]  20       36800 3746          9.82     
   beta[3]  20       37410 3746          9.99     
   beta[4]  20       38030 3746         10.20     
   deviance 20       36800 3746          9.82     
   sigma    20       36200 3746          9.66     
  

  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]      deviance         sigma
  Lag 0    1.0000000000  1.0000000000  1.000000000  1.000000000  1.0000000000  1.0000000000
  Lag 10  -0.0086789407  0.0002409218 -0.012834845 -0.014876564 -0.0155077394 -0.0009823887
  Lag 50  -0.0133637002 -0.0146591155 -0.022266022  0.007955582 -0.0131845968 -0.0066436459
  Lag 100  0.0007452402  0.0128511117 -0.001800657 -0.007853232  0.0008253937 -0.0101001901
  Lag 500 -0.0029682531  0.0011976990  0.001848702 -0.005150388 -0.0016034261 -0.0118212907
  

  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]
  Lag 0   1.000000000  1.0000000000  1.000000000
  Lag 1   0.064013692  0.0605531343  0.026397084
  Lag 5   0.007321517  0.0009337118  0.006726571
  Lag 10 -0.016376224 -0.0098527765 -0.009376572
  Lag 50 -0.017683148 -0.0290390559  0.010863075
  

  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)     AGroup B     AGroup C            B
  Lag 0   1.0000000000  1.000000000  1.000000000  1.000000000
  Lag 1   0.0290559273  0.001379705 -0.023247646 -0.018452103
  Lag 5  -0.0035230620 -0.003667234  0.004820692  0.008991159
  Lag 10  0.0180539117 -0.016292304  0.028751812  0.024889288
  Lag 50  0.0008045354 -0.013900306  0.008315425  0.009651236
  

  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

    library(bayesplot)
    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.6e-05 seconds
    1000 transitions using 10 leapfrog steps per transition would take 0.36 seconds.
    Adjust your expectations accordingly!
    
    
    
     Elapsed Time: 0.318686 seconds (Warm-up)
                   0.056968 seconds (Sampling)
                   0.375654 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.371641 seconds (Warm-up)
                   0.056263 seconds (Sampling)
                   0.427904 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(~A + B, newdata)
## get median parameter estimates
head(mcmc)

  (Intercept)  AGroup B  AGroup C          B   sigma2
1    50.18839 -18.23379 -22.87414 -0.4394640 14.11917
2    49.39556 -15.64877 -20.98086 -0.4528605 13.22004
3    49.44854 -16.36367 -23.66886 -0.3887300 15.57976
4    51.46891 -16.26642 -24.26313 -0.4709924 12.01597
5    50.55617 -17.54613 -22.44262 -0.4252965 16.35805
6    50.73340 -17.54078 -23.75846 -0.4598987 10.84731

coefs = apply(mcmc[, 1:4], 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(~A + B, newdata)
## get median parameter estimates
coefs = apply(mcmc[, 1:4], 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 = A))

ggplot(newdata) + geom_point(aes(y = resid, x = B))

And now for studentized residuals

mcmc = as.data.frame(data.mcmcpack)
# generate a model matrix
newdata = data
Xmat = model.matrix(~A + B, newdata)
## get median parameter estimates
coefs = apply(mcmc[, 1:4], 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(~A + B, data)
## get median parameter estimates
coefs = mcmc[, 1:4]
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(A = data$A, B = data$B, yRep) %>% gather(key = Sample,
    value = Value, -A, -B)
ggplot(newdata) + geom_violin(aes(y = Value, x = A, fill = "Model"),
    alpha = 0.5) + geom_violin(data = data, aes(y = Y, x = A,
    fill = "Obs"), alpha = 0.5) + geom_point(data = data, aes(y = Y,
    x = A), position = position_jitter(width = 0.1, height = 0),
    color = "black")

ggplot(newdata) + geom_violin(aes(y = Value, x = B, fill = "Model",
    group = B, color = A), alpha = 0.5) + geom_point(data = data,
    aes(y = Y, x = B, group = B, color = A))

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|^A|^B|sigma")

mcmc_areas(as.matrix(data.mcmcpack), regex_pars = "Intercept|^A|^B|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(~A + B, newdata)
## get median parameter estimates
coefs = apply(mcmc[, 1:4], 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(~A + B, newdata)
## get median parameter estimates
coefs = apply(mcmc[, 1:4], 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 = A))

ggplot(newdata) + geom_point(aes(y = resid, x = B))

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(~A + B, newdata)
## get median parameter estimates
coefs = apply(mcmc[, 1:4], 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(~A + B, data)
## get median parameter estimates
coefs = mcmc[, 1:4]
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(A = data$A, B = data$B, yRep) %>% gather(key = Sample,
    value = Value, -A, -B)
ggplot(newdata) + geom_violin(aes(y = Value, x = A, fill = "Model"),
    alpha = 0.5) + geom_violin(data = data, aes(y = Y, x = A,
    fill = "Obs"), alpha = 0.5) + geom_point(data = data, aes(y = Y,
    x = A), position = position_jitter(width = 0.1, height = 0),
    color = "black")

ggplot(newdata) + geom_violin(aes(y = Value, x = B, fill = "Model",
    group = B, color = A), alpha = 0.5) + geom_point(data = data,
    aes(y = Y, x = B, group = B, color = A))

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(~A + B, newdata)
## get median parameter estimates
coefs = apply(mcmc[, 1:4], 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(~A + B, newdata)
## get median parameter estimates
coefs = apply(mcmc[, 1:4], 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 = A))

ggplot(newdata) + geom_point(aes(y = resid, x = B))

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(~A + B, newdata)
## get median parameter estimates
coefs = apply(mcmc[, 1:4], 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(~A + B, data)
## get median parameter estimates
coefs = mcmc[, 1:4]
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(A = data$A, B = data$B, yRep) %>% gather(key = Sample,
    value = Value, -A, -B)
ggplot(newdata) + geom_violin(aes(y = Value, x = A, fill = "Model"),
    alpha = 0.5) + geom_violin(data = data, aes(y = Y, x = A,
    fill = "Obs"), alpha = 0.5) + geom_point(data = data, aes(y = Y,
    x = A), position = position_jitter(width = 0.1, height = 0),
    color = "black")

ggplot(newdata) + geom_violin(aes(y = Value, x = B, fill = "Model",
    group = B, color = A), alpha = 0.5) + geom_point(data = data,
    aes(y = Y, x = B, group = B, color = A))

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 = A))

ggplot(dat) + geom_point(aes(y = resid, x = B))

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:-B)
ggplot(newdata) + geom_violin(aes(y = Value, x = A, fill = "Model"),
    alpha = 0.5) + geom_violin(data = data, aes(y = Y, x = A,
    fill = "Obs"), alpha = 0.5) + geom_point(data = data, aes(y = Y,
    x = A), position = position_jitter(width = 0.1, height = 0),
    color = "black")

Error: Aesthetics must be either length 1 or the same as the data (67530): x, y, fill

ggplot(newdata) + geom_violin(aes(y = Value, x = B, fill = "Model",
    group = B, color = A), alpha = 0.5) + geom_point(data = data,
    aes(y = Y, x = B, group = B, color = A))

Error: Aesthetics must be either length 1 or the same as the data (67530): x, y, fill, group, colour

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|^A|^B|sigma")

mcmc_areas(as.matrix(data.rstanarm), regex_pars = "Intercept|^A|^B|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 = A))

ggplot(dat) + geom_point(aes(y = resid, x = B))

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:-B)
ggplot(newdata) + geom_violin(aes(y = Value, x = A, fill = "Model"),
    alpha = 0.5) + geom_violin(data = data, aes(y = Y, x = A,
    fill = "Obs"), alpha = 0.5) + geom_point(data = data, aes(y = Y,
    x = A), position = position_jitter(width = 0.1, height = 0),
    color = "black")

Error: Aesthetics must be either length 1 or the same as the data (67530): x, y, fill

ggplot(newdata) + geom_violin(aes(y = Value, x = B, fill = "Model",
    group = B, color = A), alpha = 0.5) + geom_point(data = data,
    aes(y = Y, x = B, group = B, color = A))

Error: Aesthetics must be either length 1 or the same as the data (67530): x, y, fill, group, colour

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)  50.2015 1.7846 0.017846      0.0178458
AGroup B    -17.1433 1.8782 0.018782      0.0187821
AGroup C    -22.9133 1.8682 0.018682      0.0186817
B            -0.4138 0.0637 0.000637      0.0006543
sigma2       17.0079 5.0889 0.050889      0.0581569

2. Quantiles for each variable:

                2.5%      25%      50%      75%    97.5%
(Intercept)  46.7714  49.0274  50.2019  51.3437  53.7977
AGroup B    -20.8797 -18.3726 -17.1445 -15.8884 -13.4552
AGroup C    -26.6231 -24.1304 -22.9030 -21.6955 -19.1900
B            -0.5397  -0.4556  -0.4134  -0.3714  -0.2885
sigma2        9.7749  13.4022  16.1334  19.6823  29.3438

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

         term    estimate  std.error    conf.low   conf.high
1 (Intercept)  50.2014506 1.78458442  46.7660798  53.7766054
2    AGroup B -17.1433463 1.87821127 -20.7503649 -13.3485818
3    AGroup C -22.9132573 1.86817189 -26.6020778 -19.1852142
4           B  -0.4137879 0.06369625  -0.5369301  -0.2860994
5      sigma2  17.0078626 5.08890812   8.9411967  27.1577316

• the intercept of the first group (Group A) is 50.2014506
• the mean of the second group (Group B) is -17.1433463 units greater than (A)
• the mean of the third group (Group C) is -22.9132573 units greater than (A)
• a one unit increase in B in Group A is associated with a -0.4137879 units increase in Y

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 = 0)

[1] 0

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

[1] 0

mcmcpvalue(data.mcmcpack[, 4])  # effect of (slope = 0)

[1] 0

mcmcpvalue(data.mcmcpack[, 2:4])  # effect of (model)

[1] 0

There is evidence that the reponse differs between the groups.

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]   50.201   1.803  46.617  49.021  50.207  51.404  53.741 1.001  8400
beta[2]  -17.141   1.879 -20.832 -18.393 -17.149 -15.910 -13.402 1.001  8400
beta[3]  -22.917   1.882 -26.640 -24.161 -22.922 -21.682 -19.181 1.001  5900
beta[4]   -0.414   0.065  -0.543  -0.457  -0.414  -0.372  -0.285 1.001 15000
sigma      4.175   0.619   3.178   3.735   4.108   4.531   5.587 1.001  7200
deviance 169.232   3.548 164.517 166.636 168.502 171.075 178.013 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 = 6.3 and DIC = 175.5
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]  50.2007407 1.80285405  46.5994039  53.6923687
2  beta[2] -17.1407725 1.87934256 -20.8989661 -13.5038133
3  beta[3] -22.9165431 1.88243441 -26.5907610 -19.1480417
4  beta[4]  -0.4139048 0.06529536  -0.5434197  -0.2853424
5 deviance 169.2319386 3.54753309 163.9376050 176.1387130
6    sigma   4.1749569 0.61875758   3.0490746   5.3897463

• the intercept of the first group (Group A) is 50.2007407
• the mean of the second group (Group B) is -17.1407725 units greater than (A)
• the mean of the third group (Group C) is -22.9165431 units greater than (A)
• a one unit increase in B in Group A is associated with a -0.4139048 units increase in Y

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 = 0)

[1] 0

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

[1] 0

mcmcpvalue(data.r2jags$BUGSoutput$sims.matrix[, "beta[4]"])  # effect of (slope = 0)

[1] 0

mcmcpvalue(data.r2jags$BUGSoutput$sims.matrix[, 2:4])  # effect of (model)

[1] 0

There is evidence that the reponse differs between the groups.

Matrix model (RSTAN)

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

Inference for Stan model: 3d2414c9dcf4b5e12be870eadd2c894a.
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]  50.22    0.06 1.79  46.60  49.03  50.24  51.40  53.83  1027    1
beta[2] -17.12    0.06 1.92 -20.88 -18.31 -17.14 -15.96 -13.19  1210    1
beta[3] -22.88    0.05 1.93 -26.78 -24.11 -22.86 -21.58 -19.24  1407    1
beta[4]  -0.42    0.00 0.06  -0.54  -0.46  -0.42  -0.37  -0.29  1282    1
sigma     4.11    0.02 0.58   3.12   3.70   4.05   4.45   5.37  1353    1

Samples were drawn using NUTS(diag_e) at Fri Nov  3 15:05:11 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]  50.2161150 1.78767398  46.4168478  53.5737947
2 beta[2] -17.1206004 1.91600172 -20.7216779 -13.0233814
3 beta[3] -22.8794576 1.92793878 -26.8184381 -19.3130173
4 beta[4]  -0.4157079 0.06337683  -0.5361929  -0.2891039
5   sigma   4.1094504 0.58022288   3.0779992   5.3039276

• the intercept of the first group (Group A) is 50.216115
• the mean of the second group (Group B) is -17.1206004 units greater than (A)
• the mean of the third group (Group C) is -22.8794576 units greater than (A)
• a one unit increase in B in Group A is associated with a -0.4157079 units increase in Y

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 = 0)

[1] 0

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

[1] 0

mcmcpvalue(as.matrix(data.rstan)[, "beta[4]"])  # effect of (slope = 0)

[1] 0

mcmcpvalue(as.matrix(data.rstan)[, 2:4])  # effect of (model)

[1] 0

There is evidence that the reponse differs between the groups.

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

Computed from 1500 by 30 log-likelihood matrix

         Estimate  SE
elpd_loo    -87.3 3.5
p_loo         4.6 1.0
looic       174.5 7.1

Pareto k diagnostic values:
                         Count  Pct 
(-Inf, 0.5]   (good)     29    96.7%
 (0.5, 0.7]   (ok)        1     3.3%
   (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.

# 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)

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 file12361ebf982d.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.2e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.12 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.018509 seconds (Warm-up)
               0.041875 seconds (Sampling)
               0.060384 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.017436 seconds (Warm-up)
               0.046029 seconds (Sampling)
               0.063465 seconds (Total)


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

Gradient evaluation took 6e-06 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.06 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.012316 seconds (Warm-up)
               0.026582 seconds (Sampling)
               0.038898 seconds (Total)

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

Computed from 1500 by 30 log-likelihood matrix

         Estimate  SE
elpd_loo   -116.0 3.2
p_loo         1.6 0.4
looic       232.0 6.4

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 ~ A + B
 algorithm: sampling
 priors:    see help('prior_summary')
 sample:    2250 (posterior sample size)
 num obs:   30

Estimates:
                mean   sd     2.5%   25%    50%    75%    97.5%
(Intercept)     50.2    1.9   46.4   49.0   50.2   51.5   54.1 
AGroup B       -17.2    1.9  -21.0  -18.5  -17.2  -15.9  -13.3 
AGroup C       -22.9    2.0  -26.6  -24.2  -22.9  -21.6  -19.0 
B               -0.4    0.1   -0.5   -0.5   -0.4   -0.4   -0.3 
sigma            4.2    0.6    3.2    3.8    4.1    4.5    5.6 
mean_PPD        30.0    1.1   27.8   29.3   30.0   30.7   32.2 
log-posterior  -98.3    1.8 -103.1  -99.1  -97.9  -96.9  -95.9 

Diagnostics:
              mcse Rhat n_eff
(Intercept)   0.0  1.0  1980 
AGroup B      0.0  1.0  2145 
AGroup C      0.0  1.0  2250 
B             0.0  1.0  2250 
sigma         0.0  1.0  1200 
mean_PPD      0.0  1.0  1720 
log-posterior 0.1  1.0   960 

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)  50.2286902 1.9365801   46.4428738  54.1263008
2      AGroup B -17.2114440 1.9386491  -21.0066484 -13.4338089
3      AGroup C -22.9037572 1.9679361  -26.6757211 -19.0189639
4             B  -0.4138748 0.0681432   -0.5451839  -0.2852376
5         sigma   4.1990238 0.6317351    3.0581470   5.4328442
6      mean_PPD  30.0266842 1.1171901   27.8860300  32.2422092
7 log-posterior -98.2605466 1.8364590 -101.7802647 -95.7415132

• the intercept of the first group (Group A) is 50.2286902
• the mean of the second group (Group B) is -17.211444 units greater than (A)
• the mean of the third group (Group C) is -22.9037572 units greater than (A)
• a one unit increase in B in Group A is associated with a -0.4138748 units increase in Y

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
head(as.matrix(data.rstanarm))

          parameters
iterations (Intercept)  AGroup B  AGroup C          B    sigma
      [1,]    49.25161 -18.13674 -20.71374 -0.3963549 4.384061
      [2,]    50.28149 -20.18446 -23.60224 -0.3823662 4.601748
      [3,]    47.22865 -13.74348 -20.42610 -0.3439257 3.559008
      [4,]    50.36383 -18.72387 -23.75338 -0.3776768 4.211066
      [5,]    49.33028 -18.88190 -20.94786 -0.3634085 3.830714
      [6,]    49.08087 -14.40159 -22.87480 -0.4144611 3.793997

mcmcpvalue(as.matrix(data.rstanarm)[, "AGroup B"])  # effect of (B-A = 0)

[1] 0

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

[1] 0

mcmcpvalue(as.matrix(data.rstanarm)[, "B"])  # effect of (slope = 0)

[1] 0

mcmcpvalue(as.matrix(data.rstanarm)[, 2:4])  # effect of (model)

[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 30 log-likelihood matrix

         Estimate  SE
elpd_loo    -87.2 3.3
p_loo         4.3 0.9
looic       174.4 6.7

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

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

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



 Elapsed Time: 0.027383 seconds (Warm-up)
               0.039991 seconds (Sampling)
               0.067374 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.0218 seconds (Warm-up)
               0.039113 seconds (Sampling)
               0.060913 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.02501 seconds (Warm-up)
               0.043676 seconds (Sampling)
               0.068686 seconds (Total)

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

Computed from 2250 by 30 log-likelihood matrix

         Estimate  SE
elpd_loo   -116.0 3.1
p_loo         1.5 0.4
looic       232.0 6.1

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 
    -28.8       4.5 

Matrix model (BRMS)

summary(data.brms)

 Family: gaussian(identity) 
Formula: Y ~ A + B 
   Data: data (Number of observations: 30) 
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    50.20      1.76    46.79    53.75       2177    1
AGroupB     -17.17      1.83   -20.81   -13.68       2158    1
AGroupC     -22.89      1.85   -26.43   -19.25       2127    1
B            -0.41      0.06    -0.54    -0.29       2129    1

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     4.11      0.59     3.14     5.41       1854    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  50.1970335 1.75540856  46.6017008  53.5440980
2   b_AGroupB -17.1728328 1.83135058 -20.6803108 -13.5582289
3   b_AGroupC -22.8915861 1.85449337 -26.3846325 -19.2178870
4         b_B  -0.4122714 0.06370931  -0.5416812  -0.2961189
5       sigma   4.1081576 0.58644550   3.1075609   5.3375706

• the intercept of the first group (Group A) is 50.1970335
• the mean of the second group (Group B) is -17.1728328 units greater than (A)
• the mean of the third group (Group C) is -22.8915861 units greater than (A)
• a one unit increase in B in Group A is associated with a -0.4122714 units increase in Y

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
head(as.matrix(data.brms))

          parameters
iterations b_Intercept b_AGroupB b_AGroupC        b_B    sigma      lp__
      [1,]    46.14429 -17.89799 -23.20626 -0.1747196 4.234569 -112.8789
      [2,]    50.46466 -17.00736 -21.94628 -0.4189975 3.964478 -105.7635
      [3,]    49.51423 -16.74688 -22.42555 -0.4566537 3.784797 -106.7988
      [4,]    52.03476 -20.63035 -26.45616 -0.3795796 4.637355 -108.6522
      [5,]    49.58570 -18.21573 -23.70766 -0.4119922 4.153369 -106.9716
      [6,]    51.87136 -19.76110 -26.12697 -0.4639165 3.485067 -109.3362

mcmcpvalue(as.matrix(data.brms)[, "b_AGroupB"])  # effect of (B-A = 0)

[1] 0

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

[1] 0

mcmcpvalue(as.matrix(data.brms)[, "b_B"])  # effect of (slope = 0)

[1] 0

mcmcpvalue(as.matrix(data.brms)[, 2:4])  # effect of (model)

[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
 174.34 7.05

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

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

Gradient evaluation took 1e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.1 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.019211 seconds (Warm-up)
               0.014546 seconds (Sampling)
               0.033757 seconds (Total)


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

Gradient evaluation took 4e-06 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.04 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.015666 seconds (Warm-up)
               0.017803 seconds (Sampling)
               0.033469 seconds (Total)


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

Gradient evaluation took 4e-06 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.04 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.016449 seconds (Warm-up)
               0.016509 seconds (Sampling)
               0.032958 seconds (Total)

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

  LOOIC   SE
 232.15 6.38

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 = expand.grid(A = levels(data$A), B = seq(min(data$B), max(data$B),
    len = 100))
Xmat = model.matrix(~A + B, newdata)
coefs = mcmc[, 1:4]
fit = coefs %*% t(Xmat)
newdata = newdata %>% cbind(tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval"))

ggplot(newdata, aes(y = estimate, x = B, fill = A)) + geom_ribbon(aes(ymin = conf.low,
    ymax = conf.high), alpha = 0.2) + geom_line() + scale_y_continuous("Y") +
    scale_x_continuous("B") + 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(~A + B, 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 = B, fill = A)) + geom_point(data = rdata,
    aes(y = partial.resid, x = B, color = A)) + geom_ribbon(aes(ymin = conf.low,
    ymax = conf.high), alpha = 0.2) + geom_line() + scale_y_continuous("Y") +
    scale_x_continuous("B") + theme_classic()

Matrix model (JAGS)

mcmc = data.r2jags$BUGSoutput$sims.matrix
## Calculate the fitted values
newdata = expand.grid(A = levels(data$A), B = seq(min(data$B), max(data$B),
    len = 100))
Xmat = model.matrix(~A + B, newdata)
coefs = mcmc[, c("beta[1]", "beta[2]", "beta[3]", "beta[4]")]
fit = coefs %*% t(Xmat)
newdata = newdata %>% cbind(tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval"))

ggplot(newdata, aes(y = estimate, x = B, fill = A)) + geom_ribbon(aes(ymin = conf.low,
    ymax = conf.high), alpha = 0.2) + geom_line() + scale_y_continuous("Y") +
    scale_x_continuous("B") + 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(~A + B, 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 = B, fill = A)) + geom_point(data = rdata,
    aes(y = partial.resid, x = B, color = A)) + geom_ribbon(aes(ymin = conf.low,
    ymax = conf.high), alpha = 0.2) + geom_line() + scale_y_continuous("Y") +
    scale_x_continuous("B") + theme_classic()

Matrix model (RSTAN)

mcmc = as.matrix(data.rstan)
## Calculate the fitted values
newdata = expand.grid(A = levels(data$A), B = seq(min(data$B), max(data$B),
    len = 100))
Xmat = model.matrix(~A + B, newdata)
coefs = mcmc[, c("beta[1]", "beta[2]", "beta[3]", "beta[4]")]
fit = coefs %*% t(Xmat)
newdata = newdata %>% cbind(tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval"))

ggplot(newdata, aes(y = estimate, x = B, fill = A)) + geom_ribbon(aes(ymin = conf.low,
    ymax = conf.high), alpha = 0.2) + geom_line() + scale_y_continuous("Y") +
    scale_x_continuous("B") + 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(~A + B, 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 = B, fill = A)) + geom_point(data = rdata,
    aes(y = partial.resid, x = B, color = A)) + geom_ribbon(aes(ymin = conf.low,
    ymax = conf.high), alpha = 0.2) + geom_line() + scale_y_continuous("Y") +
    scale_x_continuous("B") + theme_classic()

Matrix model (RSTANARM)

## Calculate the fitted values
newdata = expand.grid(A = levels(data$A), B = seq(min(data$B),
    max(data$B), len = 100))
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 = B, fill = A)) + geom_ribbon(aes(ymin = conf.low,
    ymax = conf.high), alpha = 0.2) + geom_line() + scale_y_continuous("Y") +
    scale_x_continuous("B") + 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 = B, fill = A)) + geom_point(data = rdata,
    aes(y = partial.resid, x = B, color = A)) + geom_ribbon(aes(ymin = conf.low,
    ymax = conf.high), alpha = 0.2) + geom_line() + scale_y_continuous("Y") +
    scale_x_continuous("B") + 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

$A

$B

## Calculate the fitted values
newdata = expand.grid(A = levels(data$A), B = seq(min(data$B),
    max(data$B), len = 100))
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 = B, fill = A)) + geom_ribbon(aes(ymin = conf.low,
    ymax = conf.high), alpha = 0.2) + geom_line() + scale_y_continuous("Y") +
    scale_x_continuous("B") + 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 = B, fill = A)) + geom_point(data = rdata,
    aes(y = partial.resid, x = B, color = A)) + geom_ribbon(aes(ymin = conf.low,
    ymax = conf.high), alpha = 0.2) + geom_line() + scale_y_continuous("Y") +
    scale_x_continuous("B") + theme_classic()

In order to allow direct comparison to the frequentist equivalents, I will explore the same set of planned and "Tukey's" test comparisons described here. For the "planned comparison" we defined two contrasts:

1. group B vs group C
2. group A vs the average of groups B and C

Lets start by comparing each group to each other group in a pairwise manner. Arguably the most elegant way to do this is to generate a Tukey's contrast matrix. This is a model matrix specific to comparing each group to each other group. Again, since the lines are parallel, it does not really matter what level of B we estimate these efffects at - so lets use the mean B.

mcmc = data.mcmcpack
coefs <- as.matrix(mcmc)[, 1:4]
newdata = expand.grid(A = levels(data$A), B = mean(data$B))
# A Tukeys contrast matrix
library(multcomp)
tuk.mat <- contrMat(n = table(newdata$A), type = "Tukey")
tuk.mat

	 Multiple Comparisons of Means: Tukey Contrasts

                  Group A Group B Group C
Group B - Group A      -1       1       0
Group C - Group A      -1       0       1
Group C - Group B       0      -1       1

Xmat <- model.matrix(~A + B, data = newdata)
Xmat

  (Intercept) AGroup B AGroup C        B
1           1        0        0 16.50175
2           1        1        0 16.50175
3           1        0        1 16.50175
attr(,"assign")
[1] 0 1 1 2
attr(,"contrasts")
attr(,"contrasts")$A
[1] "contr.treatment"

pairwise.mat <- tuk.mat %*% Xmat
pairwise.mat

                  (Intercept) AGroup B AGroup C B
Group B - Group A           0        1        0 0
Group C - Group A           0        0        1 0
Group C - Group B           0       -1        1 0

mcmc_areas(coefs %*% t(pairwise.mat))

(comps = tidyMCMC(coefs %*% t(pairwise.mat), conf.int = TRUE, conf.method = "HPDinterval"))

               term   estimate std.error   conf.low  conf.high
1 Group B - Group A -17.143346  1.878211 -20.750365 -13.348582
2 Group C - Group A -22.913257  1.868172 -26.602078 -19.185214
3 Group C - Group B  -5.769911  1.821653  -9.552934  -2.271914

ggplot(comps, aes(y = estimate, x = term)) + geom_pointrange(aes(ymin = conf.low,
    ymax = conf.high)) + geom_hline(yintercept = 0, linetype = "dashed") +
    scale_y_continuous("Effect size") + scale_x_discrete("") + coord_flip() +
    theme_classic()

With a couple of modifications, we could also express this as percentage changes. A percentage change represents the change (difference between groups) divided by one of the groups (determined by which group you want to express the percentage change to). Hence, we generate an additional mcmc matrix that represents the cell means for the divisor group (group we want to express change relative to).

Since the tuk.mat defines comparisons as -1 and 1 pairs, if we simply replace all the -1 with 0, the eventual matrix multiplication will result in estimates of the divisor cell means instread of the difference. We can then divide the original mcmc matrix above with this new divisor mcmc matrix to yeild a mcmc matrix of percentage change.

# Modify the tuk.mat to replace -1 with 0.  This will allow us to get a
# mcmc matrix of ..
tuk.mat[tuk.mat == -1] = 0
comp.mat <- tuk.mat %*% Xmat
comp.mat

                  (Intercept) AGroup B AGroup C        B
Group B - Group A           1        1        0 16.50175
Group C - Group A           1        0        1 16.50175
Group C - Group B           1        0        1 16.50175

comp.mcmc = 100 * (coefs %*% t(pairwise.mat))/coefs %*% t(comp.mat)
(comps = tidyMCMC(comp.mcmc, conf.int = TRUE, conf.method = "HPDinterval"))

               term   estimate std.error   conf.low  conf.high
1 Group B - Group A  -65.78138  9.862131  -84.99925 -46.617429
2 Group C - Group A -112.84960 15.200206 -144.80969 -84.547452
3 Group C - Group B  -28.70578 10.358636  -50.05632  -9.004512

ggplot(comps, aes(y = estimate, x = term)) + geom_pointrange(aes(ymin = conf.low,
    ymax = conf.high)) + geom_hline(yintercept = 0, linetype = "dashed") +
    scale_y_continuous("Effect size (%)") + scale_x_discrete("") + coord_flip() +
    theme_classic()

And now for the specific planned comparisons (Group B vs Group C as well as Group A vs the average of Groups B and C). This is achieved by generating our own contrast matrix (defining the contributions of each group to each contrast).

c.mat = rbind(c(0, 1, -1), c(1/2, -1/3, -1/3))
c.mat

     [,1]       [,2]       [,3]
[1,]  0.0  1.0000000 -1.0000000
[2,]  0.5 -0.3333333 -0.3333333

mcmc = data.mcmcpack
coefs <- as.matrix(mcmc)[, 1:4]
newdata <- data.frame(A = levels(data$A), B = mean(data$B))
Xmat <- model.matrix(~A + B, data = newdata)
c.mat = c.mat %*% Xmat
c.mat

     (Intercept)   AGroup B   AGroup C         B
[1,]   0.0000000  1.0000000 -1.0000000  0.000000
[2,]  -0.1666667 -0.3333333 -0.3333333 -2.750292

(comps = tidyMCMC(as.mcmc(coefs %*% t(c.mat)), conf.int = TRUE, conf.method = "HPDinterval"))

  term estimate std.error conf.low conf.high
1 var1 5.769911 1.8216528 2.271914  9.552934
2 var2 6.123330 0.9187204 4.344498  7.965326

Lets start by comparing each group to each other group in a pairwise manner. Arguably the most elegant way to do this is to generate a Tukey's contrast matrix. This is a model matrix specific to comparing each group to each other group. Again, since the lines are parallel, it does not really matter what level of B we estimate these efffects at - so lets use the mean B.

mcmc = data.r2jags$BUGSoutput$sims.matrix
coefs <- as.matrix(mcmc)[, 1:4]
newdata <- data.frame(A = levels(data$A), B = mean(data$B))
# A Tukeys contrast matrix
library(multcomp)
tuk.mat <- contrMat(n = table(newdata$A), type = "Tukey")
Xmat <- model.matrix(~A + B, data = newdata)
pairwise.mat <- tuk.mat %*% Xmat
pairwise.mat

                  (Intercept) AGroup B AGroup C B
Group B - Group A           0        1        0 0
Group C - Group A           0        0        1 0
Group C - Group B           0       -1        1 0

mcmc_areas(coefs %*% t(pairwise.mat))

(comps = tidyMCMC(coefs %*% t(pairwise.mat), conf.int = TRUE, conf.method = "HPDinterval"))

               term   estimate std.error   conf.low  conf.high
1 Group B - Group A -17.140773  1.879343 -20.898966 -13.503813
2 Group C - Group A -22.916543  1.882434 -26.590761 -19.148042
3 Group C - Group B  -5.775771  1.880763  -9.402033  -1.979244

ggplot(comps, aes(y = estimate, x = term)) + geom_pointrange(aes(ymin = conf.low,
    ymax = conf.high)) + geom_hline(yintercept = 0, linetype = "dashed") +
    scale_y_continuous("Effect size") + scale_x_discrete("") + coord_flip() +
    theme_classic()

With a couple of modifications, we could also express this as percentage changes. A percentage change represents the change (difference between groups) divided by one of the groups (determined by which group you want to express the percentage change to). Hence, we generate an additional mcmc matrix that represents the cell means for the divisor group (group we want to express change relative to).

Since the tuk.mat defines comparisons as -1 and 1 pairs, if we simply replace all the -1 with 0, the eventual matrix multiplication will result in estimates of the divisor cell means instread of the difference. We can then divide the original mcmc matrix above with this new divisor mcmc matrix to yeild a mcmc matrix of percentage change.

# Modify the tuk.mat to replace -1 with 0.  This will allow us to get a
# mcmc matrix of ..
tuk.mat[tuk.mat == -1] = 0
comp.mat <- tuk.mat %*% Xmat
comp.mat

                  (Intercept) AGroup B AGroup C        B
Group B - Group A           1        1        0 16.50175
Group C - Group A           1        0        1 16.50175
Group C - Group B           1        0        1 16.50175

comp.mcmc = 100 * (coefs %*% t(pairwise.mat))/coefs %*% t(comp.mat)
(comps = tidyMCMC(comp.mcmc, conf.int = TRUE, conf.method = "HPDinterval"))

               term   estimate std.error   conf.low  conf.high
1 Group B - Group A  -65.77396  9.861029  -85.38307 -46.803881
2 Group C - Group A -112.95081 15.533036 -144.44128 -83.512153
3 Group C - Group B  -28.79007 10.740187  -50.07490  -8.241323

ggplot(comps, aes(y = estimate, x = term)) + geom_pointrange(aes(ymin = conf.low,
    ymax = conf.high)) + geom_hline(yintercept = 0, linetype = "dashed") +
    scale_y_continuous("Effect size (%)") + scale_x_discrete("") + coord_flip() +
    theme_classic()

And now for the specific planned comparisons (Group B vs Group C as well as Group A vs the average of Groups B and C). This is achieved by generating our own contrast matrix (defining the contributions of each group to each contrast).

c.mat = rbind(c(0, 1, -1), c(1/2, -1/3, -1/3))
c.mat

     [,1]       [,2]       [,3]
[1,]  0.0  1.0000000 -1.0000000
[2,]  0.5 -0.3333333 -0.3333333

mcmc = data.r2jags$BUGSoutput$sims.matrix
coefs <- as.matrix(mcmc)[, 1:4]
newdata <- data.frame(A = levels(data$A), B = mean(data$B))
Xmat <- model.matrix(~A + B, data = newdata)
c.mat = c.mat %*% Xmat
c.mat

     (Intercept)   AGroup B   AGroup C         B
[1,]   0.0000000  1.0000000 -1.0000000  0.000000
[2,]  -0.1666667 -0.3333333 -0.3333333 -2.750292

(comps = tidyMCMC(as.mcmc(coefs %*% t(c.mat)), conf.int = TRUE, conf.method = "HPDinterval"))

  term estimate std.error conf.low conf.high
1 var1 5.775771  1.880763 1.979244  9.402033
2 var2 6.124007  0.913995 4.278228  7.879558

Lets start by comparing each group to each other group in a pairwise manner. Arguably the most elegant way to do this is to generate a Tukey's contrast matrix. This is a model matrix specific to comparing each group to each other group. Again, since the lines are parallel, it does not really matter what level of B we estimate these efffects at - so lets use the mean B.

mcmc = data.rstan
coefs <- as.matrix(mcmc)[, 1:4]
newdata <- data.frame(A = levels(data$A), B = mean(data$B))
# A Tukeys contrast matrix
library(multcomp)
tuk.mat <- contrMat(n = table(newdata$A), type = "Tukey")
Xmat <- model.matrix(~A + B, data = newdata)
pairwise.mat <- tuk.mat %*% Xmat
pairwise.mat

                  (Intercept) AGroup B AGroup C B
Group B - Group A           0        1        0 0
Group C - Group A           0        0        1 0
Group C - Group B           0       -1        1 0

mcmc_areas(coefs %*% t(pairwise.mat))

(comps = tidyMCMC(coefs %*% t(pairwise.mat), conf.int = TRUE, conf.method = "HPDinterval"))

               term   estimate std.error   conf.low  conf.high
1 Group B - Group A -17.120600  1.916002 -20.721678 -13.023381
2 Group C - Group A -22.879458  1.927939 -26.818438 -19.313017
3 Group C - Group B  -5.758857  1.908307  -9.465031  -1.972144

ggplot(comps, aes(y = estimate, x = term)) + geom_pointrange(aes(ymin = conf.low,
    ymax = conf.high)) + geom_hline(yintercept = 0, linetype = "dashed") +
    scale_y_continuous("Effect size") + scale_x_discrete("") + coord_flip() +
    theme_classic()

With a couple of modifications, we could also express this as percentage changes. A percentage change represents the change (difference between groups) divided by one of the groups (determined by which group you want to express the percentage change to). Hence, we generate an additional mcmc matrix that represents the cell means for the divisor group (group we want to express change relative to).

Since the tuk.mat defines comparisons as -1 and 1 pairs, if we simply replace all the -1 with 0, the eventual matrix multiplication will result in estimates of the divisor cell means instread of the difference. We can then divide the original mcmc matrix above with this new divisor mcmc matrix to yeild a mcmc matrix of percentage change.

# Modify the tuk.mat to replace -1 with 0.  This will allow us to get a
# mcmc matrix of ..
tuk.mat[tuk.mat == -1] = 0
comp.mat <- tuk.mat %*% Xmat
comp.mat

                  (Intercept) AGroup B AGroup C        B
Group B - Group A           1        1        0 16.50175
Group C - Group A           1        0        1 16.50175
Group C - Group B           1        0        1 16.50175

comp.mcmc = 100 * (coefs %*% t(pairwise.mat))/coefs %*% t(comp.mat)
(comps = tidyMCMC(comp.mcmc, conf.int = TRUE, conf.method = "HPDinterval"))

               term   estimate std.error   conf.low  conf.high
1 Group B - Group A  -65.68005  9.938618  -86.57138 -46.874022
2 Group C - Group A -112.65860 15.767824 -145.18904 -84.337691
3 Group C - Group B  -28.69047 10.825916  -48.80312  -6.685448

ggplot(comps, aes(y = estimate, x = term)) + geom_pointrange(aes(ymin = conf.low,
    ymax = conf.high)) + geom_hline(yintercept = 0, linetype = "dashed") +
    scale_y_continuous("Effect size (%)") + scale_x_discrete("") + coord_flip() +
    theme_classic()

And now for the specific planned comparisons (Group B vs Group C as well as Group A vs the average of Groups B and C). This is achieved by generating our own contrast matrix (defining the contributions of each group to each contrast).

c.mat = rbind(c(0, 1, -1), c(1/2, -1/3, -1/3))
c.mat

     [,1]       [,2]       [,3]
[1,]  0.0  1.0000000 -1.0000000
[2,]  0.5 -0.3333333 -0.3333333

mcmc = data.rstan
coefs <- as.matrix(mcmc)[, 1:4]
newdata <- data.frame(A = levels(data$A), B = mean(data$B))
Xmat <- model.matrix(~A + B, data = newdata)
c.mat = c.mat %*% Xmat
c.mat

     (Intercept)   AGroup B   AGroup C         B
[1,]   0.0000000  1.0000000 -1.0000000  0.000000
[2,]  -0.1666667 -0.3333333 -0.3333333 -2.750292

(comps = tidyMCMC(as.mcmc(coefs %*% t(c.mat)), conf.int = TRUE, conf.method = "HPDinterval"))

  term estimate std.error conf.low conf.high
1 var1 5.758857 1.9083068 1.972144  9.465031
2 var2 6.107318 0.9312412 4.373877  8.159077

Lets start by comparing each group to each other group in a pairwise manner. Arguably the most elegant way to do this is to generate a Tukey's contrast matrix. This is a model matrix specific to comparing each group to each other group. Again, since the lines are parallel, it does not really matter what level of B we estimate these efffects at - so lets use the mean B.

mcmc = data.rstanarm
coefs <- as.matrix(mcmc)[, 1:4]
newdata <- data.frame(A = levels(data$A), B = mean(data$B))
# A Tukeys contrast matrix
library(multcomp)
tuk.mat <- contrMat(n = table(newdata$A), type = "Tukey")
Xmat <- model.matrix(~A + B, data = newdata)
pairwise.mat <- tuk.mat %*% Xmat
pairwise.mat

                  (Intercept) AGroup B AGroup C B
Group B - Group A           0        1        0 0
Group C - Group A           0        0        1 0
Group C - Group B           0       -1        1 0

mcmc_areas(coefs %*% t(pairwise.mat))

(comps = tidyMCMC(coefs %*% t(pairwise.mat), conf.int = TRUE, conf.method = "HPDinterval"))

               term   estimate std.error   conf.low  conf.high
1 Group B - Group A -17.211444  1.938649 -21.006648 -13.433809
2 Group C - Group A -22.903757  1.967936 -26.675721 -19.018964
3 Group C - Group B  -5.692313  1.892174  -9.159957  -1.799432

ggplot(comps, aes(y = estimate, x = term)) + geom_pointrange(aes(ymin = conf.low,
    ymax = conf.high)) + geom_hline(yintercept = 0, linetype = "dashed") +
    scale_y_continuous("Effect size") + scale_x_discrete("") + coord_flip() +
    theme_classic()

With a couple of modifications, we could also express this as percentage changes. A percentage change represents the change (difference between groups) divided by one of the groups (determined by which group you want to express the percentage change to). Hence, we generate an additional mcmc matrix that represents the cell means for the divisor group (group we want to express change relative to).

Since the tuk.mat defines comparisons as -1 and 1 pairs, if we simply replace all the -1 with 0, the eventual matrix multiplication will result in estimates of the divisor cell means instread of the difference. We can then divide the original mcmc matrix above with this new divisor mcmc matrix to yeild a mcmc matrix of percentage change.

# Modify the tuk.mat to replace -1 with 0.  This will allow us to get a
# mcmc matrix of ..
tuk.mat[tuk.mat == -1] = 0
comp.mat <- tuk.mat %*% Xmat
comp.mat

                  (Intercept) AGroup B AGroup C        B
Group B - Group A           1        1        0 16.50175
Group C - Group A           1        0        1 16.50175
Group C - Group B           1        0        1 16.50175

comp.mcmc = 100 * (coefs %*% t(pairwise.mat))/coefs %*% t(comp.mat)
(comps = tidyMCMC(comp.mcmc, conf.int = TRUE, conf.method = "HPDinterval"))

               term   estimate std.error   conf.low  conf.high
1 Group B - Group A  -66.16030  10.09871  -86.76507 -46.675359
2 Group C - Group A -112.70201  15.95929 -144.71321 -82.068118
3 Group C - Group B  -28.33604  10.75415  -50.15437  -8.304887

ggplot(comps, aes(y = estimate, x = term)) + geom_pointrange(aes(ymin = conf.low,
    ymax = conf.high)) + geom_hline(yintercept = 0, linetype = "dashed") +
    scale_y_continuous("Effect size (%)") + scale_x_discrete("") + coord_flip() +
    theme_classic()

And now for the specific planned comparisons (Group B vs Group C as well as Group A vs the average of Groups B and C). This is achieved by generating our own contrast matrix (defining the contributions of each group to each contrast).

c.mat = rbind(c(0, 1, -1), c(1/2, -1/3, -1/3))
c.mat

     [,1]       [,2]       [,3]
[1,]  0.0  1.0000000 -1.0000000
[2,]  0.5 -0.3333333 -0.3333333

mcmc = data.rstanarm
coefs <- as.matrix(mcmc)[, 1:4]
newdata <- data.frame(A = levels(data$A), B = mean(data$B))
Xmat <- model.matrix(~A + B, data = newdata)
c.mat = c.mat %*% Xmat
c.mat

     (Intercept)   AGroup B   AGroup C         B
[1,]   0.0000000  1.0000000 -1.0000000  0.000000
[2,]  -0.1666667 -0.3333333 -0.3333333 -2.750292

(comps = tidyMCMC(as.mcmc(coefs %*% t(c.mat)), conf.int = TRUE, conf.method = "HPDinterval"))

  term estimate std.error conf.low conf.high
1 var1 5.692313 1.8921744 1.799432  9.159957
2 var2 6.138562 0.9558204 4.235113  7.985014

Lets start by comparing each group to each other group in a pairwise manner. Arguably the most elegant way to do this is to generate a Tukey's contrast matrix. This is a model matrix specific to comparing each group to each other group. Again, since the lines are parallel, it does not really matter what level of B we estimate these efffects at - so lets use the mean B.

mcmc = data.brms
coefs <- as.matrix(mcmc)[, 1:4]
newdata <- data.frame(A = levels(data$A), B = mean(data$B))
# A Tukeys contrast matrix
library(multcomp)
tuk.mat <- contrMat(n = table(newdata$A), type = "Tukey")
Xmat <- model.matrix(~A + B, data = newdata)
pairwise.mat <- tuk.mat %*% Xmat
pairwise.mat

                  (Intercept) AGroup B AGroup C B
Group B - Group A           0        1        0 0
Group C - Group A           0        0        1 0
Group C - Group B           0       -1        1 0

mcmc_areas(coefs %*% t(pairwise.mat))

(comps = tidyMCMC(coefs %*% t(pairwise.mat), conf.int = TRUE, conf.method = "HPDinterval"))

               term   estimate std.error   conf.low  conf.high
1 Group B - Group A -17.172833  1.831351 -20.680311 -13.558229
2 Group C - Group A -22.891586  1.854493 -26.384632 -19.217887
3 Group C - Group B  -5.718753  1.846079  -9.169838  -2.065845

ggplot(comps, aes(y = estimate, x = term)) + geom_pointrange(aes(ymin = conf.low,
    ymax = conf.high)) + geom_hline(yintercept = 0, linetype = "dashed") +
    scale_y_continuous("Effect size") + scale_x_discrete("") + coord_flip() +
    theme_classic()

With a couple of modifications, we could also express this as percentage changes. A percentage change represents the change (difference between groups) divided by one of the groups (determined by which group you want to express the percentage change to). Hence, we generate an additional mcmc matrix that represents the cell means for the divisor group (group we want to express change relative to).

Since the tuk.mat defines comparisons as -1 and 1 pairs, if we simply replace all the -1 with 0, the eventual matrix multiplication will result in estimates of the divisor cell means instread of the difference. We can then divide the original mcmc matrix above with this new divisor mcmc matrix to yeild a mcmc matrix of percentage change.

# Modify the tuk.mat to replace -1 with 0.  This will allow us to get a
# mcmc matrix of ..
tuk.mat[tuk.mat == -1] = 0
comp.mat <- tuk.mat %*% Xmat
comp.mat

                  (Intercept) AGroup B AGroup C        B
Group B - Group A           1        1        0 16.50175
Group C - Group A           1        0        1 16.50175
Group C - Group B           1        0        1 16.50175

comp.mcmc = 100 * (coefs %*% t(pairwise.mat))/coefs %*% t(comp.mat)
(comps = tidyMCMC(comp.mcmc, conf.int = TRUE, conf.method = "HPDinterval"))

               term   estimate std.error   conf.low  conf.high
1 Group B - Group A  -65.90339  9.650503  -84.05345 -46.582017
2 Group C - Group A -112.55948 15.341280 -143.66383 -84.530670
3 Group C - Group B  -28.43318 10.498519  -47.72283  -7.691681

ggplot(comps, aes(y = estimate, x = term)) + geom_pointrange(aes(ymin = conf.low,
    ymax = conf.high)) + geom_hline(yintercept = 0, linetype = "dashed") +
    scale_y_continuous("Effect size (%)") + scale_x_discrete("") + coord_flip() +
    theme_classic()

And now for the specific planned comparisons (Group B vs Group C as well as Group A vs the average of Groups B and C). This is achieved by generating our own contrast matrix (defining the contributions of each group to each contrast).

c.mat = rbind(c(0, 1, -1), c(1/2, -1/3, -1/3))
c.mat

     [,1]       [,2]       [,3]
[1,]  0.0  1.0000000 -1.0000000
[2,]  0.5 -0.3333333 -0.3333333

mcmc = data.brms
coefs <- as.matrix(mcmc)[, 1:4]
newdata <- data.frame(A = levels(data$A), B = mean(data$B))
Xmat <- model.matrix(~A + B, data = newdata)
c.mat = c.mat %*% Xmat
c.mat

     (Intercept)   AGroup B   AGroup C         B
[1,]   0.0000000  1.0000000 -1.0000000  0.000000
[2,]  -0.1666667 -0.3333333 -0.3333333 -2.750292

(comps = tidyMCMC(as.mcmc(coefs %*% t(c.mat)), conf.int = TRUE, conf.method = "HPDinterval"))

  term estimate std.error conf.low conf.high
1 var1 5.718753 1.8460785 2.065845  9.169838
2 var2 6.122501 0.8984429 4.344599  7.808558

Variance components, the amount of added variance attributed to each influence, are traditionally estimated for so called random effects. These are the effects for which the levels employed in the design are randomly selected to represent a broader range of possible levels. For such effects, effect sizes (differences between each level and a reference level) are of little value. Instead, the 'importance' of the variables are measured in units of variance components.

On the other hand, regular variance components for fixed factors (those whose measured levels represent the only levels of interest) are not logical - since variance components estimate variance as if the levels are randomly selected from a larger population. Nevertheless, in order to compare and contrast the scale of variability of both fixed and random factors, it is necessary to measure both on the same scale (sample or population based variance).

Finite-population variance components (Gelman, 2005) assume that the levels of all factors (fixed and random) in the design are all the possible levels available. In other words, they are assumed to represent finite populations of levels. Sample (rather than population) statistics are then used to calculate these finite-population variances (or standard deviations).

Since standard deviation (and variance) are bound at zero, standard deviation posteriors are typically non-normal. Consequently, medians and HPD intervals are more robust estimates.

library(broom)
mcmc = data.mcmcpack
head(mcmc)

Markov Chain Monte Carlo (MCMC) output:
Start = 1001 
End = 1007 
Thinning interval = 1 
     (Intercept)  AGroup B  AGroup C          B   sigma2
[1,]    50.18839 -18.23379 -22.87414 -0.4394640 14.11917
[2,]    49.39556 -15.64877 -20.98086 -0.4528605 13.22004
[3,]    49.44854 -16.36367 -23.66886 -0.3887300 15.57976
[4,]    51.46891 -16.26642 -24.26313 -0.4709924 12.01597
[5,]    50.55617 -17.54613 -22.44262 -0.4252965 16.35805
[6,]    50.73340 -17.54078 -23.75846 -0.4598987 10.84731
[7,]    51.16774 -17.41150 -22.41515 -0.3996571 17.57170

# A
wch = grep("^A", colnames(mcmc))
# Get the rowwise standard deviations between effects parameters
sd.A = apply(mcmc[, wch], 1, sd)
# B
wch = grep("^B", colnames(mcmc))
sd.B = sd(data$B) * abs(mcmc[, wch])
# generate a model matrix
newdata = data
Xmat = model.matrix(~A + B, newdata)
## get median parameter estimates
wch = grep("(Intercept)|^A|^B", colnames(mcmc))
coefs = mcmc[, wch]
fit = coefs %*% t(Xmat)
resid = sweep(fit, 2, data$Y, "-")
sd.resid = apply(resid, 1, sd)

sd.all = cbind(sd.A, sd.B, sd.resid)
(fpsd = tidyMCMC(sd.all, conf.int = TRUE, conf.method = "HPDinterval"))

      term estimate std.error conf.low conf.high
1     sd.A 4.080665 1.2858143 1.606486  6.754944
2     sd.B 5.030454 0.7743607 3.478134  6.527505
3 sd.resid 3.983878 0.1906139 3.761980  4.356994

# OR expressed as a percentage
(fpsd.p = tidyMCMC(100 * sd.all/rowSums(sd.all), estimate.method = "median",
    conf.int = TRUE, conf.method = "HPDinterval"))

      term estimate std.error conf.low conf.high
1     sd.A 31.28801  7.123578 16.55312  44.15266
2     sd.B 38.61193  5.120548 28.29820  48.32862
3 sd.resid 30.05025  4.065060 24.60465  38.59845

## we can even plot this as a Bayesian ANOVA table
ggplot(fpsd, aes(y = estimate, x = term)) + geom_hline(yintercept = 0,
    linetype = "dashed") + geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
    geom_text(aes(label = sprintf("%.2f%%", fpsd.p$estimate), vjust = -1)) +
    scale_y_continuous("Finite population standard deviation") + scale_x_discrete() +
    coord_flip() + theme_classic()

Conclusions: Approximately 31.3% of the total finite population standard deviation is due to x.

library(broom)
mcmc = data.r2jags$BUGSoutput$sims.matrix
head(mcmc)

      beta[1]   beta[2]   beta[3]    beta[4] deviance    sigma
[1,] 52.60917 -18.82408 -23.38489 -0.4880702 167.0879 4.283265
[2,] 53.75233 -18.89312 -24.92630 -0.5191816 171.4824 5.124059
[3,] 51.15722 -16.81356 -23.25841 -0.5336914 170.4625 3.705732
[4,] 50.32853 -16.98682 -23.41316 -0.4149829 165.9603 4.545555
[5,] 48.78829 -15.18438 -24.19762 -0.4395482 172.1889 4.458376
[6,] 51.55891 -17.64261 -22.37656 -0.5040106 166.7893 4.155680

# A
wch = grep("beta.[2-3]", colnames(mcmc))
# Get the rowwise standard deviations between effects parameters
sd.A = apply(mcmc[, wch], 1, sd)
# B
wch = grep("beta.4", colnames(mcmc))
sd.B = sd(data$B) * abs(mcmc[, wch])
# generate a model matrix
newdata = data
Xmat = model.matrix(~A + B, newdata)
## get median parameter estimates
wch = grep("beta.[1-4]", colnames(mcmc))
coefs = mcmc[, wch]
fit = coefs %*% t(Xmat)
resid = sweep(fit, 2, data$Y, "-")
sd.resid = apply(resid, 1, sd)

sd.all = cbind(sd.A, sd.B, sd.resid)
(fpsd = tidyMCMC(sd.all, conf.int = TRUE, conf.method = "HPDinterval"))

      term estimate std.error conf.low conf.high
1     sd.A 4.085929 1.3242287 1.399537  6.648241
2     sd.B 5.031876 0.7938012 3.468931  6.606399
3 sd.resid 3.991447 0.1970424 3.762837  4.379866

# OR expressed as a percentage
(fpsd.p = tidyMCMC(100 * sd.all/rowSums(sd.all), estimate.method = "median",
    conf.int = TRUE, conf.method = "HPDinterval"))

      term estimate std.error conf.low conf.high
1     sd.A 31.33381  7.376887 15.57164  43.97278
2     sd.B 38.65910  5.313671 27.75097  48.62132
3 sd.resid 29.97544  4.188131 24.58803  39.19370

## we can even plot this as a Bayesian ANOVA table
ggplot(fpsd, aes(y = estimate, x = term)) + geom_hline(yintercept = 0,
    linetype = "dashed") + geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
    geom_text(aes(label = sprintf("%.2f%%", fpsd.p$estimate), vjust = -1)) +
    scale_y_continuous("Finite population standard deviation") + scale_x_discrete() +
    coord_flip() + theme_classic()

Conclusions: Approximately 31.3% of the total finite population standard deviation is due to x.