Tutorial 7.3b - Multiple linear regression (Bayesian)

Multiple linear regression models can include predictors (terms) that are incorporated additively (no interactions) or multiplicatively (with interactions). As such we will explore these separately for each modelling tool.

The observed response ($y_i$) are assumed to be drawn from a normal distribution with a given mean ($\mu$) and standard deviation ($\sigma$). The expected values ($\mu$) are themselves determined by the linear predictor ($\beta_0 + \beta X_i$). In this case, $\beta_0$ represents the y-intercept (value of $y$ when all of the $x$'s are equal to zero) and the set of $\beta$'s represent the rates of change in $y$ for every unit change in each $x$ (the effect) holding each other $x$ constant.

Note that since we should always center all predictors (by subtracting the mean of each $x$ from the repective values of each $x$), the y-intercept represents the value of $y$ at the average value of each $x$.

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

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

Additive model

library(MCMCpack)
data.mcmcpack.add <- MCMCregress(y ~ cx1 + cx2, data = data)

Multiplicative model

library(MCMCpack)
data.mcmcpack.mult <- MCMCregress(y ~ cx1 * cx2, data = data)

If we define the model in terms of matrices, the JAGS model definition is identical for both additive and multiplicative models. $$\begin{align} y_i&\sim{}N(\mu_i, \tau)\\ \mu_i &= \beta_0 + \beta X_i\\ \beta_0&\sim{}N(0,1.0{E-6}) \hspace{1cm}\mathsf{non-informative~prior~for~interept}\\ \beta_j&\sim{}N(0,1.0{E-6}) \hspace{1cm}\mathsf{non-informative~prior~for~partial~slopes}\\ \tau &= 1/\sigma^2\\ \sigma&\sim{}U(0,100)\\ \end{align} $$

Define the model

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

Additive model

Define the data list

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 matrix (X)
• the number of predictor variables (nX)
• the total number of observed items (n)

X = model.matrix(~cx1 + cx2, data = data)
data.list <- with(data, list(y = y, X = X[, -1], nX = ncol(X) -
    1, n = nrow(data)))

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("beta0", "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.

## load the R2jags package
library(R2jags)

data.r2jags.add <- 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: 100
   Unobserved stochastic nodes: 4
   Total graph size: 618

Initializing model

print(data.r2jags.add)

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]    3.028   0.501   2.037   2.692   3.032   3.365   4.009 1.001 15000
beta[2]    1.389   0.432   0.545   1.101   1.394   1.675   2.235 1.001 15000
beta0      3.823   0.115   3.597   3.745   3.822   3.901   4.048 1.001  8000
sigma      1.146   0.083   1.001   1.088   1.140   1.198   1.322 1.001 15000
deviance 309.526   2.883 305.925 307.427 308.888 310.918 316.903 1.002  3100

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 = 4.2 and DIC = 313.7
DIC is an estimate of expected predictive error (lower deviance is better).

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

Multiplicative model

Define the data list

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 matrix (X)
• the number of predictor terms (nX)
• the total number of observed items (n)

X = model.matrix(~cx1 * cx2, data = data)
data.list <- with(data, list(y = y, X = X[, -1], nX = ncol(X) - 1, n = nrow(data)))

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("beta0", "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.

## load the R2jags package
library(R2jags)

data.r2jags.mult <- 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: 100
   Unobserved stochastic nodes: 5
   Total graph size: 721

Initializing model

print(data.r2jags.mult)

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]    2.931   0.499   1.969   2.593   2.929   3.265   3.902 1.001 15000
beta[2]    1.344   0.426   0.508   1.053   1.343   1.631   2.181 1.001 15000
beta[3]    2.675   1.256   0.182   1.845   2.669   3.509   5.158 1.001 15000
beta0      3.671   0.134   3.413   3.580   3.671   3.760   3.932 1.001 12000
sigma      1.126   0.082   0.981   1.070   1.120   1.176   1.307 1.001  8000
deviance 305.828   3.280 301.532 303.441 305.137 307.457 314.045 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 = 5.4 and DIC = 311.2
DIC is an estimate of expected predictive error (lower deviance is better).

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

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.

Define the model

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

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

					  data { 
					  int<lower=1> n;   // total number of observations 
					  vector[n] Y;      // response variable 
					  int<lower=1> nX;  // number of effects 
					  matrix[n, nX] X;   // model matrix 
					  } 
					  transformed data { 
					  matrix[n, nX - 1] Xc;  // centered version of X 
					  vector[nX - 1] means_X;  // column means of X before centering 
					  
					  for (i in 2:nX) { 
					  means_X[i - 1] = mean(X[, i]); 
					  Xc[, i - 1] = X[, i] - means_X[i - 1]; 
					  }  
					  } 
					  parameters { 
					  vector[nX-1] beta;  // population-level effects 
					  real cbeta0;  // center-scale intercept 
					  real<lower=0> sigma;  // residual SD 
					  } 
					  transformed parameters { 
					  } 
					  model { 
					  vector[n] mu; 
					  mu = Xc * beta + cbeta0; 
					  // prior specifications 
					  beta ~ normal(0, 100); 
					  cbeta0 ~ normal(0, 100); 
					  sigma ~ cauchy(0, 5); 
					  // likelihood contribution 
					  Y ~ normal(mu, sigma); 
					  } 
					  generated quantities {
					  real beta0;  // population-level intercept 
					  vector[n] log_lik;
					  beta0 = cbeta0 - dot_product(means_X, beta);
					  for (i in 1:n) {
					  log_lik[i] = normal_lpdf(Y[i] | Xc[i] * beta + cbeta0, sigma);
					  } 
					  }

Additive model

Define the data list

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 matrix (X)
• the number of predictor variables (nX)
• the total number of observed items (n)

X = model.matrix(~cx1 + cx2, data = data)
data.list <- with(data, list(Y = y, X = X, nX = ncol(X), n = nrow(data)))

Define the MCMC chain parameters

Next we should define the behavioural parameters of the No-U-Turn 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

nChains = 3
burnInSteps = 1000
thinSteps = 3
numSavedSteps = 3000  #across all chains
nIter = ceiling(burnInSteps + (numSavedSteps * thinSteps)/nChains)
nIter

[1] 4000

Fit the model

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

During the warmup stage, the No-U-Turn sampler (NUTS) attempts to determine the optimum stepsize - the stepsize that achieves the target acceptance rate (0.8 or 80% by default) without divergence (occurs when the stepsize is too large relative to the curvature of the log posterior - and results in approximations that are likely to diverge and be biased) and without hitting the maximum treedepth (10). At each iteration of the NUTS algorithm, the number of leapfrog steps doubles (as it increases the treedepth) and only terminates when either the NUTS criterion are satisfied or the tree depth reaches the maximum (10 by default).

## load the rstan package
library(rstan)

data.rstan.add <- stan(data = data.list, model_code = modelString, chains = nChains,
    iter = nIter, warmup = burnInSteps, thin = thinSteps, save_dso = TRUE)

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 file2b18778a1a1f.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 '4e0a54cd22440d9847d7705bdb1ff803' NOW (CHAIN 1).

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


Iteration:    1 / 4000 [  0%]  (Warmup)
Iteration:  400 / 4000 [ 10%]  (Warmup)
Iteration:  800 / 4000 [ 20%]  (Warmup)
Iteration: 1001 / 4000 [ 25%]  (Sampling)
Iteration: 1400 / 4000 [ 35%]  (Sampling)
Iteration: 1800 / 4000 [ 45%]  (Sampling)
Iteration: 2200 / 4000 [ 55%]  (Sampling)
Iteration: 2600 / 4000 [ 65%]  (Sampling)
Iteration: 3000 / 4000 [ 75%]  (Sampling)
Iteration: 3400 / 4000 [ 85%]  (Sampling)
Iteration: 3800 / 4000 [ 95%]  (Sampling)
Iteration: 4000 / 4000 [100%]  (Sampling)

 Elapsed Time: 0.073416 seconds (Warm-up)
               0.143677 seconds (Sampling)
               0.217093 seconds (Total)


SAMPLING FOR MODEL '4e0a54cd22440d9847d7705bdb1ff803' NOW (CHAIN 2).

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


Iteration:    1 / 4000 [  0%]  (Warmup)
Iteration:  400 / 4000 [ 10%]  (Warmup)
Iteration:  800 / 4000 [ 20%]  (Warmup)
Iteration: 1001 / 4000 [ 25%]  (Sampling)
Iteration: 1400 / 4000 [ 35%]  (Sampling)
Iteration: 1800 / 4000 [ 45%]  (Sampling)
Iteration: 2200 / 4000 [ 55%]  (Sampling)
Iteration: 2600 / 4000 [ 65%]  (Sampling)
Iteration: 3000 / 4000 [ 75%]  (Sampling)
Iteration: 3400 / 4000 [ 85%]  (Sampling)
Iteration: 3800 / 4000 [ 95%]  (Sampling)
Iteration: 4000 / 4000 [100%]  (Sampling)

 Elapsed Time: 0.060339 seconds (Warm-up)
               0.177215 seconds (Sampling)
               0.237554 seconds (Total)


SAMPLING FOR MODEL '4e0a54cd22440d9847d7705bdb1ff803' NOW (CHAIN 3).

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 / 4000 [  0%]  (Warmup)
Iteration:  400 / 4000 [ 10%]  (Warmup)
Iteration:  800 / 4000 [ 20%]  (Warmup)
Iteration: 1001 / 4000 [ 25%]  (Sampling)
Iteration: 1400 / 4000 [ 35%]  (Sampling)
Iteration: 1800 / 4000 [ 45%]  (Sampling)
Iteration: 2200 / 4000 [ 55%]  (Sampling)
Iteration: 2600 / 4000 [ 65%]  (Sampling)
Iteration: 3000 / 4000 [ 75%]  (Sampling)
Iteration: 3400 / 4000 [ 85%]  (Sampling)
Iteration: 3800 / 4000 [ 95%]  (Sampling)
Iteration: 4000 / 4000 [100%]  (Sampling)

 Elapsed Time: 0.0688 seconds (Warm-up)
               0.149029 seconds (Sampling)
               0.217829 seconds (Total)

data.rstan.add

Inference for Stan model: 4e0a54cd22440d9847d7705bdb1ff803.
3 chains, each with iter=4000; warmup=1000; thin=3; 
post-warmup draws per chain=1000, total post-warmup draws=3000.

               mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
beta[1]        3.02    0.01 0.51   2.03   2.69   3.02   3.34   4.05  2747    1
beta[2]        1.38    0.01 0.43   0.52   1.10   1.39   1.68   2.22  2854    1
cbeta0         3.83    0.00 0.12   3.59   3.75   3.83   3.91   4.06  2317    1
sigma          1.15    0.00 0.08   0.99   1.09   1.14   1.20   1.32  2806    1
beta0          3.83    0.00 0.12   3.59   3.75   3.83   3.91   4.06  2317    1
log_lik[1]    -1.21    0.00 0.11  -1.46  -1.28  -1.20  -1.13  -1.02  2200    1
log_lik[2]    -1.11    0.00 0.09  -1.31  -1.17  -1.11  -1.05  -0.96  2916    1
log_lik[3]    -1.45    0.00 0.14  -1.74  -1.54  -1.44  -1.35  -1.21  2829    1
log_lik[4]    -1.16    0.00 0.09  -1.35  -1.22  -1.16  -1.10  -1.00  2219    1
log_lik[5]    -1.37    0.00 0.10  -1.59  -1.44  -1.37  -1.30  -1.19  2327    1
log_lik[6]    -1.12    0.00 0.08  -1.28  -1.17  -1.12  -1.07  -0.98  2313    1
log_lik[7]    -1.27    0.00 0.16  -1.63  -1.36  -1.24  -1.15  -1.02  2806    1
log_lik[8]    -2.43    0.01 0.25  -2.95  -2.59  -2.42  -2.26  -1.99  2248    1
log_lik[9]    -1.42    0.00 0.16  -1.77  -1.51  -1.40  -1.31  -1.15  2796    1
log_lik[10]   -1.09    0.00 0.08  -1.24  -1.14  -1.09  -1.04  -0.95  2402    1
log_lik[11]   -1.21    0.00 0.08  -1.39  -1.27  -1.21  -1.15  -1.05  2839    1
log_lik[12]   -1.16    0.00 0.10  -1.37  -1.22  -1.15  -1.09  -0.98  2774    1
log_lik[13]   -1.22    0.00 0.10  -1.44  -1.28  -1.21  -1.14  -1.03  2784    1
log_lik[14]   -1.81    0.00 0.20  -2.24  -1.94  -1.80  -1.67  -1.46  2828    1
log_lik[15]   -1.14    0.00 0.11  -1.40  -1.20  -1.13  -1.07  -0.96  2934    1
log_lik[16]   -1.12    0.00 0.09  -1.32  -1.17  -1.11  -1.06  -0.96  2917    1
log_lik[17]   -1.07    0.00 0.08  -1.24  -1.12  -1.07  -1.02  -0.93  2768    1
log_lik[18]   -1.07    0.00 0.08  -1.23  -1.12  -1.07  -1.02  -0.93  2581    1
log_lik[19]   -1.07    0.00 0.08  -1.24  -1.12  -1.07  -1.02  -0.92  2737    1
log_lik[20]   -1.06    0.00 0.07  -1.21  -1.11  -1.06  -1.01  -0.92  2809    1
log_lik[21]   -1.98    0.00 0.22  -2.45  -2.12  -1.96  -1.82  -1.58  2405    1
log_lik[22]   -1.16    0.00 0.12  -1.44  -1.22  -1.15  -1.08  -0.97  2216    1
log_lik[23]   -1.17    0.00 0.11  -1.42  -1.23  -1.16  -1.09  -0.98  2231    1
log_lik[24]   -1.07    0.00 0.08  -1.25  -1.13  -1.07  -1.02  -0.92  2452    1
log_lik[25]   -1.24    0.00 0.10  -1.45  -1.30  -1.23  -1.17  -1.07  2104    1
log_lik[26]   -1.37    0.00 0.12  -1.63  -1.45  -1.36  -1.29  -1.16  2896    1
log_lik[27]   -1.09    0.00 0.08  -1.24  -1.14  -1.09  -1.04  -0.94  2487    1
log_lik[28]   -1.27    0.00 0.13  -1.55  -1.35  -1.26  -1.18  -1.06  2881    1
log_lik[29]   -2.30    0.01 0.31  -2.98  -2.50  -2.28  -2.08  -1.75  2682    1
log_lik[30]   -1.40    0.00 0.12  -1.67  -1.48  -1.39  -1.32  -1.18  2285    1
log_lik[31]   -1.22    0.00 0.09  -1.40  -1.27  -1.21  -1.16  -1.06  2306    1
log_lik[32]   -1.06    0.00 0.07  -1.21  -1.11  -1.06  -1.01  -0.92  2518    1
log_lik[33]   -2.19    0.00 0.27  -2.77  -2.36  -2.18  -2.01  -1.73  2877    1
log_lik[34]   -2.26    0.00 0.24  -2.78  -2.42  -2.25  -2.09  -1.83  2785    1
log_lik[35]   -1.07    0.00 0.08  -1.22  -1.12  -1.07  -1.02  -0.93  2802    1
log_lik[36]   -3.41    0.01 0.37  -4.20  -3.64  -3.39  -3.15  -2.74  2146    1
log_lik[37]   -1.53    0.00 0.21  -2.00  -1.65  -1.51  -1.38  -1.18  2829    1
log_lik[38]   -1.22    0.00 0.13  -1.50  -1.29  -1.20  -1.13  -1.01  2842    1
log_lik[39]   -1.86    0.00 0.15  -2.17  -1.95  -1.85  -1.75  -1.59  2213    1
log_lik[40]   -1.32    0.00 0.11  -1.55  -1.39  -1.31  -1.24  -1.12  2662    1
log_lik[41]   -1.58    0.00 0.13  -1.85  -1.67  -1.58  -1.49  -1.35  2326    1
log_lik[42]   -2.93    0.01 0.36  -3.71  -3.15  -2.91  -2.68  -2.30  2290    1
log_lik[43]   -1.43    0.00 0.14  -1.73  -1.51  -1.42  -1.32  -1.18  2520    1
log_lik[44]   -1.09    0.00 0.08  -1.26  -1.14  -1.08  -1.03  -0.94  2925    1
log_lik[45]   -1.31    0.00 0.13  -1.60  -1.39  -1.30  -1.22  -1.09  2624    1
log_lik[46]   -1.16    0.00 0.09  -1.36  -1.22  -1.16  -1.10  -1.00  2348    1
log_lik[47]   -1.47    0.00 0.17  -1.84  -1.57  -1.45  -1.34  -1.18  2530    1
log_lik[48]   -1.23    0.00 0.13  -1.52  -1.31  -1.22  -1.14  -1.03  2654    1
log_lik[49]   -2.42    0.00 0.24  -2.94  -2.57  -2.40  -2.24  -2.00  2770    1
log_lik[50]   -1.07    0.00 0.08  -1.23  -1.12  -1.07  -1.02  -0.93  2752    1
log_lik[51]   -1.78    0.00 0.21  -2.21  -1.91  -1.76  -1.63  -1.42  2456    1
log_lik[52]   -4.28    0.01 0.65  -5.68  -4.69  -4.24  -3.82  -3.13  2861    1
log_lik[53]   -1.33    0.00 0.18  -1.75  -1.43  -1.31  -1.21  -1.05  2878    1
log_lik[54]   -1.30    0.00 0.15  -1.64  -1.38  -1.28  -1.19  -1.05  2271    1
log_lik[55]   -1.59    0.00 0.19  -2.01  -1.70  -1.58  -1.46  -1.28  2922    1
log_lik[56]   -1.07    0.00 0.08  -1.23  -1.12  -1.07  -1.02  -0.92  2743    1
log_lik[57]   -1.07    0.00 0.08  -1.22  -1.12  -1.07  -1.02  -0.92  2799    1
log_lik[58]   -1.29    0.00 0.15  -1.62  -1.38  -1.27  -1.19  -1.05  2840    1
log_lik[59]   -1.52    0.00 0.16  -1.86  -1.62  -1.51  -1.41  -1.25  2484    1
log_lik[60]   -3.72    0.01 0.59  -5.02  -4.09  -3.68  -3.30  -2.69  2836    1
log_lik[61]   -1.14    0.00 0.11  -1.40  -1.20  -1.13  -1.07  -0.96  2936    1
log_lik[62]   -1.63    0.00 0.20  -2.07  -1.75  -1.61  -1.49  -1.30  2649    1
log_lik[63]   -1.30    0.00 0.11  -1.55  -1.37  -1.29  -1.23  -1.10  2260    1
log_lik[64]   -2.03    0.01 0.31  -2.71  -2.23  -2.01  -1.81  -1.52  2685    1
log_lik[65]   -2.47    0.01 0.33  -3.19  -2.67  -2.44  -2.24  -1.90  2972    1
log_lik[66]   -1.21    0.00 0.10  -1.43  -1.27  -1.20  -1.14  -1.03  2514    1
log_lik[67]   -1.81    0.00 0.27  -2.43  -1.97  -1.78  -1.61  -1.36  2935    1
log_lik[68]   -1.75    0.00 0.18  -2.15  -1.87  -1.74  -1.62  -1.44  2830    1
log_lik[69]   -1.18    0.00 0.10  -1.40  -1.24  -1.18  -1.11  -1.01  2569    1
log_lik[70]   -1.85    0.00 0.22  -2.33  -1.99  -1.84  -1.69  -1.47  2684    1
log_lik[71]   -1.07    0.00 0.07  -1.22  -1.12  -1.07  -1.02  -0.93  2542    1
log_lik[72]   -1.30    0.00 0.14  -1.61  -1.39  -1.28  -1.20  -1.07  2787    1
log_lik[73]   -1.11    0.00 0.09  -1.30  -1.16  -1.11  -1.05  -0.95  2444    1
log_lik[74]   -1.09    0.00 0.09  -1.28  -1.14  -1.08  -1.03  -0.94  2834    1
log_lik[75]   -1.25    0.00 0.10  -1.47  -1.32  -1.24  -1.17  -1.06  2795    1
log_lik[76]   -1.46    0.00 0.10  -1.67  -1.52  -1.46  -1.39  -1.28  2265    1
log_lik[77]   -2.93    0.01 0.37  -3.73  -3.17  -2.91  -2.66  -2.29  2806    1
log_lik[78]   -2.23    0.01 0.30  -2.86  -2.42  -2.21  -2.02  -1.70  2431    1
log_lik[79]   -2.47    0.01 0.32  -3.19  -2.67  -2.44  -2.25  -1.91  2812    1
log_lik[80]   -1.54    0.00 0.14  -1.84  -1.62  -1.52  -1.44  -1.29  2938    1
log_lik[81]   -1.07    0.00 0.08  -1.23  -1.12  -1.07  -1.02  -0.92  2751    1
log_lik[82]   -1.54    0.00 0.13  -1.82  -1.62  -1.53  -1.45  -1.31  2380    1
log_lik[83]   -1.11    0.00 0.09  -1.30  -1.17  -1.11  -1.06  -0.96  2357    1
log_lik[84]   -1.45    0.00 0.22  -1.97  -1.58  -1.42  -1.29  -1.10  2814    1
log_lik[85]   -1.22    0.00 0.13  -1.51  -1.29  -1.20  -1.13  -1.02  2712    1
log_lik[86]   -1.27    0.00 0.12  -1.52  -1.34  -1.26  -1.19  -1.07  2698    1
log_lik[87]   -1.84    0.00 0.24  -2.36  -2.00  -1.82  -1.67  -1.42  2501    1
log_lik[88]   -1.17    0.00 0.10  -1.37  -1.23  -1.16  -1.10  -1.00  2797    1
log_lik[89]   -2.54    0.00 0.26  -3.07  -2.71  -2.52  -2.36  -2.08  2900    1
log_lik[90]   -1.19    0.00 0.11  -1.45  -1.26  -1.18  -1.11  -1.01  2699    1
log_lik[91]   -1.42    0.00 0.17  -1.80  -1.53  -1.41  -1.31  -1.14  2488    1
log_lik[92]   -1.31    0.00 0.13  -1.59  -1.39  -1.29  -1.21  -1.09  2640    1
log_lik[93]   -1.17    0.00 0.08  -1.34  -1.22  -1.17  -1.12  -1.02  2268    1
log_lik[94]   -1.29    0.00 0.12  -1.55  -1.36  -1.28  -1.21  -1.08  2884    1
log_lik[95]   -4.17    0.01 0.62  -5.49  -4.55  -4.13  -3.75  -3.06  2957    1
log_lik[96]   -1.72    0.00 0.16  -2.05  -1.82  -1.70  -1.60  -1.43  2421    1
log_lik[97]   -1.07    0.00 0.08  -1.22  -1.12  -1.07  -1.02  -0.93  2860    1
log_lik[98]   -1.11    0.00 0.09  -1.31  -1.17  -1.11  -1.05  -0.96  2847    1
log_lik[99]   -1.54    0.00 0.15  -1.86  -1.63  -1.53  -1.44  -1.28  2492    1
log_lik[100]  -1.38    0.00 0.12  -1.63  -1.45  -1.37  -1.30  -1.17  2828    1
lp__         -62.83    0.03 1.45 -66.53 -63.52 -62.53 -61.76 -61.01  2476    1

Samples were drawn using NUTS(diag_e) at Thu Aug 17 15:32:58 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).

