ANOVA example
library(brms)
library(bayesplot)Let’s look at a repeated measures ANOVA. I’ll use a built-in dataset showing differences in weights between chicks on different diets. This is repeated measures data so we’ll be using a multilevel model.Technically these are repeated measures data across time, but I’m going to ignore time for this example (not recommended - you’d want to examine time as a predictor and diet effect moderator).
summary(ChickWeight) # again, I hate numbered conditions!## weight Time Chick Diet
## Min. : 35.0 Min. : 0.00 13 : 12 1:220
## 1st Qu.: 63.0 1st Qu.: 4.00 9 : 12 2:120
## Median :103.0 Median :10.00 20 : 12 3:120
## Mean :121.8 Mean :10.72 10 : 12 4:118
## 3rd Qu.:163.8 3rd Qu.:16.00 17 : 12
## Max. :373.0 Max. :21.00 19 : 12
## (Other):506plot(weight~Diet, data=ChickWeight)
Other than the intercept, the other regression weights correspond to condition differences. Those differences certainly can’t be more than 100, so we’ll use N(0,50) for a default prior.
mycw <- brm(weight~Diet+(1|Chick), data=ChickWeight, prior=c(set_prior("normal(0,50)", class= "b")), cores=4)## Compiling the C++ model## Start samplingsummary(mycw)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: weight ~ Diet + (1 | Chick)
## Data: ChickWeight (Number of observations: 578)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~Chick (Number of levels: 50)
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept) 17.47 4.67 7.94 26.60 1268 1.00
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept 102.84 6.01 91.00 114.55 4000 1.00
## Diet2 19.27 9.95 -0.43 38.63 3443 1.00
## Diet3 38.85 9.96 19.25 58.06 4000 1.00
## Diet4 31.28 10.47 9.92 51.91 4000 1.00
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 67.47 2.10 63.51 71.63 4000 1.00
##
## 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).contrasts(ChickWeight$Diet) # Shows the coding of the categorical variable## 2 3 4
## 1 0 0 0
## 2 1 0 0
## 3 0 1 0
## 4 0 0 1plot(mycw)
Based on background knowledge, I was aware that weights would tend to have a long upper tail. Perhaps we should assume that weight is normally distributed and we should be using a Gamma distribution. We’ll start by looking at the residuals, and then do a Bayesian Gamma regression and compare the results to that from the original Gaussian model.
hist(resid(mycw)) # Check residual distribution
So, I’m going to test an alternative model assuming gamma distributions (long-tailed) weights. I’ll use a log-link function. Moving into log world means that my regression weights will be smaller. If I leave the range at normal(0,50), I’m essentially using flat priors in log space. So, the prior range is now N(0,1).
mycwGamma <- brm(weight~Diet+(1|Chick), data=ChickWeight, prior=c(set_prior("normal(0,1)", class= "b")), family=Gamma(link="log"), cores=4)## Compiling the C++ model## Start samplingIn my experience, some generalized models throw errors due to randomized starting values. You can try adding “inits=0” to the brm if they do. Not done here, but it’s a tip in case you encounter it.
Let’s check the two models using WAIC and a posterior parameter check.
waic(mycw, mycwGamma) # gamma is much better## WAIC SE
## mycw 6534.85 33.77
## mycwGamma 6343.25 32.66
## mycw - mycwGamma 191.61 17.70pp_check(mycw) # posterior parameter check. Checks if density plot for observed (y) is similar to that predicted (yrep)## Using 10 posterior samples for ppc type 'dens_overlay' by default.
pp_check(mycwGamma)## Using 10 posterior samples for ppc type 'dens_overlay' by default.
Yep, good evidence that the Gamma-based model provides a much better fit than the Gaussian-based model.
Hypothesis Testing
Testing for specific comparisons by adding parameter values. I’m doing two sets of comparisons - one in log-transformed space and one backtransformed to the original scale. Because the link function was “log”, the “exp” command is the opposite function and thus used for backtransforming.
hypGamma = hypothesis(
mycwGamma,
hypothesis = c(
# Look only at Diet1 which is coded as Intercept in our dummy coding
"Intercept = 0",
"exp(Intercept) = 0",
# compare two diets, Diet 1 vs. Diet4 then Diet2 vs. Diet4
# Note for Diet 1 vs. Diet 4, the difference is the Diet4 dummy variable
"Diet4 = 0",
# But in backtransformed space, I have to add intercept back in
"exp(Intercept + Diet4) - exp(Intercept) = 0",
# Now for Diet2 vs. 4. Intercept is common to both so I don't need it in the log-based comparison...
"Diet2 - Diet4 = 0",
# but I have to add it back in when backtransforming.
"exp(Intercept + Diet2) - exp(Intercept + Diet4) = 0"
))
hypGamma## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Intercept) = 0 4.60 0.06 4.49 4.70 NA
## 2 (exp(Intercept)) = 0 99.46 5.52 88.91 110.34 NA
## 3 (Diet4) = 0 0.31 0.09 0.12 0.50 NA
## 4 (exp(Intercept+Di... = 0 35.78 11.73 13.53 60.20 NA
## 5 (Diet2-Diet4) = 0 -0.11 0.11 -0.34 0.10 NA
## 6 (exp(Intercept+Di... = 0 -14.64 13.65 -42.74 12.29 NA
## Post.Prob Star
## 1 NA *
## 2 NA *
## 3 NA *
## 4 NA *
## 5 NA
## 6 NA
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.plot(hypGamma)
mcmc_dens(hypGamma$samples, pars=c("H2", "H4", "H6"))
mcmc_areas(hypGamma$samples, pars=c("H2", "H4", "H6"), prob=.95)
mcmc_intervals(hypGamma$samples, pars=c("H2", "H4", "H6"))
Now, let’s just plot the four groups using this method. Note that marginal_effects works, too, so we’ll start there with the 95% credible interval and then use our bayesplot functions to show it in other ways.
marginal_effects(mycwGamma)
hypGamma = hypothesis(
mycwGamma,
hypothesis = c(
"exp(Intercept) = 0",
"exp(Intercept + Diet2) = 0",
"exp(Intercept + Diet3) = 0",
"exp(Intercept + Diet4) = 0"
))
toplot = hypGamma$samples # extracts the distributions for the hypotheses
# change labels to the columns which defaulted to "H1", "H2", "H3", "H4" for the four hypotheses
names(toplot) = c("Diet 1", "Diet 2", "Diet 3", "Diet 4")
mcmc_dens(toplot)
mcmc_areas(toplot, prob=.95)
mcmc_intervals(toplot, prob_outer=.95, prob=.68)