Tutorial 7.2b - Simple linear regression (Bayesian)

The purpose of fitting a model in this case is to explore the relationship between y and x. Since both y and x are continuous, a simple regression line is a good start.

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 $x$ is equal to zero) and $\beta$ represents the rate of change in $y$ for every unit change in $x$ (the effect).

Note that in this form, the y-intercept is of little interest. Indeed for many applications, a value of $x$ would be outside the domain of the collected data, outside the logical bounds of the actual variable or else outside the domain of interest. If however, we center the predictor variable (by subtracting the mean of $x$ from each $x$, then the y-intercept represents the value of $y$ at the average value of $x$. This certainly has more meaning. Note that centering the predictor does not effect the estimate of slope.

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 (1000) for both the intercept and the treatment effect and a wide half-cauchy (scale=25) for the standard deviation. $$ \begin{align} y_i &\sim{} N(\mu, \sigma)\\ \mu &= \beta_0 + \beta x_i\\[1em] \beta_0 &\sim{} N(0,1000)\\ \beta &\sim{} N(0,1000)\\ \sigma &\sim{} cauchy(0,25)\\ \end{align} $$

Structure of the JAGS model

Define the model

We now translate the likelihood model into BUGS/JAGS code and store the code in an external file.
$y_i\sim{}N(\mu_i, \tau)\\ \mu_i = \beta_0+\beta_1 x_i\\ \beta_0\sim{}N(0,1.0{E-6}) \hspace{1cm}\mathsf{non-informative~prior}\\ \beta_1\sim{}N(0,1.0{E-6}) \hspace{1cm}\mathsf{non-informative~prior}\\ \tau = 1/\sigma^2\\ \sigma\sim{}U(0,100)\\ $

• the percentage decline ($100*\frac{((max(x)-min(x))*\beta)+min(y)}{min{y}}$) and the probability that $y$ decline by more than 25%.
• the finite-population variance components

modelString = "
  model {
  #Likelihood
  for (i in 1:n) {
  y[i]~dnorm(mu[i],tau)
  mu[i] <- beta0+beta1*x[i]
  y.err[i] <- y[i] - mu[i]
  }
  
  #Priors
  beta0 ~ dnorm(0.01,1.0E-6)
  beta1 ~ dnorm(0,1.0E-6)
  tau <- 1 / (sigma * sigma)
  sigma~dunif(0,100)
  }
  "
## write the model to <a href=''></a> text file (I
## suggest you alter the path to somewhere more
## relevant to your system!)
writeLines(modelString, con = "../downloads/BUGSscripts/tut7.2bS4.1.txt")

modelString = "
  model {
  #Likelihood
  for (i in 1:n) {
  y[i]~dnorm(mu[i],tau)
  mu[i] <- beta0+beta1*x[i]
  y.err[i] <- y[i] - mu[i]
  }
  
  #Priors
  beta0 ~ dnorm(0.01,1.0E-6)
  beta1 ~ dnorm(0,1.0E-6)
  tau <- 1 / (sigma * sigma)
  sigma~dunif(0,100)
  
  #Other Derived parameters 
  p.decline <- 1-step(beta1)
  ymin<-beta0+beta1*min(x)                  
  xrange <- max(x) - min(x)       
  decline <- 100*((xrange*beta1)+ymin)/ymin 
  p.decline25 <- step(decline-25)
  
  #finite-population variance components
  sd.x <- abs(beta1)*sd(x[])
  sd.resid <- sd(y.err)
  }
  "
## write the model to <a href=''></a> text file (I
## suggest you alter the path to somewhere more
## relevant to your system!)
writeLines(modelString, con = "../downloads/BUGSscripts/tut7.2bS4.1a.txt")

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 variable (x)
• the total number of observed items (n)

data.list <- with(data, list(y = y, x = x, n = nrow(data)))
data.list

$y
 [1] 35.367731 37.918217 31.321857 41.976404
 [5] 34.147539 26.897658 31.937145 31.691624
 [9] 29.378907 23.473058 31.058906 23.949216
[13] 17.393797  7.926501 23.124655 15.775332

$x
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
[16] 16

$n
[1] 16

Define the initial values

Define the initial values ($\beta_0$, $\beta_1$ and $\sigma^2$ for the chain. Reasonable starting points can be gleaned from the data themselves. Note, this step is not absolutely necessary for simple models as the R2jags interface will automatically create sensible initial values for the parameters based on simple data summaries.

inits <- rep(list(list(beta0 = mean(data$y), beta1 = diff(tapply(data$y,
    data$x, mean)), sigma = sd(data$y))), 3)

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", "beta1", "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 <- jags(data = data.list, inits = NULL, parameters.to.save = params, model.file = "../downloads/BUGSscripts/tut7.2bS4.1.txt",
    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: 16
   Unobserved stochastic nodes: 3
   Total graph size: 119

Initializing model

print(data.r2jags)

Inference for Bugs model at "../downloads/BUGSscripts/tut7.2bS4.1.txt", 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
beta0     40.776   3.012 34.771 38.854 40.754  42.700  46.732 1.001  6500
beta1     -1.536   0.309 -2.155 -1.739 -1.534  -1.339  -0.917 1.001  8800
sigma      5.612   1.198  3.825  4.765  5.433   6.244   8.454 1.001 15000
deviance  98.913   2.916 95.586 96.795 98.131 100.243 106.547 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 = 4.3 and DIC = 103.2
DIC is an estimate of expected predictive error (lower deviance is better).

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

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

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

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

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

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

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

Structure of a stan model

Note the following important characteristics of stan code:

• A stan model file comprises a number of blocks (not all of which are compulsory).
• The stan language is an intermediary between (R/BUGS and c++), stan requires all types (integers, vectors, matrices etc) to be declared prior to use and it uses c++ commenting (// and /* */)
• Code order is important, objects must be declared before they are used. When a type is declared in one block, it is available in subsequent blocks.

					  data {
					  // declare the input data / parameters
					  }
					  transformed data {
					  // optional - for transforming/scaling input data
					  }
					  parameters {
					  // define model parameters
					  }
					  transformed parameters {
					  // optional - for deriving additional non-model parameters
					  //            note however, as they are part of the sampling chain
					  //            transformed parameters slow sampling down.
					  }
					  model {
					  // specifying priors and likelihood as well as the linear predictor
					  }
					  generated quantities {
					  // optional - derivatives (posteriors) of the samples
					  }
				  

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.
$y_i\sim{}N(\mu_i, \sigma)\\ \mu_i = \beta_0+\beta_1 x_i\\ \beta_0\sim{}N(0,1000)\\ \beta_1\sim{}N(0,1000)\\ \sigma\sim{}Cauchy(0,5)\\ $

modelString = "
  data {
  int<lower=0> n;
  vector [n] y;
  vector [n] x;
  }
  parameters {
  real beta0;
  real beta;
  real<lower=0> sigma;
  }
  model {
  vector [n] mu;
  #Priors
  beta0 ~ normal(0,10000);
  beta ~ normal(0,10000);
  sigma ~ cauchy(0,5);
 
  mu = beta0+beta*x;
  
  #Likelihood
  y~normal(mu,sigma);
  }
  "
## write the model to a text file (I suggest you alter the path to somewhere more relevant to your
## system!)
writeLines(modelString, con = "../downloads/BUGSscripts/tut7.2bS11.1.txt")

modelString1 = " 
  data {
  int<lower=0> n;
  real y[n];
  real x[n];
  }
  parameters {
  real beta0;
  real beta;
  real<lower=0> sigma;
  } 
  transformed parameters {
  real mu[n];
  for (i in 1:n)
  mu[i] = beta0+beta*x[i];
  }
  model {
  #Likelihood
  y~normal(mu,sigma);
  
  #Priors
  beta0 ~ normal(0,10000);
  beta ~ normal(0,10000);
  sigma~uniform(0,100);
  }
  "

The No-U-Turn sampler operates much more efficiently if all predictors are centered. Although it is possible to pre-center all predictors that are passed to STAN, it is then often necessary to later convert back to the original scale for graphing and further analyses. Since centering is a routine procedure, arguably it should be built into the STAN we generate. Furthermore, we should also include the back-scaling as well. The following code is inspired by the code generated by the BRMS package.

					  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 
					  beta0 = cbeta0 - dot_product(means_X, beta); 
					  }

In this version, the data are to be supplied as a model matrix (so as to leverage various vectorized and matrix multiplier routines). The transformed data block is used to center the non-intercept columns of the predictor model matrix. The model is fit on centered data thereby generating a slope and intercept. This intercept parameter is also expressed back on the non-centered scale (generated properties block).

Define the data list

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

• the response variable (y)
• the predictor variable (x)
• the total number of observed items (n)

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

$Y
 [1] 35.367731 37.918217 31.321857 41.976404 34.147539 26.897658 31.937145 31.691624 29.378907
[10] 23.473058 31.058906 23.949216 17.393797  7.926501 23.124655 15.775332

$X
   (Intercept)  x
1            1  1
2            1  2
3            1  3
4            1  4
5            1  5
6            1  6
7            1  7
8            1  8
9            1  9
10           1 10
11           1 11
12           1 12
13           1 13
14           1 14
15           1 15
16           1 16
attr(,"assign")
[1] 0 1

$nX
[1] 2

$n
[1] 16

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 <- 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 file6ff34d65b11c.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 'e7d4e08d9f1bbcf9d6d79f126f9e56c7' NOW (CHAIN 1).

Gradient evaluation took 1.4e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.14 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.018545 seconds (Warm-up)
               0.045305 seconds (Sampling)
               0.06385 seconds (Total)


SAMPLING FOR MODEL 'e7d4e08d9f1bbcf9d6d79f126f9e56c7' NOW (CHAIN 2).

Gradient evaluation took 5e-06 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.05 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.017913 seconds (Warm-up)
               0.04733 seconds (Sampling)
               0.065243 seconds (Total)


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

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


Iteration:    1 / 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.018588 seconds (Warm-up)
               0.043631 seconds (Sampling)
               0.062219 seconds (Total)

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 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, 1000)$
• weakly informative Gaussian prior for the treatment effect $\beta_1 \sim{} N(0, 1000)$
• half-cauchy prior for the variance $\sigma \sim{} Cauchy(0, 25)$

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

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

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: 1.80628 seconds (Warm-up)
               0.807009 seconds (Sampling)
               2.61329 seconds (Total)


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



 Elapsed Time: 1.98795 seconds (Warm-up)
               1.66692 seconds (Sampling)
               3.65488 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: 1.52452 seconds (Warm-up)
               0.845179 seconds (Sampling)
               2.36969 seconds (Total)

print(data.rstanarm)

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

Estimates:
            Median MAD_SD
(Intercept) 40.8    2.9  
x           -1.5    0.3  
sigma        5.4    1.1  

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

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

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

         term  estimate std.error  conf.low  conf.high
1 (Intercept) 40.734015 3.0668570 34.412462 46.3425366
2           x -1.531248 0.3148917 -2.123538 -0.8981083
3       sigma  5.619764 1.1740461  3.604894  7.9675118

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, 1000)$
• weakly informative Gaussian prior for the treatment effect $\beta_1 \sim{} N(0, 1000)$
• half-cauchy prior for the variance $\sigma \sim{} Cauchy(0, 25)$

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

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

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



 Elapsed Time: 0.006342 seconds (Warm-up)
               0.02702 seconds (Sampling)
               0.033362 seconds (Total)


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



 Elapsed Time: 0.006317 seconds (Warm-up)
               0.026596 seconds (Sampling)
               0.032913 seconds (Total)


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



 Elapsed Time: 0.005054 seconds (Warm-up)
               0.028559 seconds (Sampling)
               0.033613 seconds (Total)

print(data.brms)

 Family: gaussian(identity) 
Formula: y ~ x 
   Data: data (Number of observations: 16) 
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    40.85      3.07    34.90    47.24       2310    1
x            -1.55      0.32    -2.18    -0.90       2408    1

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     5.63      1.19      3.9     8.48       2091    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).

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

         term  estimate std.error  conf.low  conf.high