Multiplicative model

Define the data list

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 matrix (X)
• the number of predictor variables (nX)
• the total number of observed items (n)

X = model.matrix(~cx1 * cx2, data = data)
data.list <- with(data, list(Y = y, X = X, nX = ncol(X), n = nrow(data)))

Define the MCMC chain parameters

Next we should define the behavioural parameters of the No-U-Turn 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

nChains = 3
burnInSteps = 1000
thinSteps = 3
numSavedSteps = 3000  #across all chains
nIter = ceiling(burnInSteps + (numSavedSteps * thinSteps)/nChains)
nIter

[1] 4000

Fit the model

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

During the warmup stage, the No-U-Turn sampler (NUTS) attempts to determine the optimum stepsize - the stepsize that achieves the target acceptance rate (0.8 or 80% by default) without divergence (occurs when the stepsize is too large relative to the curvature of the log posterior - and results in approximations that are likely to diverge and be biased) and without hitting the maximum treedepth (10). At each iteration of the NUTS algorithm, the number of leapfrog steps doubles (as it increases the treedepth) and only terminates when either the NUTS criterion are satisfied or the tree depth reaches the maximum (10 by default).

## load the rstan package
library(rstan)

data.rstan.mult <- stan(data = data.list, model_code = modelString, chains = nChains,
    iter = nIter, warmup = burnInSteps, thin = thinSteps, save_dso = TRUE)

