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.
> 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"
> 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"
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"
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.