1 b_Intercept 40.853729 3.0714471 34.407211 46.6638946
2         b_x -1.545095 0.3205355 -2.141672 -0.8693169
3       sigma  5.631356 1.1888245  3.779939  8.2153548

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

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

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

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

• Trace plots
  View trace plots

  library(MCMCpack)
  plot(data.mcmcpack)
  

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

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

  Quantile (q) = 0.025
  Accuracy (r) = +/- 0.005
  Probability (s) = 0.95 
                                                     
               Burn-in  Total Lower bound  Dependence
               (M)      (N)   (Nmin)       factor (I)
   (Intercept) 3        4028  3746         1.080     
   x           2        3851  3746         1.030     
   sigma2      2        3680  3746         0.982     
  

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

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

          (Intercept)            x       sigma2
  Lag 0   1.000000000  1.000000000  1.000000000
  Lag 1   0.011287106 -0.000965701  0.162987189
  Lag 5  -0.004110531  0.004773571  0.008242103
  Lag 10 -0.015174763 -0.012219328  0.002405411
  Lag 50 -0.003762082 -0.007068066 -0.009647390
  

  A lag of 1 appears to be mainly sufficient to avoid autocorrelation (except for sigma², which is over 0.1). The diagnostic suggests that a lag of 5 would potentially correct this.

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

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

• Trace plots

  plot(data.mcmc)
  

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

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

  preds <- c("beta0", "beta1")
  plot(as.mcmc(data.r2jags)[, preds])
  

  
• Raftery diagnostic

  raftery.diag(data.mcmc)
  

  [[1]]
  
  Quantile (q) = 0.025
  Accuracy (r) = +/- 0.005
  Probability (s) = 0.95 
                                                  
            Burn-in  Total Lower bound  Dependence
            (M)      (N)   (Nmin)       factor (I)
   beta0    20       38660 3746         10.30     
   beta1    20       36200 3746          9.66     
   deviance 20       36800 3746          9.82     
   sigma    20       36200 3746          9.66     
  
  
  [[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       37410 3746          9.99     
   beta1    20       38660 3746         10.30     
   deviance 20       36800 3746          9.82     
   sigma    20       36800 3746          9.82     
  
  
  [[3]]
  
  Quantile (q) = 0.025
  Accuracy (r) = +/- 0.005
  Probability (s) = 0.95 
                                                  
            Burn-in  Total Lower bound  Dependence
            (M)      (N)   (Nmin)       factor (I)
   beta0    20       38030 3746         10.20     
   beta1    20       38030 3746         10.20     
   deviance 20       36200 3746          9.66     
   sigma    20       35610 3746          9.51     
  

  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        beta1     deviance         sigma
  Lag 0    1.000000000  1.000000000  1.000000000  1.0000000000
  Lag 10   0.002343329 -0.006499521 -0.005759350  0.0008347558
  Lag 50  -0.009224921 -0.010877289 -0.003269510  0.0022530505
  Lag 100 -0.004536817  0.002460483  0.007360297  0.0108817122
  Lag 500 -0.002154409  0.003107960 -0.004786948 -0.0052312783
  

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

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

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

• via coda
  • Traceplots

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

              beta[1]       cbeta0        sigma       beta0
  Lag 0   1.000000000  1.000000000  1.000000000  1.00000000
  Lag 1   0.016316577  0.020651289  0.036121993  0.00965383
  Lag 5  -0.006799616 -0.039385336 -0.006821465 -0.02197035
  Lag 10  0.012913933  0.006045351 -0.018816878  0.01012661
  Lag 50 -0.028343488  0.015861053  0.007679815 -0.00919722
  

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

    stan_trace(data.rstan)
    

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

    raftery.diag(data.rstan)
    

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

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

    stan_ac(data.rstan)
    

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

    stan_rhat(data.rstan)
    

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

    stan_ess(data.rstan)
    

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

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

    

    library(bayesplot)
    mcmc_combo(as.array(data.rstan))
    

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

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

  						  library(shinystan)    
  						  launch_shinystan(data.rstan))      
  

• It is worth exploring the influence of our priors.

Again, prior to examining the summaries, we should have explored the convergence diagnostics. There are numerous ways of working with 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)
  mcmc <- do.call(mcmc.list, plyr:::alply(s[, , -(length(s[1, 1, ]))], 2, as.mcmc))
  plot(mcmc)
  

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

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

          (Intercept)           x
  Lag 0   1.000000000  1.00000000
  Lag 1   0.053211483 -0.02917861
  Lag 5   0.004975124  0.01941260
  Lag 10 -0.033296048 -0.04010251
  Lag 50 -0.032866968 -0.01348437
  

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

    stan_trace(data.rstanarm)
    

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

    raftery.diag(data.rstanarm)
    

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

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

    stan_ac(data.rstanarm)
    

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

    stan_rhat(data.rstanarm)
    

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

    stan_ess(data.rstanarm)
    

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

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

    

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

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

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

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

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

    Gradient evaluation took 2.8e-05 seconds
    1000 transitions using 10 leapfrog steps per transition would take 0.28 seconds.
    Adjust your expectations accordingly!
    
    
    
     Elapsed Time: 1.96278 seconds (Warm-up)
                   0.039246 seconds (Sampling)
                   2.00202 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: 2.48823 seconds (Warm-up)
                   0.054196 seconds (Sampling)
                   2.54242 seconds (Total)
    

    