SAMPLING FOR MODEL '4e0a54cd22440d9847d7705bdb1ff803' NOW (CHAIN 1).

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


Iteration:    1 / 4000 [  0%]  (Warmup)
Iteration:  400 / 4000 [ 10%]  (Warmup)
Iteration:  800 / 4000 [ 20%]  (Warmup)
Iteration: 1001 / 4000 [ 25%]  (Sampling)
Iteration: 1400 / 4000 [ 35%]  (Sampling)
Iteration: 1800 / 4000 [ 45%]  (Sampling)
Iteration: 2200 / 4000 [ 55%]  (Sampling)
Iteration: 2600 / 4000 [ 65%]  (Sampling)
Iteration: 3000 / 4000 [ 75%]  (Sampling)
Iteration: 3400 / 4000 [ 85%]  (Sampling)
Iteration: 3800 / 4000 [ 95%]  (Sampling)
Iteration: 4000 / 4000 [100%]  (Sampling)

 Elapsed Time: 0.075668 seconds (Warm-up)
               0.14075 seconds (Sampling)
               0.216418 seconds (Total)


SAMPLING FOR MODEL '4e0a54cd22440d9847d7705bdb1ff803' NOW (CHAIN 2).

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


Iteration:    1 / 4000 [  0%]  (Warmup)
Iteration:  400 / 4000 [ 10%]  (Warmup)
Iteration:  800 / 4000 [ 20%]  (Warmup)
Iteration: 1001 / 4000 [ 25%]  (Sampling)
Iteration: 1400 / 4000 [ 35%]  (Sampling)
Iteration: 1800 / 4000 [ 45%]  (Sampling)
Iteration: 2200 / 4000 [ 55%]  (Sampling)
Iteration: 2600 / 4000 [ 65%]  (Sampling)
Iteration: 3000 / 4000 [ 75%]  (Sampling)
Iteration: 3400 / 4000 [ 85%]  (Sampling)
Iteration: 3800 / 4000 [ 95%]  (Sampling)
Iteration: 4000 / 4000 [100%]  (Sampling)

 Elapsed Time: 0.058992 seconds (Warm-up)
               0.139282 seconds (Sampling)
               0.198274 seconds (Total)


SAMPLING FOR MODEL '4e0a54cd22440d9847d7705bdb1ff803' NOW (CHAIN 3).

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 / 4000 [  0%]  (Warmup)
Iteration:  400 / 4000 [ 10%]  (Warmup)
Iteration:  800 / 4000 [ 20%]  (Warmup)
Iteration: 1001 / 4000 [ 25%]  (Sampling)
Iteration: 1400 / 4000 [ 35%]  (Sampling)
Iteration: 1800 / 4000 [ 45%]  (Sampling)
Iteration: 2200 / 4000 [ 55%]  (Sampling)
Iteration: 2600 / 4000 [ 65%]  (Sampling)
Iteration: 3000 / 4000 [ 75%]  (Sampling)
Iteration: 3400 / 4000 [ 85%]  (Sampling)
Iteration: 3800 / 4000 [ 95%]  (Sampling)
Iteration: 4000 / 4000 [100%]  (Sampling)

 Elapsed Time: 0.062121 seconds (Warm-up)
               0.147839 seconds (Sampling)
               0.20996 seconds (Total)

data.rstan.mult

Inference for Stan model: e7d4e08d9f1bbcf9d6d79f126f9e56c7.
3 chains, each with iter=4000; warmup=1000; thin=3; 
post-warmup draws per chain=1000, total post-warmup draws=3000.

          mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
beta[1]   2.94    0.01 0.51   1.96   2.60   2.94   3.29   3.95  2830    1
beta[2]   1.34    0.01 0.43   0.51   1.05   1.33   1.62   2.19  2755    1
beta[3]   2.70    0.02 1.24   0.28   1.88   2.70   3.52   5.17  2701    1
cbeta0    3.82    0.00 0.11   3.60   3.75   3.82   3.90   4.05  2854    1
sigma     1.13    0.00 0.08   0.98   1.07   1.12   1.18   1.31  3000    1
beta0     3.67    0.00 0.14   3.40   3.58   3.67   3.76   3.94  2515    1
lp__    -60.99    0.03 1.67 -65.06 -61.83 -60.64 -59.76 -58.82  2502    1

Samples were drawn using NUTS(diag_e) at Tue Aug 15 15:21:02 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.

Additive model

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

data.rstanarm.add = stan_glm(y ~ cx1 + cx2, data = data, iter = 2000, warmup = 200,
    chains = 3, thin = 2, refresh = 0, prior_intercept = normal(0, 100),
    prior = normal(0, 100), prior_aux = cauchy(0, 2))

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



 Elapsed Time: 0.198223 seconds (Warm-up)
               0.678626 seconds (Sampling)
               0.876849 seconds (Total)


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



 Elapsed Time: 0.202665 seconds (Warm-up)
               0.536397 seconds (Sampling)
               0.739062 seconds (Total)


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



 Elapsed Time: 0.160257 seconds (Warm-up)
               0.56306 seconds (Sampling)
               0.723317 seconds (Total)

print(data.rstanarm.add)

stan_glm
 family:  gaussian [identity]
 formula: y ~ cx1 + cx2
------

Estimates:
            Median MAD_SD
(Intercept) 3.8    0.1   
cx1         3.0    0.5   
cx2         1.4    0.4   
sigma       1.1    0.1   

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

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

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

         term estimate  std.error  conf.low conf.high
1 (Intercept) 3.825270 0.11429998 3.6128940  4.066603
2         cx1 3.021411 0.50685407 1.9964692  3.978835
3         cx2 1.381486 0.43201840 0.5227742  2.206361
4       sigma 1.148550 0.08466719 0.9874055  1.319120

Multiplicative model

data.rstanarm.mult = stan_glm(y ~ cx1 * cx2, data = data, iter = 2000,
    warmup = 200, chains = 3, thin = 2, refresh = 0, prior_intercept = normal(0,
        100), prior = normal(0, 100), prior_aux = cauchy(0, 2))

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



 Elapsed Time: 0.112049 seconds (Warm-up)
               0.580318 seconds (Sampling)
               0.692367 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.059281 seconds (Warm-up)
               0.420321 seconds (Sampling)
               0.479602 seconds (Total)


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



 Elapsed Time: 0.210373 seconds (Warm-up)
               0.70337 seconds (Sampling)
               0.913743 seconds (Total)

print(data.rstanarm.mult)

stan_glm
 family:  gaussian [identity]
 formula: y ~ cx1 * cx2
------

Estimates:
            Median MAD_SD
(Intercept) 3.7    0.1   
cx1         2.9    0.5   
cx2         1.3    0.4   
cx1:cx2     2.7    1.2   
sigma       1.1    0.1   

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

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

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

         term estimate  std.error  conf.low conf.high
1 (Intercept) 3.668833 0.13766864 3.3846348  3.921644
2         cx1 2.933656 0.49365987 1.9342004  3.875018
3         cx2 1.335035 0.41932624 0.5170813  2.187972
4     cx1:cx2 2.703573 1.23876604 0.3719395  5.211188
5       sigma 1.128580 0.08523699 0.9605405  1.293866

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.

Additive model

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

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

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



 Elapsed Time: 0.013242 seconds (Warm-up)
               0.072156 seconds (Sampling)
               0.085398 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.012669 seconds (Warm-up)
               0.068201 seconds (Sampling)
               0.08087 seconds (Total)


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



 Elapsed Time: 0.013731 seconds (Warm-up)
               0.065669 seconds (Sampling)
               0.0794 seconds (Total)

print(data.brms.add)

 Family: gaussian(identity) 
Formula: y ~ cx1 + cx2 
   Data: data (Number of observations: 100) 
Samples: 3 chains, each with iter = 2000; warmup = 200; thin = 2; 
         total post-warmup samples = 2700
    ICs: LOO = Not computed; WAIC = Not computed
 
Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept     3.82      0.12     3.59     4.04       2700    1
cx1           3.04      0.52     2.04     4.08       2427    1
cx2           1.38      0.44     0.48     2.24       2463    1

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     1.15      0.08     0.99     1.33       2405    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.add, conf.int = TRUE, conf.method = "HPDinterval")

         term estimate  std.error  conf.low conf.high
1 b_Intercept 3.818513 0.11580124 3.6098266  4.062734
2       b_cx1 3.042654 0.51964636 2.0187823  4.042107
3       b_cx2 1.380825 0.44437714 0.5271996  2.258048
4       sigma 1.147234 0.08407413 0.9784836  1.302088

Multiplicative model

data.brms.mult = brm(y ~ cx1 * cx2, data = data, iter = 2000, warmup = 200,
    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.9e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.19 seconds.
Adjust your expectations accordingly!



 Elapsed Time: 0.023432 seconds (Warm-up)
               0.10201 seconds (Sampling)
               0.125442 seconds (Total)


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



 Elapsed Time: 0.021031 seconds (Warm-up)
               0.075103 seconds (Sampling)
               0.096134 seconds (Total)


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



 Elapsed Time: 0.019726 seconds (Warm-up)
               0.07707 seconds (Sampling)
               0.096796 seconds (Total)

print(data.brms.mult)

 Family: gaussian(identity) 
Formula: y ~ cx1 * cx2 
   Data: data (Number of observations: 100) 
Samples: 3 chains, each with iter = 2000; warmup = 200; thin = 2; 
         total post-warmup samples = 2700
    ICs: LOO = Not computed; WAIC = Not computed
 
Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept     3.67      0.13     3.42     3.94       2534    1
cx1           2.92      0.49     1.94     3.91       2700    1
cx2           1.35      0.43     0.54     2.19       2528    1
cx1:cx2       2.66      1.27     0.11     5.16       2411    1

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     1.12      0.08     0.98      1.3       2585    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.mult, conf.int = TRUE, conf.method = "HPDinterval")

         term estimate  std.error   conf.low conf.high
1 b_Intercept 3.671568 0.13285306 3.43612342  3.943035
2       b_cx1 2.923696 0.49283728 1.89910831  3.857202
3       b_cx2 1.353976 0.42606897 0.59537098  2.231514
4   b_cx1:cx2 2.659722 1.26670106 0.04200689  5.067596
5       sigma 1.123292 0.08287744 0.97447809  1.290851

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. Rather than dublicate this for both additive and multiplicative models, we will only explore the multiplicative model.

