Chapter 23 Bayesian Analysis of Linear Contrasts

Linear contrasts are extremely important in the analysis of experimental data. The linear.contrast() function in bcfcdae can also handle Bayesian fits from bglmm(). Please read Standard Analysis of Linear Contrasts for the background on the contrasts we will be demonstrating.

All of the functions to do contrasts work on a linear model with one or more factor variables as predictors. So let’s set up a couple of model fits using the rat liver weight data and the resin lifetime data.
> library(cfcdae)
> library(bcfcdae)
> library(rstan)
> data(RatLiverWeight)
> bfit.RLW <- bglmm(weight~diet,data=RatLiverWeight,
+                   quiet=TRUE,seed=10024,adapt_delta = .99)

Note that we needed to increase adapt_delta to avoid divergent transitions.

We need to check our diagnostics; we don’t show the results here, but everything looks fine:

plot(bfit.RLW,"trace")
plot(bfit.RLW,"autocor")
summary(bfit.RLW)[,9:10]
The usage of linear.contrast() is the same whether or not the model is Bayesian, but the output is different for Bayesian fits.
> linear.contrast(bfit.RLW,diet,c(1/3,1/3,1/3,-1))
  post. mean   post. sd      lower       upper approx BF
1 -0.2245658 0.08922923 -0.3958321 -0.04152788  11.04806
attr(,"class")
[1] "linear.contrast"

The output is now the posterior mean and standard deviation for the contrast, a Bayesian credible interval for the contrast, and an approximate Bayes factor for the full model relative to the same model with the model coefficients constrained so that the contrast is zero. Here, the interval does not contain zero and the Bayes factor prefers the full model.

You can change the coverage for the interval in the same way we did for non-Bayesian inference:
> linear.contrast(bfit.RLW,diet,c(1/3,1/3,1/3,-1),conf=.99)
  post. mean   post. sd    lower       upper approx BF
1 -0.2245658 0.08922923 -0.45495 0.003107889  11.04806
attr(,"class")
[1] "linear.contrast"
You can also request an exact Bayes factor for the contrast; this is slower because it must fit another Bayesian model.
> linear.contrast(bfit.RLW,diet,c(1/3,1/3,1/3,-1),exactBF = TRUE)
Iteration: 1
Iteration: 2
Iteration: 3
Iteration: 4
Warning: There were 2 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Iteration: 1
Iteration: 2
Iteration: 3
Iteration: 4
Iteration: 5
  post. mean   post. sd      lower       upper       BF
1 -0.2245658 0.08922923 -0.3958321 -0.04152788 7.621167
attr(,"class")
[1] "linear.contrast"
Divergent steps in the MCMC are bad. One potential solution is to increase adapt_delta (must be less than 1).
> linear.contrast(bfit.RLW,diet,c(1/3,1/3,1/3,-1),exactBF = TRUE,adapt_delta=.999)
Iteration: 1
Iteration: 2
Iteration: 3
Iteration: 4
Iteration: 5
Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Iteration: 1
Iteration: 2
Iteration: 3
Iteration: 4
Iteration: 5
Iteration: 6
  post. mean   post. sd      lower       upper       BF
1 -0.2245658 0.08922923 -0.3958321 -0.04152788 7.694299
attr(,"class")
[1] "linear.contrast"
Just to emphasize what linear.contrast() is doing when you ask for an exact Bayes factor, it is refitting the Bayesian model with the additional argument add.constraints=list(diet=c(1/3,1/3,1/3,-1)). In this model, you will always have the average of the first three diet effects equal to the last diet effect. To see that, look at
> bfit2.RLW <- bglmm(weight~diet,data=RatLiverWeight,quiet=TRUE,
+                   add.constraints=list(diet=c(1/3,1/3,1/3,-1)),
+                   adapt_delta=.999,show.redundant.effects = TRUE,
+                   seed=10025)
Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
> summary(bfit2.RLW)[,1]
    (Intercept)           diet1           diet2           diet3           diet4          sigma0 
       3.72e+00        4.14e-02       -2.95e-02       -1.19e-02        3.01e-18        2.45e-01 
sigma.Intercept      sigma.diet 
       4.35e+00        8.54e-02 

The four diet effects must always add to zero. If we impose the constraint that the sum of the first three must equal the last one, then the last one and the sum of the first three must also be zero. We see that in the coefficients.

As usual, we can do more than one contrast in the command:
> linear.contrast(bfit.RLW,diet,cbind(c(1/3,1/3,1/3,-1),c(1,-.5,-.5,0)))
  post. mean   post. sd       lower       upper approx BF
1 -0.2245658 0.08922923 -0.39583210 -0.04152788 11.048056
2  0.1246026 0.08986251 -0.04880934  0.30312266  1.217267
attr(,"class")
[1] "linear.contrast"

Here we also compared the first mean to the average of the second and third means. There is little to choose between constraining this contrast to be zero or letting it vary freely.