• via shinystan

  						  library(shinystan) 
  						  launch_shinystan(data.rstanarm))      
  

• It is worth exploring the influence of our priors.

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

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

• via coda
  • Traceplots

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

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

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

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

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

    stan_trace(data.brms$fit)
    

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

    raftery.diag(data.brms)
    

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

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

    stan_ac(data.brms$fit)
    

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

    stan_rhat(data.brms$fit)
    

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

    stan_ess(data.brms$fit)
    

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

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

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

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

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

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

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

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

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

Residuals against predictors

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

And now for studentized residuals

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

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

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

mcmc = as.matrix(data.mcmcpack)
# generate a model matrix
Xmat = model.matrix(~x, data)
## get median parameter estimates
coefs = mcmc[, 1:2]
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)

Although residuals can be computed directly within R2jags, we can calculate them manually from the posteriors to be consistent across other approaches.

mcmc = data.r2jags$BUGSoutput$sims.matrix[, c("beta0", "beta1")]
# generate a model matrix
newdata = data.frame(x = data$x)
Xmat = model.matrix(~x, newdata)
## get median parameter estimates
coefs = apply(mcmc, 2, median)
fit = as.vector(coefs %*% t(Xmat))
resid = data$y - fit
ggplot() + geom_point(data = NULL, aes(y = resid, x = fit))

Residuals against predictors

mcmc = data.r2jags$BUGSoutput$sims.matrix[, c("beta0", "beta1")]
# generate a model matrix
newdata = data.frame(x = data$x)
Xmat = model.matrix(~x, newdata)
## get median parameter estimates
coefs = apply(mcmc, 2, median)
fit = as.vector(coefs %*% t(Xmat))
resid = data$y - fit
ggplot() + geom_point(data = NULL, aes(y = resid, x = data$x))

And now for studentized residuals

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

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

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

mcmc = data.r2jags$BUGSoutput$sims.matrix
# generate a model matrix
Xmat = model.matrix(~x, data)
## get median parameter estimates
coefs = mcmc[, c("beta0", "beta1")]
fit = coefs %*% t(Xmat)
## draw samples from this model
yRep = sapply(1:nrow(mcmc), function(i) rnorm(nrow(data), fit[i,
    ], mcmc[i, "sigma2"]))

Error in mcmc[i, "sigma2"]: subscript out of bounds

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)

Although residuals can be computed directly within Rstan, we can calculate them manually from the posteriors to be consistent across other approaches.

mcmc = as.matrix(data.rstan)[, c("beta0", "beta[1]")]
# generate a model matrix
newdata = data.frame(x = data$x)
Xmat = model.matrix(~x, newdata)
## get median parameter estimates
coefs = apply(mcmc, 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.matrix(data.rstan)[, c("beta0", "beta[1]")]
# generate a model matrix
newdata = data.frame(x = data$x)
Xmat = model.matrix(~x, newdata)
## get median parameter estimates
coefs = apply(mcmc, 2, median)
fit = as.vector(coefs %*% t(Xmat))
resid = data$y - fit
ggplot() + geom_point(data = NULL, aes(y = resid, x = data$x))

And now for studentized residuals

mcmc = as.matrix(data.rstan)[, c("beta0", "beta[1]")]
# generate a model matrix
newdata = data.frame(x = data$x)
Xmat = model.matrix(~x, newdata)
## get median parameter estimates
coefs = apply(mcmc, 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.rstan)
# generate a model matrix
Xmat = model.matrix(~x, data)
## get median parameter estimates
coefs = mcmc[, c("beta0", "beta[1]")]
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)

Residuals can be computed directly within RSTANARM.

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

Residuals against predictors

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

And now for studentized residuals

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

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

Lets 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)
newdata = data %>% cbind(t(y_pred)) %>% gather(key = "Rep", value = "Value",
    -y, -x)
ggplot(newdata, aes(Value, x = x)) + geom_violin(color = "blue",
    fill = "blue", alpha = 0.5) + geom_violin(data = data, aes(y = y,
    x = x), fill = "red", color = "red", alpha = 0.5)

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

## Calculate the fitted values
newdata = data.frame(x = seq(min(data$x, na.rm = TRUE), max(data$x,
    na.rm = TRUE), len = 1000))
fit = posterior_predict(data.rstanarm, newdata = newdata)
newdata = newdata %>% cbind(tidyMCMC(as.mcmc(fit), conf.int = TRUE,
    conf.method = "HPDinterval"))

ggplot(newdata, aes(y = estimate, x = x)) + geom_point(data = data,
    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()

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.

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)[, "Estimate"]
fit = fitted(data.brms)[, "Estimate"]
ggplot() + geom_point(data = NULL, aes(y = resid, x = fit))

Residuals against predictors

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

And now for studentized residuals

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

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

Lets 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)
newdata = data %>% cbind(t(y_pred)) %>% gather(key = "Rep", value = "Value",
    -y, -x)
ggplot(newdata, aes(Value, x = x)) + geom_violin(color = "blue",
    fill = "blue", alpha = 0.5) + geom_violin(data = data, aes(y = y,
    x = x), fill = "red", color = "red", alpha = 0.5)

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

## Calculate the fitted values
newdata = data.frame(x = seq(min(data$x, na.rm = TRUE), max(data$x,
    na.rm = TRUE), len = 1000))
fit = posterior_predict(data.brms, newdata = newdata)
newdata = newdata %>% cbind(tidyMCMC(as.mcmc(fit), conf.int = TRUE,
    conf.method = "HPDinterval"))

ggplot(newdata, aes(y = estimate, x = x)) + geom_point(data = data,
    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()

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.

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.

summary(data.mcmcpack)

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

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

              Mean      SD Naive SE Time-series SE
(Intercept) 40.788  2.9157 0.029157       0.028679
x           -1.538  0.3026 0.003026       0.003026
sigma2      30.453 14.0296 0.140296       0.162065

2. Quantiles for each variable:

              2.5%    25%    50%    75%   97.5%
(Intercept) 35.087 38.949 40.802 42.620 46.6754
x           -2.148 -1.729 -1.537 -1.349 -0.9395
sigma2      14.015 21.183 27.281 35.757 64.8614

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

         term  estimate  std.error  conf.low  conf.high
1 (Intercept) 40.787840  2.9156846 34.694040 46.2035365
2           x -1.538321  0.3026426 -2.133656 -0.9323802
3      sigma2 30.452600 14.0295572 10.955592 56.7338303

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.

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

}
## since values are less than zero
mcmcpvalue(data.mcmcpack[, 2])