• Trace plots
  View trace plots

  library(MCMCpack)
  plot(data.mcmcpack.mult)
  

  

  
  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.mult)
  

  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     
   cx1         2        3834  3746         1.020     
   cx2         2        3650  3746         0.974     
   cx1:cx2     2        3680  3746         0.982     
   sigma2      2        3710  3746         0.990     
  

  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.mult)
  

          (Intercept)           cx1         cx2       cx1:cx2       sigma2
  Lag 0   1.000000000  1.0000000000  1.00000000  1.0000000000  1.000000000
  Lag 1   0.005947655  0.0005726118 -0.03747573 -0.0009129132  0.026538088
  Lag 5   0.004657780  0.0004847322  0.01289516 -0.0066712267 -0.001533080
  Lag 10 -0.014139624 -0.0055672420 -0.01006061  0.0089793972  0.006885994
  Lag 50 -0.005874175 -0.0037488027  0.01129924 -0.0125086637  0.018083774
  

  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. Rather than dublicate this for both additive and multiplicative models, we will only explore the multiplicative model.

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

• 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("beta0", "beta[1]", "beta[2]")
  plot(as.mcmc(data.r2jags.mult)[, 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)
   beta0    20       35610 3746          9.51     
   beta[1]  20       36200 3746          9.66     
   beta[2]  20       38660 3746         10.30     
   beta[3]  20       37410 3746          9.99     
   deviance 20       38660 3746         10.30     
   sigma    20       37410 3746          9.99     
  
  
  [[2]]
  
  Quantile (q) = 0.025
  Accuracy (r) = +/- 0.005
  Probability (s) = 0.95 
                                                  
            Burn-in  Total Lower bound  Dependence
            (M)      (N)   (Nmin)       factor (I)
   beta0    20       38030 3746         10.20     
   beta[1]  20       39950 3746         10.70     
   beta[2]  20       36800 3746          9.82     
   beta[3]  20       38030 3746         10.20     
   deviance 20       38030 3746         10.20     
   sigma    10       37410 3746          9.99     
  
  
  [[3]]
  
  Quantile (q) = 0.025
  Accuracy (r) = +/- 0.005
  Probability (s) = 0.95 
                                                  
            Burn-in  Total Lower bound  Dependence
            (M)      (N)   (Nmin)       factor (I)
   beta0    20       37410 3746          9.99     
   beta[1]  20       37410 3746          9.99     
   beta[2]  20       38660 3746         10.30     
   beta[3]  20       38030 3746         10.20     
   deviance 20       37410 3746          9.99     
   sigma    20       36800 3746          9.82     
  

  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)
  

                 beta0      beta[1]       beta[2]      beta[3]     deviance       sigma
  Lag 0    1.000000000  1.000000000  1.000000e+00 1.0000000000  1.000000000 1.000000000
  Lag 10  -0.003559426  0.008596973 -6.002651e-04 0.0009684002 -0.011764207 0.003113188
  Lag 50   0.007319792 -0.001489343 -1.080306e-02 0.0012955417 -0.008450594 0.005073925
  Lag 100 -0.015299928  0.001915914  5.393375e-05 0.0082757392  0.001525074 0.017887984
  Lag 500  0.002792247  0.015376731  8.762674e-03 0.0065778151 -0.009101738 0.010133970
  

  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. 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)
  s = as.array(data.rstan.mult)
  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.mult)
  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]        cbeta0
  Lag 0   1.000000000  1.00000000  1.000000000  1.0000000000
  Lag 1   0.006483455  0.02858605  0.033585670  0.0182224466
  Lag 5  -0.008761706  0.01281171 -0.013686626  0.0147104379
  Lag 10  0.005060245 -0.01900523 -0.015181308 -0.0212298669
  Lag 50 -0.028188168 -0.01532768  0.005026116  0.0007438316
  

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

    stan_trace(data.rstan.mult)
    

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

    raftery.diag(data.rstan.mult)
    

    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.mult)
    

    
    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.mult)
    

    
    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.mult)
    

    
    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.mult), regex_pars = "beta|sigma")
    

    

    library(bayesplot)
    mcmc_combo(as.matrix(data.rstan.mult), 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.mult), 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.mult)
  

• It is worth exploring the influence of our priors.

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)
  s = as.array(data.rstanarm.mult)
  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.mult)
  mcmc <- do.call(mcmc.list, plyr:::alply(s[, , -(length(s[1, 1, ]))], 2, as.mcmc))
  autocorr.diag(mcmc)
  

         (Intercept)         cx1          cx2      cx1:cx2
  Lag 0   1.00000000 1.000000000  1.000000000  1.000000000
  Lag 1   0.20527489 0.003647371  0.039860083 -0.020218444
  Lag 5  -0.04746531 0.008339269 -0.016665899 -0.024428764
  Lag 10  0.01484256 0.033294835  0.002029290  0.027973299
  Lag 50 -0.01554411 0.012833299 -0.005277095  0.002660752
  

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

    stan_trace(data.rstanarm.mult)
    

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

    raftery.diag(data.rstanarm.mult)
    

    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.mult)
    

    
    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.mult)
    

    
    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.mult)
    

    
    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.mult), regex_pars = "Intercept|x|sigma")
    

    

    library(bayesplot)
    mcmc_combo(as.array(data.rstanarm.mult))
    

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

    library(bayesplot)
    mcmc_dens(as.array(data.rstanarm.mult))
    

    
    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.mult, color_by = "vs", group_by = TRUE,
        facet_args = list(scales = "free_y"))
    

    Gradient evaluation took 6.5e-05 seconds
    1000 transitions using 10 leapfrog steps per transition would take 0.65 seconds.
    Adjust your expectations accordingly!
    
    
    
     Elapsed Time: 0.103395 seconds (Warm-up)
                   0.112204 seconds (Sampling)
                   0.215599 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.155166 seconds (Warm-up)
                   0.139936 seconds (Sampling)
                   0.295102 seconds (Total)
    

    
• via shinystan

  						  library(shinystan) 
  						  launch_shinystan(data.rstanarm.mult))      
  

• It is worth exploring the influence of our priors.

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.mult)
  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.mult)
  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.mult$fit)
    

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

    raftery.diag(data.brms.mult)
    

    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.mult$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.mult$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.mult$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.mult)
# generate a model matrix
newdata = data
Xmat = model.matrix(~cx1 * cx2, 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.mcmcpack.mult)
# generate a model matrix
newdata = newdata
Xmat = model.matrix(~cx1 * cx2, 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)
newdata.melt = newdata %>% gather(key = Pred, value = Value, cx1:cx2)
ggplot(newdata.melt) + geom_point(aes(y = resid, x = Value)) + facet_wrap(~Pred)

And now for studentized residuals

mcmc = as.data.frame(data.mcmcpack.mult)
# generate a model matrix
newdata = data
Xmat = model.matrix(~cx1 * cx2, 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.mult)
# generate a model matrix
Xmat = model.matrix(~cx1 * cx2, 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"])))
ggplot() + geom_density(data = NULL, aes(x = as.vector(yRep),
    fill = "Model"), alpha = 0.5) + geom_density(data = data,
    aes(x = y, fill = "Obs"), alpha = 0.5)

We can also explore the posteriors of each parameter.

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

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

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

mcmc = data.r2jags.mult$BUGSoutput$sims.matrix %>% as.data.frame %>% dplyr:::select(beta0,
    contains("beta"), sigma) %>% as.matrix
# generate a model matrix
newdata = data
Xmat = model.matrix(~cx1 * cx2, 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.mult$BUGSoutput$sims.matrix %>% as.data.frame %>% dplyr:::select(beta0,
    contains("beta"), sigma) %>% as.matrix
# generate a model matrix
newdata = newdata
Xmat = model.matrix(~cx1 * cx2, 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)
newdata.melt = newdata %>% gather(key = Pred, value = Value, cx1:cx2)
ggplot(newdata.melt) + geom_point(aes(y = resid, x = Value)) + facet_wrap(~Pred)

And now for studentized residuals

mcmc = data.r2jags.mult$BUGSoutput$sims.matrix %>% as.data.frame %>% dplyr:::select(beta0,
    contains("beta"), sigma) %>% as.matrix
# generate a model matrix
newdata = data
Xmat = model.matrix(~cx1 * cx2, 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.mult$BUGSoutput$sims.matrix %>% as.data.frame %>%
    dplyr:::select(beta0, contains("beta"), sigma) %>% as.matrix
# generate a model matrix
Xmat = model.matrix(~cx1 * cx2, 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"]))
ggplot() + geom_density(data = NULL, aes(x = as.vector(yRep),
    fill = "Model"), alpha = 0.5) + geom_density(data = data,
    aes(x = y, fill = "Obs"), alpha = 0.5)

We can also explore the posteriors of each parameter.

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

mcmc_areas(data.r2jags.mult$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.mult) %>% dplyr:::select(beta0, starts_with("beta"),
    sigma) %>% as.matrix
# generate a model matrix
newdata = data
Xmat = model.matrix(~cx1 * cx2, 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.mult) %>% dplyr:::select(beta0, starts_with("beta"),
    sigma) %>% as.matrix
# generate a model matrix
newdata = newdata
Xmat = model.matrix(~cx1 * cx2, 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)
newdata.melt = newdata %>% gather(key = Pred, value = Value, cx1:cx2)
ggplot(newdata.melt) + geom_point(aes(y = resid, x = Value)) + facet_wrap(~Pred)

And now for studentized residuals

mcmc = as.data.frame(data.rstan.mult) %>% dplyr:::select(beta0, starts_with("beta"),
    sigma) %>% as.matrix
# generate a model matrix
newdata = data
Xmat = model.matrix(~cx1 * cx2, 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.mult) %>% dplyr:::select(beta0,
    starts_with("beta"), sigma) %>% as.matrix
# generate a model matrix
Xmat = model.matrix(~cx1 * cx2, 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"]))
ggplot() + geom_density(data = NULL, aes(x = as.vector(yRep),
    fill = "Model"), alpha = 0.5) + geom_density(data = data,
    aes(x = y, fill = "Obs"), alpha = 0.5)

We can also explore the posteriors of each parameter.

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

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

Residuals can be computed directly within RSTANARM.

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

Residuals against predictors

resid = resid(data.rstanarm.mult)
dat.melt = data %>% mutate(resid = resid) %>% gather(key = Pred, value = Value, cx1:cx2)
ggplot(dat.melt) + geom_point(aes(y = resid, x = Value)) + facet_wrap(~Pred)

And now for studentized residuals

resid = resid(data.rstanarm.mult)
sresid = resid/sd(resid)
fit = fitted(data.rstanarm.mult)
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 compare draws (predictions) from the posterior (in blue) with the distributions of observed data (red) via violin plots. Violin plots are similar to boxplots except that they display more visual information about the density distribution.

y_pred = posterior_predict(data.rstanarm.mult)
newdata = data %>% cbind(t(y_pred)) %>% gather(key = "Rep", value = "Value",
    -y:-cx2)
newdata.melt = newdata %>% gather(key = "Pred", value = "Pred_val",
    cx1:cx2)
data.melt = data %>% gather(key = "Pred", value = "Pred_val",
    cx1:cx2)
ggplot(newdata.melt, aes(Value, x = Pred_val)) + geom_violin(color = "blue",
    fill = "blue", alpha = 0.5) + geom_violin(data = data.melt,
    aes(y = y, x = Pred_val), fill = "red", color = "red", alpha = 0.5) +
    facet_wrap(~Pred)

Yet another way to approach validation is to explore trends in posteriors in the context of the observed data.

## Calculate the fitted values
newdata = rbind(data.frame(cx1 = seq(min(data$cx1, na.rm = TRUE),
    max(data$cx1, na.rm = TRUE), len = 100), cx2 = 0), data.frame(cx1 = 0,
    cx2 = seq(min(data$cx2, na.rm = TRUE), max(data$cx2, na.rm = TRUE),
        len = 100)))
fit = posterior_predict(data.rstanarm.mult, newdata = newdata)
newdata = newdata %>% cbind(tidyMCMC(as.mcmc(fit), conf.int = TRUE,
    conf.method = "HPDinterval"))
newdata.melt = newdata %>% gather(key = "Pred", value = "Value",
    cx1:cx2) %>% filter(Value != 0)
data.melt = data %>% gather(key = "Pred", value = "Value", cx1:cx2)
ggplot(newdata.melt, aes(y = estimate, x = Value)) + geom_point(data = data.melt,
    aes(y = y)) + geom_line() + geom_ribbon(aes(ymin = conf.low,
    ymax = conf.high), fill = "blue", alpha = 0.3) + scale_y_continuous("Y") +
    scale_x_continuous("X") + theme_classic() + facet_wrap(~Pred)

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.mult), regex_pars = "Intercept|cx|sigma")

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

Residuals can be computed directly within BRMS. By default, the residuals and fitted extractor functions in brms return summarized versions (means, SE and credibility intervals). We are only interested in the mean (Estimate) estimates.

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

Residuals against predictors

resid = resid(data.brms.mult)[, "Estimate"]
dat.melt = data %>% mutate(resid = resid) %>% gather(key = Pred, value = Value, cx1:cx2)
ggplot(dat.melt) + geom_point(aes(y = resid, x = Value)) + facet_wrap(~Pred)

And now for studentized residuals

resid = resid(data.brms.mult)[, "Estimate"]
sresid = resid/sd(resid)
fit = fitted(data.brms.mult)[, "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 compare draws (predictions) from the posterior (in blue) with the distributions of observed data (red) via violin plots. Violin plots are similar to boxplots except that they display more visual information about the density distribution.

y_pred = posterior_predict(data.brms.mult)
newdata = data %>% cbind(t(y_pred)) %>% gather(key = "Rep", value = "Value",
    -y:-cx2)
newdata.melt = newdata %>% gather(key = "Pred", value = "Pred_val",
    cx1:cx2)
data.melt = data %>% gather(key = "Pred", value = "Pred_val",
    cx1:cx2)
ggplot(newdata.melt, aes(Value, x = Pred_val)) + geom_violin(color = "blue",
    fill = "blue", alpha = 0.5) + geom_violin(data = data.melt,
    aes(y = y, x = Pred_val), fill = "red", color = "red", alpha = 0.5) +
    facet_wrap(~Pred)

Yet another way to approach validation is to explore trends in posteriors in the context of the observed data.

## Calculate the fitted values
newdata = rbind(data.frame(cx1 = seq(min(data$cx1, na.rm = TRUE),
    max(data$cx1, na.rm = TRUE), len = 100), cx2 = 0), data.frame(cx1 = 0,
    cx2 = seq(min(data$cx2, na.rm = TRUE), max(data$cx2, na.rm = TRUE),
        len = 100)))
fit = posterior_predict(data.brms.mult, newdata = newdata)
newdata = newdata %>% cbind(tidyMCMC(as.mcmc(fit), conf.int = TRUE,
    conf.method = "HPDinterval"))
newdata.melt = newdata %>% gather(key = "Pred", value = "Value",
    cx1:cx2) %>% filter(Value != 0)
data.melt = data %>% gather(key = "Pred", value = "Value", cx1:cx2)
ggplot(newdata.melt, aes(y = estimate, x = Value)) + geom_point(data = data.melt,
    aes(y = y)) + geom_line() + geom_ribbon(aes(ymin = conf.low,
    ymax = conf.high), fill = "blue", alpha = 0.3) + scale_y_continuous("Y") +
    scale_x_continuous("X") + theme_classic() + facet_wrap(~Pred)

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.mult), regex_pars = "b_|sigma")

mcmc_areas(as.matrix(data.brms.mult), 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)
    }

}

Additive model (MCMCpack)

summary(data.mcmcpack.add)

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) 3.823 0.1140 0.001140       0.001145
cx1         3.011 0.5045 0.005045       0.005045
cx2         1.397 0.4354 0.004354       0.004419
sigma2      1.306 0.1921 0.001921       0.001968

2. Quantiles for each variable:

              2.5%   25%   50%   75% 97.5%
(Intercept) 3.6041 3.746 3.823 3.899 4.046
cx1         2.0282 2.671 3.018 3.345 4.010
cx2         0.5306 1.106 1.399 1.693 2.252
sigma2      0.9839 1.168 1.286 1.419 1.734

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

         term estimate std.error  conf.low conf.high
1 (Intercept) 3.822627 0.1139509 3.6048428  4.046831
2         cx1 3.011318 0.5044755 2.0435576  4.015099
3         cx2 1.396664 0.4353955 0.5283247  2.245911
4      sigma2 1.305520 0.1920807 0.9492838  1.677171

• when cx2 is held constant, a one unit increase in cx1 is associated with a 3.011318 change in y. That is, y increases at a rate of 3.011318 per unit increase in cx1 when standardized for cx2.
• when cx1 is held constant, a one unit increase in cx2 is associated with a 1.396664 change in y. That is, y increases at a rate of 1.396664 per unit increase in cx2 when standardized for cx1.

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.add[, 2])

[1] 0

mcmcpvalue(data.mcmcpack.add[, 3])

[1] 0.0019

With a p-value of essentially 0, we would conclude that there is almost no evidence that the slope was likely to be equal to zero, suggesting there is a relationship.

Multiplicative model (MCMCpack)

summary(data.mcmcpack.mult)

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) 3.671 0.1330 0.001330       0.001330
cx1         2.928 0.4994 0.004994       0.004994
cx2         1.342 0.4194 0.004194       0.004040
cx1:cx2     2.667 1.2398 0.012398       0.012398
sigma2      1.256 0.1843 0.001843       0.001877

2. Quantiles for each variable:

              2.5%   25%   50%   75% 97.5%
(Intercept) 3.4147 3.581 3.672 3.759 3.938
cx1         1.9555 2.597 2.929 3.260 3.910
cx2         0.4986 1.072 1.343 1.622 2.161
cx1:cx2     0.2512 1.828 2.669 3.495 5.080
sigma2      0.9460 1.125 1.238 1.369 1.660

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

         term estimate std.error  conf.low conf.high
1 (Intercept) 3.671103 0.1330200 3.4112189  3.933302
2         cx1 2.928409 0.4993563 1.9354791  3.878062
3         cx2 1.341735 0.4193631 0.4912293  2.150196
4     cx1:cx2 2.667311 1.2398248 0.3215486  5.133875
5      sigma2 1.256300 0.1843042 0.9336364  1.632568

• at the average level of cx2 (=0), a one unit increase in cx1 is associated with a 2.9284093 change in y. That is, y increases at a rate of 2.9284093 per unit increase in cx1 when standardized for cx2.
• at the average level of cx1 (=0), a one unit increase in cx2 is associated with a 1.3417354 change in y. That is, y increases at a rate of 1.3417354 per unit increase in cx2 when standardized for cx1.
• the degree to which the rate of change in response associated with a one unit change in cx1 changes over the range of cx2 (and vice versa) is 2.6673106.

The 95% confidence intervals for the interaction partial slope does not overlap with 0 implying a significant interaction between cx1 and cx2. This suggests that the nature of the relationship between y and cx1 depends on the level of cx2 (and vice versa). The estimates of the effect of cx1 are only appropriate when cx2 = 0 etc.

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.mult[, 2])

[1] 0

mcmcpvalue(data.mcmcpack.mult[, 3])

[1] 0.0022

mcmcpvalue(data.mcmcpack.mult[, 4])

[1] 0.0311

With a p-value of essentially 0, we would conclude that there is almost no evidence that the slope was likely to be equal to zero, suggesting there is a relationship.

Additive model (JAGS)

print(data.r2jags.add)

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]    3.028   0.501   2.037   2.692   3.032   3.365   4.009 1.001 15000
beta[2]    1.389   0.432   0.545   1.101   1.394   1.675   2.235 1.001 15000
beta0      3.823   0.115   3.597   3.745   3.822   3.901   4.048 1.001  8000
sigma      1.146   0.083   1.001   1.088   1.140   1.198   1.322 1.001 15000
deviance 309.526   2.883 305.925 307.427 308.888 310.918 316.903 1.002  3100

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 = 4.2 and DIC = 313.7
DIC is an estimate of expected predictive error (lower deviance is better).

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

      term   estimate  std.error    conf.low  conf.high
1    beta0   3.822567 0.11519912   3.6077546   4.057015
2  beta[1]   3.028129 0.50095089   2.0517612   4.018021
3  beta[2]   1.388584 0.43169505   0.5412069   2.229200
4 deviance 309.526490 2.88347247 305.5279931 315.076927
5    sigma   1.145686 0.08280216   0.9979548   1.318215

• when cx2 is held constant, a one unit increase in cx1 is associated with a 3.0281293 change in y. That is, y increases at a rate of 3.0281293 per unit increase in cx1 when standardized for cx2.
• when cx1 is held constant, a one unit increase in cx2 is associated with a 1.388584 change in y. That is, y increases at a rate of 1.388584 per unit increase in cx2 when standardized for cx1.

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.add$BUGSoutput$sims.matrix[, "beta[1]"])

[1] 0

mcmcpvalue(data.r2jags.add$BUGSoutput$sims.matrix[, "beta[2]"])

[1] 0.001866667

With a p-value of essentially 0, we would conclude that there is almost no evidence that the slope was likely to be equal to zero, suggesting there is a relationship.

Multiplicative model (JAGS)

print(data.r2jags.mult)

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]    2.931   0.499   1.969   2.593   2.929   3.265   3.902 1.001 15000
beta[2]    1.344   0.426   0.508   1.053   1.343   1.631   2.181 1.001 15000
beta[3]    2.675   1.256   0.182   1.845   2.669   3.509   5.158 1.001 15000
beta0      3.671   0.134   3.413   3.580   3.671   3.760   3.932 1.001 12000
sigma      1.126   0.082   0.981   1.070   1.120   1.176   1.307 1.001  8000
deviance 305.828   3.280 301.532 303.441 305.137 307.457 314.045 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 = 5.4 and DIC = 311.2
DIC is an estimate of expected predictive error (lower deviance is better).

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

      term   estimate  std.error    conf.low  conf.high
1    beta0   3.670651 0.13361307   3.4145748   3.933344
2  beta[1]   2.930546 0.49865289   1.9650810   3.893332
3  beta[2]   1.344297 0.42580446   0.5242358   2.194335
4  beta[3]   2.675152 1.25623900   0.1764777   5.139000
5 deviance 305.827579 3.27978099 301.0348960 312.277506
6    sigma   1.126052 0.08177241   0.9662253   1.287832

• at the average level of cx2 (=0), a one unit increase in cx1 is associated with a 2.9305463 change in y. That is, y increases at a rate of 2.9305463 per unit increase in cx1 when standardized for cx2.
• at the average level of cx1 (=0), a one unit increase in cx2 is associated with a 1.3442967 change in y. That is, y increases at a rate of 1.3442967 per unit increase in cx2 when standardized for cx1.
• the degree to which the rate of change in response associated with a one unit change in cx1 changes over the range of cx2 (and vice versa) is 2.6751525.

The 95% confidence intervals for the interaction partial slope does not overlap with 0 implying a significant interaction between cx1 and cx2. This suggests that the nature of the relationship between y and cx1 depends on the level of cx2 (and vice versa). The estimates of the effect of cx1 are only appropriate when cx2 = 0 etc.

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.mult$BUGSoutput$sims.matrix[, "beta[1]"])

[1] 0

mcmcpvalue(data.r2jags.mult$BUGSoutput$sims.matrix[, "beta[2]"])

[1] 0.001933333

mcmcpvalue(data.r2jags.mult$BUGSoutput$sims.matrix[, "beta[3]"])