[1] 2e-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.

print(data.r2jags)

Inference for Bugs model at "../downloads/BUGSscripts/tut7.2bS4.1.txt", 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
beta0     40.776   3.012 34.771 38.854 40.754  42.700  46.732 1.001  6500
beta1     -1.536   0.309 -2.155 -1.739 -1.534  -1.339  -0.917 1.001  8800
sigma      5.612   1.198  3.825  4.765  5.433   6.244   8.454 1.001 15000
deviance  98.913   2.916 95.586 96.795 98.131 100.243 106.547 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 = 4.3 and DIC = 103.2
DIC is an estimate of expected predictive error (lower deviance is better).

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

      term  estimate std.error  conf.low   conf.high
1    beta0 40.776236 3.0118914 34.644637  46.5544619
2    beta1 -1.535883 0.3094887 -2.117052  -0.8847727
3 deviance 98.913094 2.9160662 95.342003 104.6398121
4    sigma  5.611529 1.1977617  3.510032   7.8862896

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.

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

}
## since values are less than zero
mcmcpvalue(data.r2jags$BUGSoutput$sims.matrix[, c("beta1")])

[1] 0

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.

summary(data.rstan)

$summary
              mean     se_mean        sd       2.5%        25%        50%        75%       97.5%
beta[1]  -1.540077 0.005547721 0.2992531  -2.140083  -1.733530  -1.536049  -1.347254  -0.9578622
cbeta0   27.699678 0.026072380 1.4002154  24.952730  26.779818  27.686496  28.616701  30.4432020
sigma     5.378446 0.020515831 1.0852954   3.734084   4.614093   5.199031   5.958270   7.9476304
beta0    40.790336 0.053432423 2.9034858  35.075289  38.948498  40.737680  42.663585  46.4509997
lp__    -33.769878 0.026822294 1.3557193 -37.440064 -34.368897 -33.422015 -32.778343 -32.1970451
           n_eff      Rhat
beta[1] 2909.699 0.9994293
cbeta0  2884.220 1.0003349
sigma   2798.451 0.9995632
beta0   2952.770 0.9993640
lp__    2554.746 0.9995627

$c_summary
, , chains = chain:1

         stats
parameter       mean        sd       2.5%        25%        50%        75%       97.5%
  beta[1]  -1.531479 0.3068162  -2.109599  -1.735114  -1.539470  -1.334883  -0.9173666
  cbeta0   27.746909 1.4159460  24.954516  26.855717  27.747541  28.659675  30.4723348
  sigma     5.396365 1.0961031   3.734274   4.594885   5.232184   6.006155   7.8181447
  beta0    40.764477 3.0140699  34.540731  38.900198  40.880995  42.739393  46.5386786
  lp__    -33.809348 1.4068914 -37.550169 -34.393107 -33.462536 -32.782456 -32.1869446

, , chains = chain:2

         stats
parameter       mean        sd       2.5%        25%        50%        75%      97.5%
  beta[1]  -1.550123 0.2974626  -2.155754  -1.746184  -1.538468  -1.353810  -1.010301
  cbeta0   27.687988 1.3881674  25.134487  26.703722  27.683492  28.600342  30.383584
  sigma     5.402310 1.0729273   3.838954   4.640805   5.209596   5.950922   7.982861
  beta0    40.864033 2.8689516  35.483307  38.998427  40.801638  42.774958  46.193502
  lp__    -33.734080 1.3104041 -37.188086 -34.291847 -33.386710 -32.790364 -32.213962

, , chains = chain:3

         stats
parameter       mean        sd       2.5%        25%        50%        75%       97.5%
  beta[1]  -1.538631 0.2933243  -2.141996  -1.716454  -1.531518  -1.354423  -0.9678914
  cbeta0   27.664137 1.3964915  24.918672  26.754701  27.613705  28.591730  30.4340678
  sigma     5.336664 1.0866013   3.619856   4.595549   5.154616   5.908103   8.0011510
  beta0    40.742497 2.8255388  35.183039  38.892966  40.634539  42.437613  46.3971867
  lp__    -33.766205 1.3484235 -37.384673 -34.447354 -33.396535 -32.763341 -32.2065629

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

     term  estimate std.error  conf.low  conf.high      rhat  ess
1 beta[1] -1.540077 0.2992531 -2.136075 -0.9571987 0.9994293 2910
2  cbeta0 27.699678 1.4002154 25.186804 30.6173088 1.0003349 2884
3   sigma  5.378446 1.0852954  3.487439  7.4980289 0.9995632 2798
4   beta0 40.790336 2.9034858 34.878535 46.2187799 0.9993640 2953

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.

Also note that since our STAN model incorporated predictor centering, we have estimates of the intercept based on both centered (cbeta0) and uncentered data (beta0). Since the intercept from uncentered data is beyond the domain of our sampling data it has very little interpretability. However, the intercept based on centered data can be interpreted as the estimate of the response at the mean predictor (in this case 27.699678).

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

}
## since values are less than zero
mcmcpvalue(as.matrix(data.rstan)[, c("beta[1]")])

[1] 0

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.

summary(data.rstanarm)

Model Info:

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

Estimates:
                mean   sd    2.5%   25%   50%   75%   97.5%
(Intercept)    40.7    3.1  34.7   38.8  40.8  42.7  46.9  
x              -1.5    0.3  -2.2   -1.7  -1.5  -1.3  -0.9  
sigma           5.6    1.2   3.8    4.8   5.4   6.2   8.5  
mean_PPD       27.7    2.0  23.6   26.4  27.7  29.0  31.7  
log-posterior -65.2    1.3 -68.6  -65.8 -64.9 -64.2 -63.7  

Diagnostics:
              mcse Rhat n_eff
(Intercept)   0.1  1.0  2271 
x             0.0  1.0  2700 
sigma         0.0  1.0  1264 
mean_PPD      0.0  1.0  1808 
log-posterior 0.0  1.0  1078 

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

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

           term   estimate std.error   conf.low   conf.high      rhat  ess
1   (Intercept)  40.734015 3.0668570  34.412462  46.3425366 1.0000549 2271
2             x  -1.531248 0.3148917  -2.123538  -0.8981083 1.0001968 2700
3         sigma   5.619764 1.1740461   3.604894   7.9675118 0.9996872 1264
4      mean_PPD  27.698236 2.0410524  23.858059  31.9079814 1.0002406 1808
5 log-posterior -65.225515 1.3374104 -67.797185 -63.5678144 1.0007470 1078

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.

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

}
## since values are less than zero
mcmcpvalue(as.matrix(data.rstanarm)[, c("x")])

[1] 0

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, prob = 0.95)

                 2.5%      97.5%
(Intercept) 34.668050 46.8721257
x           -2.151347 -0.9122868
sigma        3.846747  8.4674463

An alternative way of quantifying the impact of a predictor is to compare models with and without the predictor. 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))

Computed from 2700 by 16 log-likelihood matrix

         Estimate  SE
elpd_loo    -51.2 3.2
p_loo         2.9 1.3
looic       102.4 6.3

Pareto k diagnostic values:
                         Count  Pct 
(-Inf, 0.5]   (good)     15    93.8%
 (0.5, 0.7]   (ok)        1     6.2%
   (0.7, 1]   (bad)       0     0.0%
   (1, Inf)   (very bad)  0     0.0%

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

(reduced = loo(update(data.rstanarm, formula = . ~ 1)))

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.00911 seconds (Warm-up)
               0.035428 seconds (Sampling)
               0.044538 seconds (Total)


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



 Elapsed Time: 0.008509 seconds (Warm-up)
               0.034266 seconds (Sampling)
               0.042775 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.011456 seconds (Warm-up)
               0.035874 seconds (Sampling)
               0.04733 seconds (Total)

Computed from 2700 by 16 log-likelihood matrix

         Estimate  SE