[1] 0.0366

With a p-value of essentially 0, we would conclude that there is almost no evidence that the slope was likely to be equal to zero, suggesting there is a relationship.

Additive model (RSTAN)

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

Inference for Stan model: 4e0a54cd22440d9847d7705bdb1ff803.
3 chains, each with iter=4000; warmup=1000; thin=3; 
post-warmup draws per chain=1000, total post-warmup draws=3000.

        mean se_mean   sd 2.5%  25%  50%  75% 97.5% n_eff Rhat
beta0   3.83    0.00 0.12 3.59 3.75 3.83 3.91  4.06  2317    1
beta[1] 3.02    0.01 0.51 2.03 2.69 3.02 3.34  4.05  2747    1
beta[2] 1.38    0.01 0.43 0.52 1.10 1.39 1.68  2.22  2854    1
sigma   1.15    0.00 0.08 0.99 1.09 1.14 1.20  1.32  2806    1

Samples were drawn using NUTS(diag_e) at Thu Aug 17 15:32:58 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.add, conf.int = TRUE, conf.method = "HPDinterval",
    pars = c("beta0", "beta", "sigma"))

     term estimate  std.error  conf.low conf.high
1   beta0 3.826110 0.11727211 3.5953559  4.058548
2 beta[1] 3.022736 0.50695556 2.0124349  4.027372
3 beta[2] 1.384733 0.43249804 0.5363925  2.228436
4   sigma 1.145670 0.08415868 0.9808821  1.301693

• when cx2 is held constant, a one unit increase in cx1 is associated with a 1.3847333 change in y. That is, y increases at a rate of 1.3847333 per unit increase in cx1 when standardized for cx2.
• when cx1 is held constant, a one unit increase in cx2 is associated with a 3.82611 change in y. That is, y increases at a rate of 3.82611 per unit increase in cx2 when standardized for cx1.

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.add)[, "beta[1]"])

[1] 0

mcmcpvalue(as.matrix(data.rstan.add)[, "beta[2]"])

[1] 0.0006666667

With a p-value of essentially 0, we would conclude that there is almost no evidence that the slope was likely to be equal to zero, suggesting there is a relationship.

Multiplicative model (RSTAN)

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

Inference for Stan model: 4e0a54cd22440d9847d7705bdb1ff803.
3 chains, each with iter=4000; warmup=1000; thin=3; 
post-warmup draws per chain=1000, total post-warmup draws=3000.

        mean se_mean   sd 2.5%  25%  50%  75% 97.5% n_eff Rhat
beta0   3.67    0.00 0.14 3.40 3.58 3.67 3.76  3.94  2515    1
beta[1] 2.94    0.01 0.51 1.96 2.60 2.94 3.29  3.95  2830    1
beta[2] 1.34    0.01 0.43 0.51 1.05 1.33 1.62  2.19  2755    1
beta[3] 2.70    0.02 1.24 0.28 1.88 2.70 3.52  5.17  2701    1
sigma   1.13    0.00 0.08 0.98 1.07 1.12 1.18  1.31  3000    1

Samples were drawn using NUTS(diag_e) at Wed Aug 16 16:06:53 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.mult, conf.int = TRUE, conf.method = "HPDinterval",
    pars = c("beta0", "beta", "sigma"))

     term estimate  std.error  conf.low conf.high
1   beta0 3.667709 0.13512216 3.3976772  3.923658
2 beta[1] 2.944087 0.50916320 1.9886013  3.974140
3 beta[2] 1.335580 0.43089574 0.4927849  2.175954
4 beta[3] 2.695545 1.23745247 0.2529844  5.066406
5   sigma 1.127368 0.08492246 0.9646372  1.292453

• at the average level of cx2 (=0), a one unit increase in cx1 is associated with a 1.3355795 change in y. That is, y increases at a rate of 1.3355795 per unit increase in cx1 when standardized for cx2.
• at the average level of cx1 (=0), a one unit increase in cx2 is associated with a 2.6955452 change in y. That is, y increases at a rate of 2.6955452 per unit increase in cx2 when standardized for cx1.
• the degree to which the rate of change in response associated with a one unit change in cx1 changes over the range of cx2 (and vice versa) is 3.8218391.

The 95% confidence intervals for the interaction partial slope does not overlap with 0 implying a significant interaction between cx1 and cx2. This suggests that the nature of the relationship between y and cx1 depends on the level of cx2 (and vice versa). The estimates of the effect of cx1 are only appropriate when cx2 = 0 etc.

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.mult)[, "beta[1]"])

[1] 0

mcmcpvalue(as.matrix(data.rstan.mult)[, "beta[2]"])

[1] 0.001666667

mcmcpvalue(as.matrix(data.rstan.mult)[, "beta[3]"])

[1] 0.032

With a p-value of essentially 0, we would conclude that there is almost no evidence that the slope was likely to be equal to zero, suggesting there is a relationship.

An alternative way of quantifying the impact of an interaction is to compare models with and without the interactions. In a Bayesian context, this can be achieved by comparing the leave-one-out cross-validation statistics. Leave-one-out (LOO) cross-validation explores how well a series of models can predict withheld values (Vehtari, Gelman, and Gabry, 2016b). The LOO Information Criterion (LOOIC) is analogous to the AIC except that the LOOIC takes priors into consideration, does not assume that the posterior distribution is drawn from a multivariate normal and integrates over parameter uncertainty so as to yield a distribution of looic rather than just a point estimate. The LOOIC does however assume that all observations are equally influential (it does not matter which observations are left out). This assumption can be examined via the Pareta k estimate (values greater than 0.5 or more conservatively 0.75 are considered overly influential).

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

Computed from 3000 by 100 log-likelihood matrix

         Estimate   SE
elpd_loo   -155.3  6.2
p_loo         4.4  0.7
looic       310.5 12.4

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

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

Computed from 3000 by 100 log-likelihood matrix

         Estimate   SE
elpd_loo   -156.8  6.9
p_loo         3.8  0.7
looic       313.6 13.7

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

Additive model (RSTANARM)

summary(data.rstanarm.add)

Model Info:

 function:  stan_glm
 family:    gaussian [identity]
 formula:   y ~ cx1 + cx2
 algorithm: sampling
 priors:    see help('prior_summary')
 sample:    2700 (posterior sample size)
 num obs:   100

Estimates:
                mean   sd     2.5%   25%    50%    75%    97.5%
(Intercept)      3.8    0.1    3.6    3.7    3.8    3.9    4.1 
cx1              3.0    0.5    2.0    2.7    3.0    3.4    4.0 
cx2              1.4    0.4    0.5    1.1    1.4    1.7    2.2 
sigma            1.1    0.1    1.0    1.1    1.1    1.2    1.3 
mean_PPD         3.8    0.2    3.5    3.7    3.8    3.9    4.1 
log-posterior -165.0    1.5 -168.8 -165.8 -164.6 -163.9 -163.2 

Diagnostics:
              mcse Rhat n_eff
(Intercept)   0.0  1.0  1591 
cx1           0.0  1.0  2700 
cx2           0.0  1.0  2700 
sigma         0.0  1.0  1641 
mean_PPD      0.0  1.0  2101 
log-posterior 0.0  1.0  1335 

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.add$stanfit, conf.int = TRUE, conf.method = "HPDinterval")

           term    estimate  std.error     conf.low   conf.high
1   (Intercept)    3.825270 0.11429998    3.6128940    4.066603
2           cx1    3.021411 0.50685407    1.9964692    3.978835
3           cx2    1.381486 0.43201840    0.5227742    2.206361
4         sigma    1.148550 0.08466719    0.9874055    1.319120
5      mean_PPD    3.823862 0.15921259    3.5193260    4.137062
6 log-posterior -164.988433 1.48402493 -167.9667057 -162.954149

• when cx2 is held constant, a one unit increase in cx1 is associated with a 3.021411 change in y. That is, y increases at a rate of 3.021411 per unit increase in cx1 when standardized for cx2.
• when cx1 is held constant, a one unit increase in cx2 is associated with a 1.3814856 change in y. That is, y increases at a rate of 1.3814856 per unit increase in cx2 when standardized for cx1.

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

## since values are less than zero
mcmcpvalue(as.matrix(data.rstanarm.add)[, "cx1"])

[1] 0

mcmcpvalue(as.matrix(data.rstanarm.add)[, "cx2"])

[1] 0.001481481

With a p-value of essentially 0, we would conclude that there is almost no evidence that the slope was likely to be equal to zero, suggesting there is a relationship.

Multiplicative model (RSTANARM)

summary(data.rstanarm.mult)

Model Info:

 function:  stan_glm
 family:    gaussian [identity]
 formula:   y ~ cx1 * cx2
 algorithm: sampling
 priors:    see help('prior_summary')
 sample:    2700 (posterior sample size)
 num obs:   100

Estimates:
                mean   sd     2.5%   25%    50%    75%    97.5%
(Intercept)      3.7    0.1    3.4    3.6    3.7    3.8    3.9 
cx1              2.9    0.5    1.9    2.6    2.9    3.3    3.9 
cx2              1.3    0.4    0.5    1.1    1.3    1.6    2.2 
cx1:cx2          2.7    1.2    0.3    1.9    2.7    3.5    5.1 
sigma            1.1    0.1    1.0    1.1    1.1    1.2    1.3 
mean_PPD         3.8    0.2    3.5    3.7    3.8    3.9    4.1 
log-posterior -164.1    1.7 -168.2 -164.9 -163.7 -162.8 -161.9 

Diagnostics:
              mcse Rhat n_eff
(Intercept)   0.0  1.0  1784 
cx1           0.0  1.0  2661 
cx2           0.0  1.0  2474 
cx1:cx2       0.0  1.0  2700 
sigma         0.0  1.0  1134 
mean_PPD      0.0  1.0  2086 
log-posterior 0.0  1.0  1374 

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.mult$stanfit, conf.int = TRUE, conf.method = "HPDinterval")

           term    estimate  std.error     conf.low   conf.high
1   (Intercept)    3.668833 0.13766864    3.3846348    3.921644
2           cx1    2.933656 0.49365987    1.9342004    3.875018
3           cx2    1.335035 0.41932624    0.5170813    2.187972
4       cx1:cx2    2.703573 1.23876604    0.3719395    5.211188
5         sigma    1.128580 0.08523699    0.9605405    1.293866
6      mean_PPD    3.823219 0.15952737    3.5108288    4.130149
7 log-posterior -164.059918 1.68184862 -167.3759390 -161.738628

• at the average level of cx2 (=0), a one unit increase in cx1 is associated with a 2.9336564 change in y. That is, y increases at a rate of 2.9336564 per unit increase in cx1 when standardized for cx2.
• at the average level of cx1 (=0), a one unit increase in cx2 is associated with a 1.3350352 change in y. That is, y increases at a rate of 1.3350352 per unit increase in cx2 when standardized for cx1.
• the degree to which the rate of change in response associated with a one unit change in cx1 changes over the range of cx2 (and vice versa) is 2.7035727.

The 95% confidence intervals for the interaction partial slope does not overlap with 0 implying a significant interaction between cx1 and cx2. This suggests that the nature of the relationship between y and cx1 depends on the level of cx2 (and vice versa). The estimates of the effect of cx1 are only appropriate when cx2 = 0 etc.

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

## since values are less than zero
mcmcpvalue(as.matrix(data.rstanarm.mult)[, "cx1"])

[1] 0

mcmcpvalue(as.matrix(data.rstanarm.mult)[, "cx2"])

[1] 0.002222222

mcmcpvalue(as.matrix(data.rstanarm.mult)[, "cx1:cx2"])

[1] 0.03185185

With a p-value of essentially 0, we would conclude that there is almost no evidence that the slope was likely to be equal to zero, suggesting there is a relationship.

With a p-value of essentially 0, we would conclude that there is almost no evidence that the slope was likely to be equal to zero, suggesting there is a relationship.

Alternatively we can generate posterior intervals for each parameter.

posterior_interval(data.rstanarm.mult, prob = 0.95)

                 2.5%    97.5%
(Intercept) 3.3933953 3.940445
cx1         1.9391715 3.879278
cx2         0.5196027 2.194059
cx1:cx2     0.2548959 5.109150
sigma       0.9753235 1.314384

An alternative way of quantifying the impact of an interaction is to compare models with and without the interactions. In a Bayesian context, this can be achieved by comparing the leave-one-out cross-validation statistics. Leave-one-out (LOO) cross-validation explores how well a series of models can predict withheld values (Vehtari, Gelman, and Gabry, 2016b). The LOO Information Criterion (LOOIC) is analogous to the AIC except that the LOOIC takes priors into consideration, does not assume that the posterior distribution is drawn from a multivariate normal and integrates over parameter uncertainty so as to yield a distribution of looic rather than just a point estimate. The LOOIC does however assume that all observations are equally influential (it does not matter which observations are left out). This assumption can be examined via the Pareta k estimate (values greater than 0.5 or more conservatively 0.75 are considered overly influential).

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

Computed from 2700 by 100 log-likelihood matrix

         Estimate   SE
elpd_loo   -155.3  6.2
p_loo         4.4  0.7
looic       310.5 12.4

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

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

Computed from 2700 by 100 log-likelihood matrix

         Estimate   SE
elpd_loo   -156.8  6.8
p_loo         3.8  0.7
looic       313.6 13.7

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 
     -1.5       1.8 

Additive model (BRMS)

summary(data.brms.add)

 Family: gaussian(identity) 
Formula: y ~ cx1 + cx2 
   Data: data (Number of observations: 100) 
Samples: 3 chains, each with iter = 2000; warmup = 200; thin = 2; 
         total post-warmup samples = 2700
    ICs: LOO = Not computed; WAIC = Not computed
 
Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept     3.82      0.12     3.59     4.04       2700    1
cx1           3.04      0.52     2.04     4.08       2427    1
cx2           1.38      0.44     0.48     2.24       2463    1

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     1.15      0.08     0.99     1.33       2405    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.add$fit, conf.int = TRUE, conf.method = "HPDinterval")

         term estimate  std.error  conf.low conf.high
1 b_Intercept 3.818513 0.11580124 3.6098266  4.062734
2       b_cx1 3.042654 0.51964636 2.0187823  4.042107
3       b_cx2 1.380825 0.44437714 0.5271996  2.258048
4       sigma 1.147234 0.08407413 0.9784836  1.302088

• when cx2 is held constant, a one unit increase in cx1 is associated with a 3.0426537 change in y. That is, y increases at a rate of 3.0426537 per unit increase in cx1 when standardized for cx2.
• when cx1 is held constant, a one unit increase in cx2 is associated with a 1.3808245 change in y. That is, y increases at a rate of 1.3808245 per unit increase in cx2 when standardized for cx1.

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

## since values are less than zero
mcmcpvalue(as.matrix(data.brms.add)[, "b_cx1"])

[1] 0

mcmcpvalue(as.matrix(data.brms.add)[, "b_cx2"])

[1] 0.002592593

With a p-value of essentially 0, we would conclude that there is almost no evidence that the slope was likely to be equal to zero, suggesting there is a relationship.

Multiplicative model (BRMS)

summary(data.brms.mult)

 Family: gaussian(identity) 
Formula: y ~ cx1 * cx2 
   Data: data (Number of observations: 100) 
Samples: 3 chains, each with iter = 2000; warmup = 200; thin = 2; 
         total post-warmup samples = 2700
    ICs: LOO = Not computed; WAIC = Not computed
 
Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept     3.67      0.13     3.42     3.94       2534    1
cx1           2.92      0.49     1.94     3.91       2700    1
cx2           1.35      0.43     0.54     2.19       2528    1
cx1:cx2       2.66      1.27     0.11     5.16       2411    1

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     1.12      0.08     0.98      1.3       2585    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.mult$fit, conf.int = TRUE, conf.method = "HPDinterval")

         term estimate  std.error   conf.low conf.high
1 b_Intercept 3.671568 0.13285306 3.43612342  3.943035
2       b_cx1 2.923696 0.49283728 1.89910831  3.857202
3       b_cx2 1.353976 0.42606897 0.59537098  2.231514
4   b_cx1:cx2 2.659722 1.26670106 0.04200689  5.067596
5       sigma 1.123292 0.08287744 0.97447809  1.290851

• at the average level of cx2 (=0), a one unit increase in cx1 is associated with a 2.923696 change in y. That is, y increases at a rate of 2.923696 per unit increase in cx1 when standardized for cx2.
• at the average level of cx1 (=0), a one unit increase in cx2 is associated with a 1.3539765 change in y. That is, y increases at a rate of 1.3539765 per unit increase in cx2 when standardized for cx1.
• the degree to which the rate of change in response associated with a one unit change in cx1 changes over the range of cx2 (and vice versa) is 2.6597222.

The 95% confidence intervals for the interaction partial slope does not overlap with 0 implying a significant interaction between cx1 and cx2. This suggests that the nature of the relationship between y and cx1 depends on the level of cx2 (and vice versa). The estimates of the effect of cx1 are only appropriate when cx2 = 0 etc.

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

## since values are less than zero
mcmcpvalue(as.matrix(data.brms.mult)[, "b_cx1"])

[1] 0

mcmcpvalue(as.matrix(data.brms.mult)[, "b_cx2"])

[1] 0.0007407407

mcmcpvalue(as.matrix(data.brms.mult)[, "b_cx1:cx2"])

[1] 0.04

With a p-value of essentially 0, we would conclude that there is almost no evidence that the slope was likely to be equal to zero, suggesting there is a relationship.

With a p-value of essentially 0, we would conclude that there is almost no evidence that the slope was likely to be equal to zero, suggesting there is a relationship.

Alternatively we can generate posterior intervals for each parameter.

posterior_interval(as.matrix(data.brms.mult), prob = 0.95)

                   2.5%      97.5%
b_Intercept   3.4249254   3.935567
b_cx1         1.9378137   3.907302
b_cx2         0.5370439   2.190452
b_cx1:cx2     0.1139071   5.162096
sigma         0.9772463   1.299092
lp__        -64.9204721 -58.831482

An alternative way of quantifying the impact of an interaction is to compare models with and without the interactions. In a Bayesian context, this can be achieved by comparing the leave-one-out cross-validation statistics. Leave-one-out (LOO) cross-validation explores how well a series of models can predict withheld values (Vehtari, Gelman, and Gabry, 2016b). The LOO Information Criterion (LOOIC) is analogous to the AIC except that the LOOIC takes priors into consideration, does not assume that the posterior distribution is drawn from a multivariate normal and integrates over parameter uncertainty so as to yield a distribution of looic rather than just a point estimate. The LOOIC does however assume that all observations are equally influential (it does not matter which observations are left out). This assumption can be examined via the Pareta k estimate (values greater than 0.5 or more conservatively 0.75 are considered overly influential).

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

  LOOIC    SE
 310.51 12.48

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

  LOOIC   SE
 313.62 13.7

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.

Additive model (MCMCpack)

With appropriate use of model matrices and data wrangling, it is possible to produce a single prediction data set along with ggplot() syntax to produce a multi-panel figure.

mcmc = data.mcmcpack.add
## Calculate the fitted values
newdata = rbind(data.frame(cx1 = seq(min(data$cx1, na.rm = TRUE), max(data$cx1,
    na.rm = TRUE), len = 100), cx2 = 0, Pred = 1), data.frame(cx1 = 0,
    cx2 = seq(min(data$cx2, na.rm = TRUE), max(data$cx2, na.rm = TRUE),
        len = 100), Pred = 2))