elpd_loo    -59.4 3.1
p_loo         2.1 1.0
looic       118.8 6.2

Pareto k diagnostic values:
                         Count  Pct 
(-Inf, 0.5]   (good)     15    93.8%
 (0.5, 0.7]   (ok)        0     0.0%
   (0.7, 1]   (bad)       1     6.2%
   (1, Inf)   (very bad)  0     0.0%
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 
     -8.2       1.8 

summary(data.brms)

 Family: gaussian(identity) 
Formula: y ~ x 
   Data: data (Number of observations: 16) 
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    40.85      3.07    34.90    47.24       2310    1
x            -1.55      0.32    -2.18    -0.90       2408    1

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     5.63      1.19      3.9     8.48       2091    1

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

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

         term  estimate std.error  conf.low  conf.high     rhat  ess
1 b_Intercept 40.853729 3.0714471 34.407211 46.6638946 1.001079 2310
2         b_x -1.545095 0.3205355 -2.141672 -0.8693169 1.001133 2408
3       sigma  5.631356 1.1888245  3.779939  8.2153548 1.001264 2091

Whilst 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 the two populations are identical (t=0).

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

}
## since values are less than zero
mcmcpvalue(as.matrix(data.brms)[, c("b_x")])

[1] 0.0007407407

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), prob = 0.95)

                  2.5%       97.5%
b_Intercept  34.899273  47.2408583
b_x          -2.179970  -0.8959258
sigma         3.895194   8.4787733
lp__        -36.892779 -31.5536581

An alternative way of quantifying the impact of a predictor is to compare models with and without the predictor. 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))

  LOOIC   SE
 102.43 6.22

(reduced = loo(update(data.brms, formula = . ~ 1)))

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

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


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

 Elapsed Time: 0.022215 seconds (Warm-up)
               0.029115 seconds (Sampling)
               0.05133 seconds (Total)


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

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


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

 Elapsed Time: 0.025044 seconds (Warm-up)
               0.021075 seconds (Sampling)
               0.046119 seconds (Total)


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

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


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

 Elapsed Time: 0.019155 seconds (Warm-up)
               0.022313 seconds (Sampling)
               0.041468 seconds (Total)


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

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


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

 Elapsed Time: 0.012441 seconds (Warm-up)
               0.018111 seconds (Sampling)
               0.030552 seconds (Total)


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

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


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

 Elapsed Time: 0.012443 seconds (Warm-up)
               0.016248 seconds (Sampling)
               0.028691 seconds (Total)


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

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 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1001 / 2000 [ 50%]  (Sampling)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 0.012169 seconds (Warm-up)
               0.015697 seconds (Sampling)
               0.027866 seconds (Total)

  LOOIC   SE
 118.26 6.03

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)

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.

mcmc = data.mcmcpack
## Calculate the fitted values
newdata = data.frame(x = seq(min(data$x, na.rm = TRUE), max(data$x, na.rm = TRUE),
    len = 1000))
Xmat = model.matrix(~x, newdata)
coefs = mcmc[, 1:2]
fit = coefs %*% t(Xmat)
newdata = newdata %>% cbind(tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval"))

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

If you wanted to represent sample data on the figure in such a simple example (single predictor) we could simply over- (or under-) lay the raw data.

ggplot(newdata, aes(y = estimate, x = x)) + geom_point(data = data, aes(y = y,
    x = x), color = "gray") + 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()

A more general solution would be to add the partial residuals to the figure. Partial residuals are the fitted values plus the residuals. In this simple case, that equates to exactly the same as the raw observations since $resid = obs - fitted$ and the fitted values depend only on the single predictor we are interested in.

## Calculate partial residuals fitted values
fdata = rdata = data
fMat = rMat = model.matrix(~x, fdata)
fit = as.vector(apply(coefs, 2, median) %*% t(fMat))
resid = as.vector(data$y - apply(coefs, 2, median) %*% t(rMat))
rdata = rdata %>% mutate(partial.resid = resid + fit)
ggplot(newdata, aes(y = estimate, x = 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") + scale_x_continuous("X") +
    theme_classic()

mcmc = data.r2jags$BUGSoutput$sims.matrix
## Calculate the fitted values
newdata = data.frame(x = seq(min(data$x, na.rm = TRUE), max(data$x, na.rm = TRUE),
    len = 1000))
Xmat = model.matrix(~x, newdata)
coefs = mcmc[, c("beta0", "beta1")]
fit = coefs %*% t(Xmat)
newdata = newdata %>% cbind(tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval"))
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()

If you wanted to represent sample data on the figure in such a simple example (single predictor) we could simply over- (or under-) lay the raw data.

ggplot(newdata, aes(y = estimate, x = x)) + geom_point(data = data, aes(y = y,
    x = x), color = "gray") + 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()

A more general solution would be to add the partial residuals to the figure. Partial residuals are the fitted values plus the residuals. In this simple case, that equates to exactly the same as the raw observations since $resid = obs - fitted $ and the fitted values depend only on the single predictor we are interested in.

## Calculate partial residuals fitted values
fdata = rdata = data
fMat = rMat = model.matrix(~x, fdata)
fit = as.vector(apply(coefs, 2, median) %*% t(fMat))
resid = as.vector(data$y - apply(coefs, 2, median) %*% t(rMat))
rdata = rdata %>% mutate(partial.resid = resid + fit)
ggplot(newdata, aes(y = estimate, x = 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") + scale_x_continuous("X") +
    theme_classic()

mcmc = as.matrix(data.rstan)
## Calculate the fitted values
newdata = data.frame(x = seq(min(data$x, na.rm = TRUE), max(data$x, na.rm = TRUE),
    len = 1000))
Xmat = model.matrix(~x, newdata)
coefs = mcmc[, c("beta0", "beta[1]")]
fit = coefs %*% t(Xmat)
newdata = newdata %>% cbind(tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval"))
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()

If you wanted to represent sample data on the figure in such a simple example (single predictor) we could simply over- (or under-) lay the raw data.

ggplot(newdata, aes(y = estimate, x = x)) + geom_point(data = data, aes(y = y,
    x = x), color = "gray") + 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()

A more general solution would be to add the partial residuals to the figure. Partial residuals are the fitted values plus the residuals. In this simple case, that equates to exactly the same as the raw observations since $resid = obs - fitted $ and the fitted values depend only on the single predictor we are interested in.

## Calculate partial residuals fitted values
fdata = rdata = data
fMat = rMat = model.matrix(~x, fdata)
fit = as.vector(apply(coefs, 2, median) %*% t(fMat))
resid = as.vector(data$y - apply(coefs, 2, median) %*% t(rMat))
rdata = rdata %>% mutate(partial.resid = resid + fit)
ggplot(newdata, aes(y = estimate, x = 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") + scale_x_continuous("X") +
    theme_classic()

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

## Calculate the fitted values
newdata = data.frame(x = seq(min(data$x, na.rm = TRUE), max(data$x,
    na.rm = TRUE), len = 1000))
# fit = posterior_predict(data.rstanarm, newdata=newdata)
fit = posterior_linpred(data.rstanarm, newdata = newdata)
newdata = newdata %>% cbind(tidyMCMC(as.mcmc(fit), conf.int = TRUE,
    conf.method = "HPDinterval"))

ggplot(newdata, aes(y = estimate, x = x)) + geom_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()

If you wanted to represent sample data on the figure in such a simple example (single predictor) we could simply over- (or under-) lay the raw data.

ggplot(newdata, aes(y = estimate, x = x)) + geom_point(data = data, aes(y = y,
    x = x), color = "gray") + 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()

A more general solution would be to add the partial residuals to the figure. Partial residuals are the fitted values plus the residuals. In this simple case, that equates to exactly the same as the raw observations since $resid = obs - fitted $ and the fitted values depend only on the single predictor we are interested in.

## Calculate partial residuals fitted values
fdata = rdata = data
pp = posterior_predict(data.rstanarm, newdata = fdata)
fit = as.vector(apply(pp, 2, median))
pp = posterior_predict(data.rstanarm, newdata = rdata)
resid = as.vector(data$y - apply(pp, 2, median))
rdata = rdata %>% mutate(partial.resid = resid + fit)
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") + scale_x_continuous("X") +
    theme_classic()

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

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

It is also possible to just use the marginal_effects function to generate the partial effects data which can then be plotted manually. Just note the names of the columns produced...

newdata = marginal_effects(data.brms)$x
newdata %>% head

         x  y cond__ estimate__     se__  lower__  upper__
1 1.000000 NA      1   39.33572 2.555412 33.98311 45.12585
2 1.151515 NA      1   39.10511 2.518356 33.85225 44.80221
3 1.303030 NA      1   38.86400 2.479195 33.69083 44.46928
4 1.454545 NA      1   38.62627 2.446745 33.54161 44.15459
5 1.606061 NA      1   38.39589 2.400776 33.37734 43.84009
6 1.757576 NA      1   38.17167 2.362441 33.21209 43.50193

ggplot(newdata, aes(y = estimate__, x = x)) + geom_line() + geom_ribbon(aes(ymin = lower__,
    ymax = upper__), fill = "blue", alpha = 0.3) + scale_y_continuous("Y") +
    scale_x_continuous("X") + theme_classic()

If you wanted to represent sample data on the figure in such a simple example (single predictor) we could simply over- (or under-) lay the raw data.

ggplot(newdata, aes(y = estimate__, x = x)) + geom_point(data = data, aes(y = y,
    x = x), color = "gray") + geom_line() + geom_ribbon(aes(ymin = lower__,
    ymax = upper__), fill = "blue", alpha = 0.3) + scale_y_continuous("Y") +
    scale_x_continuous("X") + theme_classic()

A more general solution would be to add the partial residuals to the figure. Partial residuals are the fitted values plus the residuals. In this simple case, that equates to exactly the same as the raw observations since $resid = obs - fitted $ and the fitted values depend only on the single predictor we are interested in.

mcmc = as.matrix(data.brms)
## Calculate the fitted values
newdata = data.frame(x = seq(min(data$x, na.rm = TRUE), max(data$x, na.rm = TRUE),
    len = 1000))
Xmat = model.matrix(~x, newdata)
coefs = mcmc[, c("b_Intercept", "b_x")]
fit = coefs %*% t(Xmat)
newdata = newdata %>% cbind(tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval"))

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

ggplot(newdata, aes(y = estimate, x = 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") + scale_x_continuous("X") +
    theme_classic()

In addition to deriving the distribution means for the slope parameter, we could make use of the Bayesian framework to derive the distribution of the effect size. In so doing, effect size could be considered as either the rate of change or alternatively, the difference between pairs of values along the predictor gradient. For the latter case, there are multiple ways of calculating an effect size, but the two most common are:

Raw effect size

the difference between two groups (as already calculated)

Cohen's D

the effect size standardized by division with the pooled standard deviation $$ D = (\mu_A - \mu_B)/\sigma $$

Percent effect size

expressing the effect size as a percent of one of the pairs. That is, whether you expressing a percentage increase or a percentage decline depends on which of the pairs of values are considered a reference value. Care must be exercised to ensure no division by zeros occur.

For simple linear models, effect size based on a rate is essentially the same as above except that it is expressed per unit of the predictor. Of course in many instances, one unit change in the predictor represents too subtle a shift in the underlying gradient to likely yield any ecologically meaningful or appreciable change in response.

Lets explore a range of effect sizes:

• Raw effect size between the largest and smallest x
• Cohen's D
• Percentage change between the largest and smallest x
• Fractional change between the largest and smallest x

Probability that a change in x is associated with greater than a 25% decline in y. Clearly, in order to explore this inference, we must first express the change in y as a percentage. This in turn requires us to calculate start and end points from which to calculate the magnitude of the effect (amount of decline in y) as well as the percentage decline. Hence, we start by predicting the distribution of y at the lowest and highest values of x.

mcmc = data.mcmcpack
newdata = data.frame(x = c(min(data$x, na.rm = TRUE), max(data$x, na.rm = TRUE)))
Xmat = model.matrix(~x, newdata)
coefs = mcmc[, c("(Intercept)", "x")]
fit = coefs %*% t(Xmat)
## Raw effect size
(RES = tidyMCMC(as.mcmc(fit[, 2] - fit[, 1]), conf.int = TRUE, conf.method = "HPDinterval"))

  term  estimate std.error  conf.low conf.high
1 var1 -23.07482   4.53964 -32.00484  -13.9857

## Cohen's D
cohenD = (fit[, 2] - fit[, 1])/sqrt(mcmc[, "sigma2"])
(cohenDES = tidyMCMC(as.mcmc(cohenD), conf.int = TRUE, conf.method = "HPDinterval"))

  term estimate std.error  conf.low conf.high
1 var1 -4.43846  1.168116 -6.752846  -2.21918

# Percentage change (relative to Group A)
ESp = 100 * (fit[, 2] - fit[, 1])/fit[, 1]
(PES = tidyMCMC(as.mcmc(ESp), conf.int = TRUE, conf.method = "HPDinterval"))

  term  estimate std.error conf.low conf.high
1 var1 -58.38874  8.549874 -75.0835  -41.3122

# Probability that the effect is greater than 25% (a decline of >25%)
sum(-1 * ESp > 25)/length(ESp)

[1] 0.9986

## fractional change
fit = fit[fit[, 2] > 0, ]
(FES = tidyMCMC(as.mcmc(fit[, 2]/fit[, 1]), conf.int = TRUE, conf.method = "HPDinterval"))

  term  estimate  std.error conf.low conf.high
1 var1 0.4161556 0.08539488  0.24904 0.5865297

Conclusions:

• On average, Y declines by -23.074816 over the observed range of x. We are 95% confident that the decline is between -32.0048416 and -13.9857023.
• The Cohen's D associated with a change over the observed range of x is -4.4384597.
• On average, Y declines by -58.389% over the observed range of x. We are 95% confident that the decline is between -75.084% and -41.312%.
• The probability that Y declines by more than 25% over the observed range of x is 0.999.
• On average, Y declines by a factor of 0.416% over the observed range of x. We are 95% confident that the decline is between a factor of 0.249% and 0.587%.

Lets explore a range of effect sizes:

• Raw effect size between the largest and smallest x
• Cohen's D
• Percentage change between the largest and smallest x
• Fractional change between the largest and smallest x

Probability that a change in x is associated with greater than a 25% decline in y. Clearly, in order to explore this inference, we must first express the change in y as a percentage. This in turn requires us to calculate start and end points from which to calculate the magnitude of the effect (amount of decline in y) as well as the percentage decline. Hence, we start by predicting the distribution of y at the lowest and highest values of x.

mcmc = data.r2jags$BUGSoutput$sims.matrix
newdata = data.frame(x = c(min(data$x, na.rm = TRUE), max(data$x, na.rm = TRUE)))
Xmat = model.matrix(~x, newdata)
coefs = mcmc[, c("beta0", "beta1")]
fit = coefs %*% t(Xmat)
## Raw effect size
(RES = tidyMCMC(as.mcmc(fit[, 2] - fit[, 1]), conf.int = TRUE, conf.method = "HPDinterval"))

  term  estimate std.error  conf.low conf.high
1 var1 -23.03824   4.64233 -31.75579 -13.27159

## Cohen's D
cohenD = (fit[, 2] - fit[, 1])/mcmc[, "sigma"]
(cohenDES = tidyMCMC(as.mcmc(cohenD), conf.int = TRUE, conf.method = "HPDinterval"))

  term  estimate std.error  conf.low conf.high
1 var1 -4.276647  1.166878 -6.610581 -2.039906

# Percentage change (relative to Group A)
ESp = 100 * (fit[, 2] - fit[, 1])/fit[, 1]
(PES = tidyMCMC(as.mcmc(ESp), conf.int = TRUE, conf.method = "HPDinterval"))

  term  estimate std.error  conf.low conf.high
1 var1 -58.28882  8.755287 -75.37096 -40.51335

# Probability that the effect is greater than 25% (a decline of >25%)
sum(-1 * ESp > 25)/length(ESp)

[1] 0.9982

## fractional change
fit = fit[fit[, 2] > 0, ]
(FES = tidyMCMC(as.mcmc(fit[, 2]/fit[, 1]), conf.int = TRUE, conf.method = "HPDinterval"))

  term  estimate  std.error  conf.low conf.high
1 var1 0.4171118 0.08755287 0.2462904 0.5948665

Conclusions:

• On average, Y declines by -23.0382448 over the observed range of x. We are 95% confident that the decline is between -31.7557852 and -13.2715912.
• The Cohen's D associated with a change over the observed range of x is -4.2766474.
• On average, Y declines by -58.289% over the observed range of x. We are 95% confident that the decline is between -75.371% and -40.513%.
• The probability that Y declines by more than 25% over the observed range of x is 0.998.
• On average, Y declines by a factor of 0.417% over the observed range of x. We are 95% confident that the decline is between a factor of 0.246% and 0.595%.

Lets explore a range of effect sizes:

• Raw effect size between the largest and smallest x
• Cohen's D
• Percentage change between the largest and smallest x
• Fractional change between the largest and smallest x

Probability that a change in x is associated with greater than a 25% decline in y. Clearly, in order to explore this inference, we must first express the change in y as a percentage. This in turn requires us to calculate start and end points from which to calculate the magnitude of the effect (amount of decline in y) as well as the percentage decline. Hence, we start by predicting the distribution of y at the lowest and highest values of x.

mcmc = as.matrix(data.rstan)
newdata = data.frame(x = c(min(data$x, na.rm = TRUE), max(data$x, na.rm = TRUE)))
Xmat = model.matrix(~x, newdata)
coefs = mcmc[, c("beta0", "beta[1]")]
fit = coefs %*% t(Xmat)
## Raw effect size
(RES = tidyMCMC(as.mcmc(fit[, 2] - fit[, 1]), conf.int = TRUE, conf.method = "HPDinterval"))

  term  estimate std.error  conf.low conf.high
1 var1 -23.10116  4.488796 -32.04112 -14.35798

## Cohen's D
cohenD = (fit[, 2] - fit[, 1])/mcmc[, "sigma"]
(cohenDES = tidyMCMC(as.mcmc(cohenD), conf.int = TRUE, conf.method = "HPDinterval"))

  term  estimate std.error  conf.low conf.high
1 var1 -4.453289  1.150329 -6.766335 -2.289957

# Percentage change (relative to Group A)
ESp = 100 * (fit[, 2] - fit[, 1])/fit[, 1]
(PES = tidyMCMC(as.mcmc(ESp), conf.int = TRUE, conf.method = "HPDinterval"))

  term  estimate std.error  conf.low conf.high
1 var1 -58.46377  8.445858 -75.83271 -42.96449

# Probability that the effect is greater than 25% (a decline of >25%)
sum(-1 * ESp > 25)/length(ESp)

[1] 0.9993333

## fractional change
fit = fit[fit[, 2] > 0, ]
(FES = tidyMCMC(as.mcmc(fit[, 2]/fit[, 1]), conf.int = TRUE, conf.method = "HPDinterval"))

  term  estimate  std.error  conf.low conf.high
1 var1 0.4153623 0.08445858 0.2416729 0.5703551

Conclusions:

• On average, Y declines by -23.101161 over the observed range of x. We are 95% confident that the decline is between -32.0411213 and -14.3579804.
• The Cohen's D associated with a change over the observed range of x is -4.4532885.
• On average, Y declines by -58.464% over the observed range of x. We are 95% confident that the decline is between -75.833% and -42.964%.
• The probability that Y declines by more than 25% over the observed range of x is 0.999.
• On average, Y declines by a factor of 0.415% over the observed range of x. We are 95% confident that the decline is between a factor of 0.242% and 0.570%.

Lets explore a range of effect sizes:

• Raw effect size between the largest and smallest x
• Cohen's D
• Percentage change between the largest and smallest x
• Fractional change between the largest and smallest x

Probability that a change in x is associated with greater than a 25% decline in y. Clearly, in order to explore this inference, we must first express the change in y as a percentage. This in turn requires us to calculate start and end points from which to calculate the magnitude of the effect (amount of decline in y) as well as the percentage decline. Hence, we start by predicting the distribution of y at the lowest and highest values of x.

mcmc = as.matrix(data.rstanarm)
newdata = data.frame(x = c(min(data$x, na.rm = TRUE), max(data$x, na.rm = TRUE)))
Xmat = model.matrix(~x, newdata)
coefs = mcmc[, c("(Intercept)", "x")]
fit = coefs %*% t(Xmat)
## Raw effect size
(RES = tidyMCMC(as.mcmc(fit[, 2] - fit[, 1]), conf.int = TRUE, conf.method = "HPDinterval"))

  term  estimate std.error  conf.low conf.high
1 var1 -22.96872  4.723375 -31.85307 -13.47162

## Cohen's D
cohenD = (fit[, 2] - fit[, 1])/mcmc[, "sigma"]
(cohenDES = tidyMCMC(as.mcmc(cohenD), conf.int = TRUE, conf.method = "HPDinterval"))

  term  estimate std.error  conf.low conf.high
1 var1 -4.252559   1.16655 -6.533503 -2.049772

# Percentage change (relative to Group A)
ESp = 100 * (fit[, 2] - fit[, 1])/fit[, 1]
(PES = tidyMCMC(as.mcmc(ESp), conf.int = TRUE, conf.method = "HPDinterval"))

  term  estimate std.error  conf.low conf.high
1 var1 -58.14203  8.945745 -74.70909 -39.84752

# Probability that the effect is greater than 25% (a decline of >25%)
sum(-1 * ESp > 25)/length(ESp)

[1] 0.9985185

## fractional change
fit = fit[fit[, 2] > 0, ]
(FES = tidyMCMC(as.mcmc(fit[, 2]/fit[, 1]), conf.int = TRUE, conf.method = "HPDinterval"))

  term  estimate  std.error  conf.low conf.high
1 var1 0.4185797 0.08945745 0.2529091 0.6015248

Conclusions:

• On average, Y declines by -22.9687244 over the observed range of x. We are 95% confident that the decline is between -31.8530656 and -13.4716247.
• The Cohen's D associated with a change over the observed range of x is -4.2525591.
• On average, Y declines by -58.142% over the observed range of x. We are 95% confident that the decline is between -74.709% and -39.848%.
• The probability that Y declines by more than 25% over the observed range of x is 0.999.
• On average, Y declines by a factor of 0.419% over the observed range of x. We are 95% confident that the decline is between a factor of 0.253% and 0.602%.

Lets explore a range of effect sizes:

• Raw effect size between the largest and smallest x
• Cohen's D
• Percentage change between the largest and smallest x
• Fractional change between the largest and smallest x

Probability that a change in x is associated with greater than a 25% decline in y. Clearly, in order to explore this inference, we must first express the change in y as a percentage. This in turn requires us to calculate start and end points from which to calculate the magnitude of the effect (amount of decline in y) as well as the percentage decline. Hence, we start by predicting the distribution of y at the lowest and highest values of x.

mcmc = as.matrix(data.brms)
newdata = data.frame(x = c(min(data$x, na.rm = TRUE), max(data$x, na.rm = TRUE)))
Xmat = model.matrix(~x, newdata)
coefs = mcmc[, c("b_Intercept", "b_x")]
fit = coefs %*% t(Xmat)
## Raw effect size
(RES = tidyMCMC(as.mcmc(fit[, 2] - fit[, 1]), conf.int = TRUE, conf.method = "HPDinterval"))

  term  estimate std.error  conf.low conf.high
1 var1 -23.17642  4.808033 -32.12509 -13.03975

## Cohen's D
cohenD = (fit[, 2] - fit[, 1])/mcmc[, "sigma"]
(cohenDES = tidyMCMC(as.mcmc(cohenD), conf.int = TRUE, conf.method = "HPDinterval"))

  term estimate std.error  conf.low conf.high
1 var1 -4.28766  1.176568 -6.596148 -1.951322

# Percentage change (relative to Group A)
ESp = 100 * (fit[, 2] - fit[, 1])/fit[, 1]
(PES = tidyMCMC(as.mcmc(ESp), conf.int = TRUE, conf.method = "HPDinterval"))

  term  estimate std.error  conf.low conf.high
1 var1 -58.50167  9.213223 -75.36007 -39.18994

# Probability that the effect is greater than 25% (a decline of >25%)
sum(-1 * ESp > 25)/length(ESp)

[1] 0.9966667

## fractional change
fit = fit[fit[, 2] > 0, ]
(FES = tidyMCMC(as.mcmc(fit[, 2]/fit[, 1]), conf.int = TRUE, conf.method = "HPDinterval"))

  term  estimate  std.error  conf.low conf.high
1 var1 0.4149833 0.09213223 0.2463993 0.6081006

Conclusions:

• On average, Y declines by -23.1764222 over the observed range of x. We are 95% confident that the decline is between -32.125086 and -13.0397529.
• The Cohen's D associated with a change over the observed range of x is -4.2876597.
• On average, Y declines by -58.502% over the observed range of x. We are 95% confident that the decline is between -75.360% and -39.190%.
• The probability that Y declines by more than 25% over the observed range of x is 0.997.
• On average, Y declines by a factor of 0.415% over the observed range of x. We are 95% confident that the decline is between a factor of 0.246% and 0.608%.

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

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

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

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

library(broom)
mcmc = data.mcmcpack
sd.x = abs(mcmc[, "x"]) * sd(data$x)
# generate a model matrix
newdata = data.frame(x = data$x)
Xmat = model.matrix(~x, newdata)
## get median parameter estimates
coefs = mcmc[, 1:2]
fit = coefs %*% t(Xmat)
resid = sweep(fit, 2, data$y, "-")
sd.resid = apply(resid, 1, sd)
sd.all = cbind(sd.x, sd.resid)
(fpsd = tidyMCMC(sd.all, conf.int = TRUE, conf.method = "HPDinterval"))

      term estimate std.error conf.low conf.high
1     sd.x 7.324011 1.4401645 4.439017 10.158235
2 sd.resid 5.113917 0.3125066 4.916653  5.689996

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

      term estimate std.error conf.low conf.high
1     sd.x 59.81066  5.410107 48.00330  64.16565
2 sd.resid 40.18934  5.410107 35.83435  51.99670

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

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

library(broom)
mcmc = data.r2jags$BUGSoutput$sims.matrix
sd.x = abs(mcmc[, "beta1"]) * sd(data$x)
# generate a model matrix
newdata = data.frame(x = data$x)
Xmat = model.matrix(~x, newdata)
## get median parameter estimates
coefs = mcmc[, 1:2]
fit = coefs %*% t(Xmat)
resid = sweep(fit, 2, data$y, "-")
sd.resid = apply(resid, 1, sd)
sd.all = cbind(sd.x, sd.resid)
(fpsd = tidyMCMC(sd.all, conf.int = TRUE, conf.method = "HPDinterval"))

      term estimate std.error conf.low conf.high
1     sd.x 7.312266 1.4734608 4.212361 10.079185
2 sd.resid 5.123532 0.3064685 4.916653  5.728441

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

      term estimate std.error conf.low conf.high
1     sd.x 59.76241  5.634456 47.19196  64.16565
2 sd.resid 40.23759  5.634456 35.83435  52.80804

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

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

library(broom)
mcmc = as.matrix(data.rstan)
sd.x = abs(mcmc[, "beta[1]"]) * sd(data$x)
# generate a model matrix
newdata = data.frame(x = data$x)
Xmat = model.matrix(~x, newdata)
## get median parameter estimates
coefs = mcmc[, c("beta0", "beta[1]")]
fit = coefs %*% t(Xmat)
resid = sweep(fit, 2, data$y, "-")
sd.resid = apply(resid, 1, sd)
sd.all = cbind(sd.x, sd.resid)
(fpsd = tidyMCMC(sd.all, conf.int = TRUE, conf.method = "HPDinterval"))

      term estimate std.error conf.low conf.high
1     sd.x 7.332235 1.4247297 4.557177 10.169750
2 sd.resid 5.110913 0.2865312 4.916653  5.666491

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

      term estimate std.error conf.low conf.high
1     sd.x 59.79741  5.318334 48.03619  64.16563
2 sd.resid 40.20259  5.318334 35.83437  51.96381

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

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

library(broom)
mcmc = as.matrix(data.rstanarm)
sd.x = abs(mcmc[, "x"]) * sd(data$x)
# generate a model matrix
newdata = data.frame(x = data$x)
Xmat = model.matrix(~x, newdata)
## get median parameter estimates
coefs = mcmc[, 1:2]
fit = coefs %*% t(Xmat)
resid = sweep(fit, 2, data$y, "-")
sd.resid = apply(resid, 1, sd)
sd.all = cbind(sd.x, sd.resid)
(fpsd = tidyMCMC(sd.all, conf.int = TRUE, conf.method = "HPDinterval"))

      term estimate std.error conf.low conf.high
1     sd.x 7.290200 1.4991843 4.275851 10.110062
2 sd.resid 5.130423 0.3148904 4.916653  5.738437

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

      term estimate std.error conf.low conf.high
1     sd.x 59.72643  5.864832 46.65895  64.16561
2 sd.resid 40.27357  5.864832 35.83439  53.34105

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

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

library(broom)
mcmc = as.matrix(data.brms)
sd.x = abs(mcmc[, "b_x"]) * sd(data$x)
# generate a model matrix
newdata = data.frame(x = data$x)
Xmat = model.matrix(~x, newdata)
## get median parameter estimates
coefs = mcmc[, 1:2]
fit = coefs %*% t(Xmat)
resid = sweep(fit, 2, data$y, "-")
sd.resid = apply(resid, 1, sd)
sd.all = cbind(sd.x, sd.resid)
(fpsd = tidyMCMC(sd.all, conf.int = TRUE, conf.method = "HPDinterval"))

      term estimate std.error conf.low conf.high
1     sd.x 7.357990 1.5170233 4.138776 10.196400
2 sd.resid 5.135858 0.3567607 4.916653  5.786968

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

      term estimate std.error conf.low conf.high
1     sd.x 59.88833   5.89935 46.88813  64.16558
2 sd.resid 40.11167   5.89935 35.83442  53.11187

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

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

In a frequentist context, the $R^2$ value is seen as a useful indicator of goodness of fit. Whilst it has long been acknowledged that this measure is not appropriate for comparing models (for such purposes information criterion such as AIC are more appropriate), it is nevertheless useful for estimating the amount (percent) of variance explained by the model.

In a frequentist context, $R^2$ is calculated as the variance in predicted values divided by the variance in the observed (response) values. Unfortunately, this classical formulation does not translate simply into a Bayesian context since the equivalently calculated numerator can be larger than the an equivalently calculated denominator - thereby resulting in an $R^2$ greater than 100%. Gelman, Goodrich, Gabry, and Ali (2017) proposed an alternative formulation in which the denominator comprises the sum of the explained variance and the variance of the residuals.

So in the standard regression model notation of: $$ \begin{align} y_i \sim{}& N(\mu_i, \sigma)\\ \mu_i =& \mathbf{X}\boldsymbol{\beta} \end{align} $$ The $R^2$ could be formulated as: $$ R^2 = \frac{\sigma^2_f}{\sigma^2_f + \sigma^2_e} $$ where $\sigma^2_f = var(\mu)$, ($\mu = \mathbf{X}\boldsymbol{\beta})$) and for Gaussian models $\sigma^2_e = var(y-\mu)$