Xmat = model.matrix(~cx1 + cx2, newdata)
coefs = mcmc[, 1:3]
fit = coefs %*% t(Xmat)
newdata = newdata %>% mutate(x1 = cx1 + mean.x1, x2 = cx2 + mean.x2) %>%
    cbind(tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval")) %>%
    mutate(x = dplyr:::recode(Pred, x1, x2))

ggplot(newdata, aes(y = estimate, x = x)) + geom_line() + geom_ribbon(aes(ymin = conf.low,
    ymax = conf.high), fill = "blue", alpha = 0.3) + scale_y_continuous("Y") +
    scale_x_continuous("X") + theme_classic() + facet_wrap(~Pred)

We cannot simply add the raw data to this figure. The reason for this is that the trends represent the effect of one predictor holding the other variable constant. Therefore, the observations we represent on the figure must likewise be standardized. We can achieve this by adding the partial residuals to the figure. Partial residuals are the fitted values plus the residuals.

## Calculate partial residuals fitted values
fdata = rdata = rbind(data.frame(cx1 = data$cx1, cx2 = 0, Pred = 1), data.frame(cx1 = 0,
    cx2 = data$cx2, Pred = 2))
fMat = rMat = model.matrix(~cx1 + cx2, 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) %>% mutate(x1 = cx1 +
    mean.x1, x2 = cx2 + mean.x2) %>% mutate(x = dplyr:::recode(Pred, x1,
    x2))

ggplot(newdata, aes(y = estimate, x = x)) + geom_point(data = rdata, aes(y = partial.resid),
    color = "gray") + geom_line() + geom_ribbon(aes(ymin = conf.low, ymax = conf.high),
    fill = "blue", alpha = 0.3) + scale_y_continuous("Y") + theme_classic() +
    facet_wrap(~Pred, strip.position = "bottom", labeller = label_bquote("x" *
        .(Pred))) + theme(axis.title.x = element_blank(), strip.background = element_blank(),
    strip.placement = "outside")

However, this method (whist partially elegant) does become overly opaque if we need more extensive axes labels since the x-axes labels are actually strip labels (which must largely be defined outside of the ggplot structure. The alternative is to simply produce each partial plot separately before arranging them together in the one figure.

mcmc = data.mcmcpack.add
## Calculate the fitted values
newdata = data.frame(cx1 = seq(min(data$cx1, na.rm = TRUE), max(data$cx1,
    na.rm = TRUE), len = 100), cx2 = 0)
Xmat = model.matrix(~cx1 + cx2, newdata)
coefs = mcmc[, 1:3]
fit = coefs %*% t(Xmat)
newdata = newdata %>% mutate(x1 = cx1 + mean.x1, x2 = cx2 + mean.x2) %>%
    cbind(tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval"))
## Now the partial residuals
fdata = rdata = data.frame(cx1 = data$cx1, cx2 = 0)
fMat = rMat = model.matrix(~cx1 + cx2, 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) %>% mutate(x1 = cx1 +
    mean.x1, x2 = cx2 + mean.x2)
g1 = ggplot(newdata, aes(y = estimate, x = x1)) + geom_point(data = rdata,
    aes(y = partial.resid), color = "grey") + geom_line() + geom_ribbon(aes(ymin = conf.low,
    ymax = conf.high), fill = "blue", alpha = 0.3) + scale_y_continuous("Y") +
    scale_x_continuous("X1") + theme_classic()

newdata = data.frame(cx2 = seq(min(data$cx2, na.rm = TRUE), max(data$cx2,
    na.rm = TRUE), len = 100), cx1 = 0)
Xmat = model.matrix(~cx1 + cx2, newdata)
coefs = mcmc[, 1:3]
fit = coefs %*% t(Xmat)
newdata = newdata %>% mutate(x1 = cx1 + mean.x1, x2 = cx2 + mean.x2) %>%
    cbind(tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval"))
## Now the partial residuals
fdata = rdata = data.frame(cx1 = 0, cx2 = data$cx2)
fMat = rMat = model.matrix(~cx1 + cx2, 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) %>% mutate(x1 = cx1 +
    mean.x1, x2 = cx2 + mean.x2)
g2 = ggplot(newdata, aes(y = estimate, x = x2)) + geom_point(data = rdata,
    aes(y = partial.resid), color = "grey") + geom_line() + geom_ribbon(aes(ymin = conf.low,
    ymax = conf.high), fill = "blue", alpha = 0.3) + scale_y_continuous("Y") +
    scale_x_continuous("X2") + theme_classic()

grid.arrange(g1, g2, ncol = 2)

Multiplicative model (MCMCpack)

For the multiplicative model, we could elect to split the trends up so as to explore the effects of one predictor at several set levels of another predictor. In this example, we will explore the effects of x1 when x2 is equal to its mean in the original data as well as one and two standard deviations below and above this mean.

mcmc = data.mcmcpack.mult
## Calculate the fitted values
newdata = expand.grid(cx1 = seq(min(data$cx1, na.rm = TRUE), max(data$cx1,
    na.rm = TRUE), len = 100), cx2 = mean(data$cx2) + sd(data$cx2) %*%
    -2:2)
Xmat = model.matrix(~cx1 * cx2, newdata)
coefs = mcmc[, 1:4]
fit = coefs %*% t(Xmat)
newdata = newdata %>% mutate(x1 = cx1 + mean.x1, x2 = cx2 + mean.x2) %>%
    cbind(tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval")) %>%
    mutate(x2 = factor(x2, labels = paste("X2:~", c(-2, -1, 0, 1, 2), "*sigma")))
## Partial residuals
fdata = rdata = expand.grid(cx1 = data$cx1, cx2 = mean(data$cx2) + sd(data$cx2) *
    -2:2)
fMat = rMat = model.matrix(~cx1 * cx2, 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) %>% mutate(x1 = cx1 +
    mean.x1, x2 = cx2 + mean.x2)
## Partition the partial residuals such that each x1 trend only includes
## x2 data that is within that range in the observed data
findNearest = function(x, y) {
    ff = fields:::rdist(x, y)
    apply(ff, 1, function(x) which(x == min(x)))
}
fn = findNearest(x = data[, c("x1", "x2")], y = rdata[, c("x1", "x2")])
rdata = rdata[fn, ] %>% mutate(x2 = factor(x2, labels = paste("X2:~", c(-2,
    -1, 0, 1, 2), "*sigma")))
ggplot(newdata, aes(y = estimate, x = x1)) + geom_line() + geom_blank(aes(y = 9)) +
    geom_point(data = rdata, aes(y = partial.resid), color = "grey") +
    geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "blue",
        alpha = 0.3) + scale_y_continuous("Y") + scale_x_continuous("X1") +
    facet_wrap(~x2, labeller = label_parsed, nrow = 1, scales = "free_y") +
    theme_classic() + theme(strip.background = element_blank())

Alternatively, we could explore the interaction by plotting a two dimensional surface as a heat map. In the example below, relative confidence (upper - lower) is indicated via contours.

mcmc = data.mcmcpack.mult
## Calculate the fitted values
newdata = expand.grid(cx1 = seq(min(data$cx1, na.rm = TRUE), max(data$cx1,
    na.rm = TRUE), len = 100), cx2 = seq(min(data$cx2, na.rm = TRUE), max(data$cx2,
    na.rm = TRUE), len = 100))
Xmat = model.matrix(~cx1 * cx2, newdata)
coefs = mcmc[, 1:4]
fit = coefs %*% t(Xmat)
newdata = newdata %>% mutate(x1 = cx1 + mean.x1, x2 = cx2 + mean.x2) %>%
    cbind(tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval"))

## Partial residuals
fdata = rdata = expand.grid(cx1 = data$cx1, cx2 = mean(data$cx2) + sd(data$cx2) *
    -2:2)
fMat = rMat = model.matrix(~cx1 * cx2, 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) %>% mutate(x1 = cx1 +
    mean.x1, x2 = cx2 + mean.x2)

ggplot(newdata, aes(y = x1, x = x2)) + geom_tile(aes(fill = estimate)) +
    geom_contour(aes(z = conf.high - conf.low)) + scale_fill_gradientn("Y",
    colors = heat.colors(10)) + geom_point(data = data, aes(size = y)) +
    scale_y_continuous("X1") + scale_x_continuous("X2") + theme_classic()

Additive model (JAGS)

With appropriate use of model matrices and data wrangling, it is possible to produce a single prediction data set along with ggplot() syntax to produce a multi-panel figure.

mcmc = data.r2jags.add$BUGSoutput$sims.matrix
## Calculate the fitted values
newdata = rbind(data.frame(cx1 = seq(min(data$cx1, na.rm = TRUE), max(data$cx1,
    na.rm = TRUE), len = 100), cx2 = 0, Pred = 1), data.frame(cx1 = 0,
    cx2 = seq(min(data$cx2, na.rm = TRUE), max(data$cx2, na.rm = TRUE),
        len = 100), Pred = 2))
Xmat = model.matrix(~cx1 + cx2, newdata)
coefs = mcmc[, c("beta0", "beta[1]", "beta[2]")]
fit = coefs %*% t(Xmat)
newdata = newdata %>% mutate(x1 = cx1 + mean.x1, x2 = cx2 + mean.x2) %>%
    cbind(tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval")) %>%
    mutate(x = dplyr:::recode(Pred, x1, x2))

ggplot(newdata, aes(y = estimate, x = x)) + geom_line() + geom_ribbon(aes(ymin = conf.low,
    ymax = conf.high), fill = "blue", alpha = 0.3) + scale_y_continuous("Y") +
    scale_x_continuous("X") + theme_classic() + facet_wrap(~Pred)

We cannot simply add the raw data to this figure. The reason for this is that the trends represent the effect of one predictor holding the other variable constant. Therefore, the observations we represent on the figure must likewise be standardized. We can achieve this by adding the partial residuals to the figure. Partial residuals are the fitted values plus the residuals.

## Calculate partial residuals fitted values
fdata = rdata = rbind(data.frame(cx1 = data$cx1, cx2 = 0, Pred = 1), data.frame(cx1 = 0,
    cx2 = data$cx2, Pred = 2))
fMat = rMat = model.matrix(~cx1 + cx2, 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) %>% mutate(x1 = cx1 +
    mean.x1, x2 = cx2 + mean.x2) %>% mutate(x = dplyr:::recode(Pred, x1,
    x2))

ggplot(newdata, aes(y = estimate, x = x)) + geom_point(data = rdata, aes(y = partial.resid),
    color = "gray") + geom_line() + geom_ribbon(aes(ymin = conf.low, ymax = conf.high),
    fill = "blue", alpha = 0.3) + scale_y_continuous("Y") + theme_classic() +
    facet_wrap(~Pred, strip.position = "bottom", labeller = label_bquote("x" *
        .(Pred))) + theme(axis.title.x = element_blank(), strip.background = element_blank(),
    strip.placement = "outside")

However, this method (whist partially elegant) does become overly opaque if we need more extensive axes labels since the x-axes labels are actually strip labels (which must largely be defined outside of the ggplot structure. The alternative is to simply produce each partial plot separately before arranging them together in the one figure.

mcmc = data.r2jags.add$BUGSoutput$sims.matrix
## Calculate the fitted values
newdata = data.frame(cx1 = seq(min(data$cx1, na.rm = TRUE), max(data$cx1,
    na.rm = TRUE), len = 100), cx2 = 0)
Xmat = model.matrix(~cx1 + cx2, newdata)
coefs = mcmc[, c("beta0", "beta[1]", "beta[2]")]
fit = coefs %*% t(Xmat)
newdata = newdata %>% mutate(x1 = cx1 + mean.x1, x2 = cx2 + mean.x2) %>%
    cbind(tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval"))
## Now the partial residuals
fdata = rdata = data.frame(cx1 = data$cx1, cx2 = 0)
fMat = rMat = model.matrix(~cx1 + cx2, 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) %>% mutate(x1 = cx1 +
    mean.x1, x2 = cx2 + mean.x2)
g1 = ggplot(newdata, aes(y = estimate, x = x1)) + geom_point(data = rdata,
    aes(y = partial.resid), color = "grey") + geom_line() + geom_ribbon(aes(ymin = conf.low,
    ymax = conf.high), fill = "blue", alpha = 0.3) + scale_y_continuous("Y") +
    scale_x_continuous("X1") + theme_classic()

newdata = data.frame(cx2 = seq(min(data$cx2, na.rm = TRUE), max(data$cx2,
    na.rm = TRUE), len = 100), cx1 = 0)
Xmat = model.matrix(~cx1 + cx2, newdata)
coefs = mcmc[, c("beta0", "beta[1]", "beta[2]")]
fit = coefs %*% t(Xmat)
newdata = newdata %>% mutate(x1 = cx1 + mean.x1, x2 = cx2 + mean.x2) %>%
    cbind(tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval"))
## Now the partial residuals
fdata = rdata = data.frame(cx1 = 0, cx2 = data$cx2)
fMat = rMat = model.matrix(~cx1 + cx2, 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) %>% mutate(x1 = cx1 +
    mean.x1, x2 = cx2 + mean.x2)
g2 = ggplot(newdata, aes(y = estimate, x = x2)) + geom_point(data = rdata,
    aes(y = partial.resid), color = "grey") + geom_line() + geom_ribbon(aes(ymin = conf.low,
    ymax = conf.high), fill = "blue", alpha = 0.3) + scale_y_continuous("Y") +
    scale_x_continuous("X2") + theme_classic()

grid.arrange(g1, g2, ncol = 2)

Multiplicative model (JAGS)

For the multiplicative model, we could elect to split the trends up so as to explore the effects of one predictor at several set levels of another predictor. In this example, we will explore the effects of x1 when x2 is equal to its mean in the original data as well as one and two standard deviations below and above this mean.

mcmc = data.r2jags.mult$BUGSoutput$sims.matrix
## Calculate the fitted values
newdata = expand.grid(cx1 = seq(min(data$cx1, na.rm = TRUE), max(data$cx1,
    na.rm = TRUE), len = 100), cx2 = mean(data$cx2) + sd(data$cx2) %*%
    -2:2)
Xmat = model.matrix(~cx1 * cx2, newdata)
coefs = mcmc[, c("beta0", "beta[1]", "beta[2]", "beta[3]")]
fit = coefs %*% t(Xmat)
newdata = newdata %>% mutate(x1 = cx1 + mean.x1, x2 = cx2 + mean.x2) %>%
    cbind(tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval")) %>%
    mutate(x2 = factor(x2, labels = paste("X2:~", c(-2, -1, 0, 1, 2), "*sigma")))
## Partial residuals
fdata = rdata = expand.grid(cx1 = data$cx1, cx2 = mean(data$cx2) + sd(data$cx2) *
    -2:2)
fMat = rMat = model.matrix(~cx1 * cx2, 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) %>% mutate(x1 = cx1 +
    mean.x1, x2 = cx2 + mean.x2)
## Partition the partial residuals such that each x1 trend only includes
## x2 data that is within that range in the observed data
findNearest = function(x, y) {
    ff = fields:::rdist(x, y)
    apply(ff, 1, function(x) which(x == min(x)))
}
fn = findNearest(x = data[, c("x1", "x2")], y = rdata[, c("x1", "x2")])
rdata = rdata[fn, ] %>% mutate(x2 = factor(x2, labels = paste("X2:~", c(-2,
    -1, 0, 1, 2), "*sigma")))
ggplot(newdata, aes(y = estimate, x = x1)) + geom_line() + geom_blank(aes(y = 9)) +
    geom_point(data = rdata, aes(y = partial.resid), color = "grey") +
    geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "blue",
        alpha = 0.3) + scale_y_continuous("Y") + scale_x_continuous("X1") +
    facet_wrap(~x2, labeller = label_parsed, nrow = 1, scales = "free_y") +
    theme_classic() + theme(strip.background = element_blank())