Last week, we learned the basics of the linear model: how we might try to predict a continuous outcome based on a single continuous or categorical predictor. This week, we’re going to tackle the much more common situation of having multiple (even many) possible predictors for a single continuous outcome.
(This week we’re not going to explore either multiple or categorical outcomes very much. The former requires us to talk about multivariate methods, and the latter requires us to talk about the general linear model, which we will talk about later in the class.)
We’re going to switch up the order slightly and start by talking about multi-way ANOVA, which extends the basic model of ANOVA you have become somewhat familiar with to the case of multiple, categorical predictor variables. As you might guess, we are then going to talk about multiple regression, which extends the linear model to predicting a single, continuous outcome from multiple continuous and/or categorical predictors.
Data example
This week we’re going to use a new data set to demonstrate the approaches we’re covering. This data comes from Amy Moore, who did her MS with Dr. Amanda Stewart here in FST. Amy graduated in 2020.
This dataset was presented and published in Amy’s thesis (Moore et al. 2020). It is data on hydrogen sulfide production from yeast used in cider making. Briefly, apple juice–the substrate for cider production–is often poor in some of the nutrients that yeast need. This results in the yeast entering stressed metabolic states, which in turn causes production of hydrogen sulfide. For those of you who are unaware, hydrogen sulfide smells rotten and eggy. It is not a good smell.
Therefore, cider producers are interested in developing production strategies that minimize hydrogen sulfide production. Yeast manufacturers have started to develop strains of yeast that are supposed to not produce as much hydrogen sulfide, and are also working on producing yeast nutrient products that should help reduce hydrogen sulfide. To produce this data set, Amy investigated several strategies to determine if they successfully reduced hydrogen sulfide production in experimental ciders.
There are three experimental factors that Amy investigated:
Yeast strain: 3 levels (different strains of yeast).
Nutrient: 5 levels (different nutrients including a Control no-nutrient treatment).
Timing: 2 levels (whether nutrients were added all at once or in 2 additions separated by a time period).
This gives us a total of \(3\times5\times2=30\) possible experimental combinations. Amy also conducted 3 replicate fermentations for each treatment, giving us a total of \(90\) observations.
# A tibble: 30 × 7
# Groups: Strain, Nutrient [15]
Strain Nutrient Timing mean_hs se_hs ll ul
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 M2 Experimental Split 717. 72.6 571. 862.
2 M2 DAP Split 622. 66.2 489. 754.
3 M2 DAP Single 583. 44.1 495. 672.
4 M2 K Split 550 76.4 397. 703.
5 M2 O Split 505 5 495 515
6 M2 K Single 500 14.4 471. 529.
7 M2 O Single 475 14.4 446. 504.
8 M2 Control Split 433. 33.3 367. 500
9 M2 Control Single 430 30 370 490
10 M2 Experimental Single 430 76.5 277. 583.
# ℹ 20 more rows
With this much data, looking at numerical summaries is a little hard. We could experiment with using table visualization tools like gt or flextable to make this easier to examine, but I think a visual exploration is going to be easier.
h2s_data|>ggplot(aes(y =HS))+geom_point(aes(x =Nutrient, color =Timing, shape =Strain), size =2)+theme_bw()+labs(y =expression(H[2]*S), color ="Nutrient addition", shape ="Yeast strain")+scale_color_brewer(palette ="Dark2")+scale_shape_manual(values =c(1, 4, 8))
We can see that there are definite patterns in the data. The M2 strain of yeast, for example, seems to produce more hydrogen sulfide across all treatments–this kind of makes sense, if we think about it: from the data we have available, the yeast strain is going to be the biggest contributor to its metabolic behavior. The other variables are a little harder to parse visually: for example, we see that the Split addition treatment seems to give both the highest and lowest values of hydrogen sulfide production.
Multi-way ANOVA
From last week, we might remember how to test whether yeast strain (a categorical variable) is related to hydrogen sulfide production (we can treat this as a continuous outcome, although of course negative values are not possible).
# How do we ask R to run the linear model for yeast strain?
But if we run this model, of course, we are implicitly representing all values of hydrogen sulfide production coming from a single strain with a single mean, ignoring (literally averaging over) the different values that might come from the different Timing or Nutrient treatments:
# Here's what we're doingh2s_data|>group_by(Strain)|>summarize(mean_yeast_hs =mean(HS))
# A tibble: 2 × 2
Timing mean_timing_hs
<chr> <dbl>
1 Single 267.
2 Split 305.
Ideally, we’d like to not run 3 different ANOVAs. Not only does this ignore information, making our models inaccurate, it also increases the chance of us finding a chance result–remember that if we stick to traditional frequentist statistics, we are using a Type I error chance to decide whether our result it improbable enough to reject the null. In English, that means that every individual experiment we have a small chance of rejecting the null when we shouldn’t, and slicing the data in multiple ways is a really good way to do this. Instead, we want to be able to find some model to describe a figure kind of like this one:
h2s_data|>group_by(Strain, Nutrient, Timing)|>summarize(mean_hs =mean(HS))|>ggplot(aes(y =mean_hs))+geom_line(aes(x =Nutrient, color =Timing, linetype =Strain, group =interaction(Timing, Strain)))+geom_point(aes(x =Nutrient, color =Timing))+theme_bw()+labs(y =expression(mean~H[2]*S), color ="Nutrient addition", linetype ="Yeast strain")+scale_color_brewer(palette ="Dark2")
Here, we can see the “cell means” (the intersection of Strain, Nutrient, and Timing), and suddenly we see a clearer picture of possible patterns…
As you have probably guessed, a linear model with multiple predictors is the answer to this dilemma. We can fit a model which will simultaneously take into account the effect of multiple predictors, and even their interactions (within the limits of our data). When we have only categorical predictors, this is typically called “multi-way ANOVA”.
The simplest version of this model only models the additive effect of the multiple treatments. We don’t allow the effect of, say, Timing to be different for the different Nutrient types.
h2s_aov<-h2s_data|>aov(HS~Nutrient+Timing+Strain, data =_)summary(h2s_aov)
We use the same function, aov(), to build multi-way ANOVA models, and the only difference is that we will add a new set of tools to our formula interface:
+ adds a main effect for a term to our formula
<variable 1> * <variable 2> interacts the two (or more) variables as well as adding all lower-level interactions and main effects
<variable 1> : <variable 2> adds only the interaction term between the two (or more) variables
In the above model we see that we have “large” effects for Timing and Strain, but no real evidence of an effect of nutrient.
An aside: this week we’re not going to use as many resampling or permutation methods. Without getting into the gory details, it is not trivial to resample experimental designs with multiple factors in a way that respects the original experimental design. For the purposes of this class we will stick with using the formula-based (least-squares) analyses.
A better alternative is Bayesian methododology–we’ll talk a little bit about this towards the end of the class. I have wanted to start teaching this class using Bayesian approaches since I began to teach, but it requires me to teach a bunch of background that a) I think will cause everyone to shut down and b) I am not really qualified to do! So I decided the most useful approach for the class as a whole would be for me to present practical methods that are standard for most of our disciplines.
But, here in the aside, I can say clearly that I think that many of these methods are not the best, and that learning to use basic Bayesian approaches–which is pretty easy in R using packages like rstanarm, brms, and tidybayes–is really a good idea!
The main-effects model is not very realistic: from our plotted data we see that there seem to be possible different effects of one variable at different levels of another.
Before we decide that we know what is important for hydrogen sulfide production, though, we should investigate possible interactions among our main effects. From the plot above, we can see that there might be different effects of, for example, Timing for different Strain values.
h2s_data|>aov(HS~Nutrient*Strain*Timing, data =_)|>summary()
Interestingly, there appears to be a very large effect for the interactionNutrient:Strain, even thoughpossibly causing the lack of consistent main effect for Nutrient. This tells us, essentially, that the effect of different nutrients is not consistent across different yeast strains. We can modify the plot we made above to visualize this interaction.
h2s_data|># We are "marginalizing" across Timing for this visualizationgroup_by(Strain, Nutrient)|>summarize(mean_hs =mean(HS))|>ggplot(aes(y =mean_hs))+geom_point(aes(x =Nutrient, color =Strain))+geom_line(aes(x =Nutrient, color =Strain, group =Strain, linetype =Strain))+theme_bw()+labs(y =expression(mean~H[2]*S), Color ="Yeast strain")+scale_color_brewer(palette ="Dark2")
We can see clearly that the M2 yeast not only produces more hydrogen sulfide overall; its hydrogen sulfide levels are relatively higher for nutrients that the other two yeasts produce less hydrogen sulfide for, and vice versa. This is a good example of an interaction effect, in that it is a clear pattern of reversed relative effect of one factor (Nutrient) across a second (yeast Strain).
Post-hoc tests
A natural question we have after this kind of model is to evaluate the specific marginal effects of different treatment levels: is the cell mean of level A of Factor 1 different from the cell mean of level B or C or D or…? Last week, we showed how to use both formula-based and computational approaches to estimate these differences; this week, with multiple factors, we’re going to stick to the formulas. A package that works well with the lm() / aov() functions is the agricolae package, which provides various posthoc tests–we looked at this last week. Let’s see what happens when we apply the agricolae::HSD.test() function to an ANOVA model with multiple predictors.
Strain emmean SE df lower.CL upper.CL
EC1118 117 16.4 82 84 149
ICV OKAY 217 16.4 82 185 250
M2 524 16.4 82 492 557
Results are averaged over the levels of: Nutrient, Timing
Confidence level used: 0.95
h2s_strain_margins|>cld(Letters =letters)
Strain emmean SE df lower.CL upper.CL .group
EC1118 117 16.4 82 84 149 a
ICV OKAY 217 16.4 82 185 250 b
M2 524 16.4 82 492 557 c
Results are averaged over the levels of: Nutrient, Timing
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 3 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
# There are some alternative ways to present these resultsplot(h2s_strain_margins, comparisons =TRUE)
EC1118 ICV OKAY M2
EC1118 [117] 0.0001 <.0001
ICV OKAY -101 [217] <.0001
M2 -408 -307 [525]
Row and column labels: Strain
Upper triangle: P values adjust = "tukey"
Diagonal: [Estimates] (emmean)
Lower triangle: Comparisons (estimate) earlier vs. later
We have to specify the trt argument to marginalize across the other factors. This means that we are estimating main-effect post-hoc tests. The multcomp::cld() function is very useful for producing easy-to-export tabular data for publications, because it will generate the useful .group labels for the different means in a way that can be output to Excel and then to Word.
How would we estimate the post-hoc tests for the Nutrient factor?
An advantage of the emmeans package is that it provides robust support for estimating interaction marginal means.
h2s_interaction_margins<-h2s_aov|>emmeans(specs =~Nutrient:Strain)# But interactions can create difficult to interpret CLDsh2s_interaction_margins|>cld()
Nutrient Strain emmean SE df lower.CL upper.CL
DAP EC1118 103 25 82 52.9 152
Control EC1118 112 25 82 62.6 162
Experimental EC1118 113 25 82 63.5 163
O EC1118 122 25 82 72.2 172
K EC1118 133 25 82 83.0 182
DAP ICV OKAY 203 25 82 153.7 253
Control ICV OKAY 213 25 82 163.4 263
Experimental ICV OKAY 214 25 82 164.3 264
O ICV OKAY 223 25 82 173.0 272
K ICV OKAY 233 25 82 183.8 283
DAP M2 511 25 82 460.9 560
Control M2 520 25 82 470.6 570
Experimental M2 521 25 82 471.5 571
O M2 530 25 82 480.2 580
K M2 541 25 82 491.0 590
.group
1234567890ABCDEF
12 4 6 8 0 B D GHIJKLMN
1 345 89 BC GH J L OPQR
123 7890A G IJK OP ST
1234567 E GHI M O Q S U
GHIJKLMNOPQRSTUV
3 5 7 9 A C EF OPQRSTUV
2 67 0A DEF I K MN STUV
456 BCDEF H LMN QR UV
890ABCD F JKL N P R T V
W
W
W
W
W
Results are averaged over the levels of: Timing
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 15 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
ANOVA is a special case of the general linear model in which all predictors are categorical–they take on only specified values, which are treated as “dummy variables” when we fit the model. If we have continuous predictors, we’re back in the land of the original linear model, which is typically called “multiple regression” (same issues as before).
Just as simple regression and ANOVA are analogues of each other, multiple regression is the equivalent of multi-way ANOVA when one or more of the predictors is continuous instead of categorical (when we have a mix of continuous and categorical variables we tend to deal with them through either multiple regression or with something called ANCOVA: Analysis of Co-Variance).
Multiple regression involves predicting the value of an outcome from many continuous or categorical predictors. Multiple regression is also a topic that deserves (and is taught as) a course-length topic in its own right (see for example STAT 6634). I will be barely scratching the surface of multiple regression here–as always, I am going to attempt to show you how to implement a method in R.
We can frame our objective in multiple regression as predicting the value of a single outcome \(y\) from a set of predictors \(x_1, x_2, ..., x_n\).
To identify variables \(x_1, x_2, ...\) etc that are actually important in predicting \(y\)
The first task is usually evaluated in the same way as we did in simple regression–by investigating goodness-of-fit statistics like \(R^2\). This is often interpreted as overall quality of the model.
The second task is often considered more important in research applications: we want to understand which \(x\)-variables significantly predict \(y\). We usually assess this by investigating statistical significance of the \(\beta\)-coefficients for each predictor.
Finally, a third and important task in multiple regression is identifying interacting predictor variables. In my usual, imprecise interpretation, an interactive effect can be stated as “the effect of \(x_1\) is different at different values of \(x_2\)”–this means that we have a multiplicative effect in our model.
The model given above omits interactive terms, which would look like:
You’ll notice that here Amy averaged across the 3 fermentation reps for each cider sample, which we had previously used to estimate the ANOVA model for hydrogen sulfide. That loss of variability would be a problem if we tried to estimate the same model here: we wouldn’t be able to estimate the 3-way interaction (if we really cared about that). But here we’re instead going to see if we can model hydrogen sulfide production from level of ethanol and titratable acidity. Let’s look at the correlations between all of our variables:
We can see that, while TA and pH are highly correlated, they still seem to have some independent information, so we’ll include them all in a model for Total_H2S. This probably isn’t a very good model for two reasons:
All of these are fermentation-related outcomes–none of them are really better \(x\) variables than another.
We have no theory-based reason for assuming that any of these \(x_n\) variables would be useful for explaining our \(y\), hydrogen-sulfide production.
Nevertheless, nothing about the linear model forbids us from using variables that have no explicit relationship to each other in a single model. Remember the sharks and ice cream?
cider_lm_interaction_model<-cider_chem_data|>lm(Total_H2S~Ethanol*TA*pH, data =_)cider_lm_interaction_model|>summary()
Call:
lm(formula = Total_H2S ~ Ethanol * TA * pH, data = cider_chem_data)
Residuals:
Min 1Q Median 3Q Max
-278.44 -92.79 28.97 105.43 287.19
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3949016 8586408 0.460 0.650
Ethanol -626468 1385070 -0.452 0.655
TA -691911 1351300 -0.512 0.614
pH -1090351 2350495 -0.464 0.647
Ethanol:TA 109747 218025 0.503 0.620
Ethanol:pH 173002 379139 0.456 0.653
TA:pH 191582 369966 0.518 0.610
Ethanol:TA:pH -30392 59689 -0.509 0.616
Residual standard error: 172.6 on 22 degrees of freedom
Multiple R-squared: 0.3622, Adjusted R-squared: 0.1593
F-statistic: 1.785 on 7 and 22 DF, p-value: 0.1412
Unsurprisingly we find none of our effects are significant in this model–but we are also trying to model 8 effects (including a 3-way interaction!) with only 30 data points. Let’s see if a simpler model might make more sense.
cider_lm_additive_model<-cider_chem_data|>lm(Total_H2S~Ethanol+TA+pH, data =_)cider_lm_additive_model|>summary()
Call:
lm(formula = Total_H2S ~ Ethanol + TA + pH, data = cider_chem_data)
Residuals:
Min 1Q Median 3Q Max
-253.94 -138.45 18.38 109.36 336.41
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6993.5 5206.4 1.343 0.19080
Ethanol -1791.6 585.2 -3.061 0.00507 **
TA -105.8 110.8 -0.955 0.34834
pH 1391.7 1112.4 1.251 0.22204
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 168.1 on 26 degrees of freedom
Multiple R-squared: 0.2855, Adjusted R-squared: 0.203
F-statistic: 3.463 on 3 and 26 DF, p-value: 0.03064
Hey, now we have a significant effect! Interesting.
This reveals one of the main weaknesses of a data-driven approach to inference. We can just keep making choices–sometimes called researcher degrees of freedom until we find a “significant” effect, and then we can say “Hey, we found something!” But remember that I started out this section by saying that this model is kind of silly?
We have just managed to write a model that makes a correlation–perhaps a chance one, perhaps a real one–seem causal. Semantically, by putting ethanol in the \(x\) position, we’re saying that “a change in ethanol of 1 unit (here, in Brix) is associated with a change of -1792 units of hydrogen sulfide (in mg/L)”. But we can’t really change ethanol in the fermentation–it’s another outcome of the fermentation metabolism that is producing both ethanol and hydrogen sulfide. What’s really happening is that yeast might produce one or the other, and we’re seeing that metabolic process outcome measured as a combination of ethanol and hydrogen sulfide.
When presented with multiple models, there are many fit criteria that have been proposed for model selection.
While not technically a fit criterion, the \(R^2\) is often informally used this way: recall that \(R^2\) is the “amount of variance” in \(y\) that is explained by the model. Higher \(R^2\) means the model is better, right? Well, a basic problem is that \(R^2\) will always increase as you add more parameters (predictors) to the linear model. We can see that in the above two models.
A very common example of a more rigorous criterion is Akaike’s Information Criterion (AIC), which takes into account a combination of the observed fit and the number of parameters–it is lowest for a model that “fits the best” with the fewest parameters. Thus, for two models that “fit” the same, the AIC will be lower (better) for the simpler model.
Call:
lm(formula = Total_H2S ~ Ethanol + TA + pH, data = cider_chem_data)
Residuals:
Min 1Q Median 3Q Max
-253.94 -138.45 18.38 109.36 336.41
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6993.5 5206.4 1.343 0.19080
Ethanol -1791.6 585.2 -3.061 0.00507 **
TA -105.8 110.8 -0.955 0.34834
pH 1391.7 1112.4 1.251 0.22204
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 168.1 on 26 degrees of freedom
Multiple R-squared: 0.2855, Adjusted R-squared: 0.203
F-statistic: 3.463 on 3 and 26 DF, p-value: 0.03064
# Our more complex modelsummary(cider_lm_interaction_model)
Call:
lm(formula = Total_H2S ~ Ethanol * TA * pH, data = cider_chem_data)
Residuals:
Min 1Q Median 3Q Max
-278.44 -92.79 28.97 105.43 287.19
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3949016 8586408 0.460 0.650
Ethanol -626468 1385070 -0.452 0.655
TA -691911 1351300 -0.512 0.614
pH -1090351 2350495 -0.464 0.647
Ethanol:TA 109747 218025 0.503 0.620
Ethanol:pH 173002 379139 0.456 0.653
TA:pH 191582 369966 0.518 0.610
Ethanol:TA:pH -30392 59689 -0.509 0.616
Residual standard error: 172.6 on 22 degrees of freedom
Multiple R-squared: 0.3622, Adjusted R-squared: 0.1593
F-statistic: 1.785 on 7 and 22 DF, p-value: 0.1412
cat("AIC for simple model:", AIC(cider_lm_additive_model), "\n")
AIC for simple model: 398.308
cat("AIC for complex model:", AIC(cider_lm_interaction_model))
AIC for complex model: 402.9006
So while \(R^2\) is higher for the complex model, AIC is lower. Note the “adjusted \(R^2\)” also captures some of this complexity.
So can we just use these criteria to select the best model? Let’s wait a minute to answer that question…
Linear models with both continuous and categorical predictors
Nothing about the linear model requires us to use only continuous (“regression) or only categorical (”ANOVA”) predictors. We can mix the two using the same set of tools. Let’s take a look at one last facet of Amy Moore’s data.
If we looked at the Excel file, we’d notice that Amy had also recorded to total fermentation duration that each cider took. Perhaps duration is related to hydrogen sulfide production: maybe if fermentations are “stuck” (taking longer), the yeast become stressed and switch metabolic pathways, producing more hydrogen sulfide.
Getting these data will require a little wrangling–we’ll need to join together two datasets. Unfortunately, Amy didn’t give unique IDs to her samples, so we will have to assume that the rows, which each refer to a fermentation replicate, are in the same order for all the data tables.
Honestly… no! But perhaps there is if we account for the different yeast, nutrient, and nutrient addition timing. We can do this using the same lm() tool we’ve been using.
h2s_anova_with_hrs<-h2s_data|>lm(HS~(Strain+Nutrient+Timing)*hrs, data =_)h2s_anova_with_hrs|>summary()
What does this model hypothesize? Note that I used a combination of + and * in the model.
Unsurprisingly, we see that hrs remains a poor predictor of hydrogen sulfide production. But we do see that there seems to be some kind of interaction between hrs and some levels of the Strain factor, and some levels of the Nutrient factor. I think we’re really fishing here (e.g., I do not think these are “real” important effects), but let’s examine one way we might look into these kinds of effects.
This is a good example of a visualization showing an interaction between a continuous and categorical predictor. We see the lines of best fit for hrs (linear models by least squares) for each level of Strain. It does actually look like there might be a different pattern for yeast M2 than for the other two yeasts. I don’t have a good explanation for this–which is why I said that, as nice as this picture is, I think we might be picking up on noise in the data. Is this pattern really there, or are we just trying all possible models, and exploiting random error in finding something that “looks good” (visually or numerically)?
Model selection
We’ve been talking about terms like “model selection”, which are jargon-y, so let’s refocus on what we’re doing: we are answering the question of “what of our experimental variables actually matter in explaining or predicting an outcome we care about?” In our example, we are asking, “if we manipulate aspects of the fermentation process, can we reduce hydrogen sulfide production?”
I’ve shown off some methods, like AIC, that seem to give us an objective answer on which model is “better”. Why don’t we just use these and make the process “objective”? While we will be getting into this question in detail when we talk about critical data literacy, here are several things to consider:
Just remember the following: we are making predictions from data, not explaining causality.
Without a theory that explains why any relationship is one we would expect to see, we run into the problem of just letting our coincidence-finding systems make up a story for us, which leads us to…
Without expertise, we are just exploiting patterns in the data that are quite possibly either accidental or are caused by a third, underlying variable. Remember the shark-attack example? I don’t know much about yeast, so I have no idea if theory would predict the kind of relationship with strain, fermentation time, and hydrogen sulfide production we just “found”.
The reason we cannot use “objective” processes is that, unfortunately, objectivity doesn’t exist! At least, not in the way we wish it did. For example, if the AIC isn’t real. This isn’t some universal constant that exists–it is a relationship derived from the Maximum Likelihood Estimate of a model with the number of parameters in the model. Essentially, this expresses a human hypothesis of what makes a model function well. But there are other criteria we could minimize that instantiate different arguments. Just because something is a number doesn’t mean it is objective.
We’ll get into this a bit more when we get into critical data literacy/thinking with data in the next weeks.
One of the motivating factors for me to stop thinking in terms of finding the “best” model fit using automatic/“objective” methods has been an increased attention to Bayesian methods. You may know Bayes Theory:
Without forcing you to go on a journey with me that’s taken me some years to get my head around, the basic point here is that what we typically use in our everyday statistics is just part of the righthand equation: \(P(Data|Hypothesis)\), where the hypothesis in question is our typicall null hypothesis–something boring and unlikely to be true. We never really think that the effect is actually 0, etc. We really care about the lefthand side: we want to know how probable our hypothesis is given our prior estimate of its probability (\(P(Hypothesis)\)) and the actual likelihood of the data we gathered given our hypothesis (\(P(Data|Hypothesis)\)).
The main point is that this approach requires us to state interesting hypotheses to start, and give a rough estimate of how probable we think those hypotheses are. As scientists, we can do that. For example, I think that it is pretty likely that choice of yeast, choice of nutrient, and how we feed the yeast have a moderately large effect on the amount of hydrogen sulfide that will be produced. I think this is true because I have heard a number of cidermakers express this general opinion (even if they don’t know the exact relationship) and it is the topic of ongoing research, in which others have found these relationships.
We can then use this prior knowledge to make a more informed ANOVA model than the one we started the class with.
# We're going to use the `rstanarm` package because it provides a simple method# for an anova-like model without the wrangling the (otherwise often nicer)# `brms` package requireslibrary(rstanarm)library(tidybayes)options(mc.cores =parallel::detectCores())h2s_bayes_aov<-stan_aov(formula =HS~Nutrient*Strain*Timing, data =h2s_data,# We need to give our prior estimate of how much# variance our model explains. 50% (0.5) seems# conservative based on the analysis we did above. prior =R2(location =0.5),# adapt_delta helps avoid some fitting problems adapt_delta =0.999)h2s_bayes_aov
Without going into too many details, a Bayesian approach simulates plausible values of our parameters of interest–the effect of our treatment levels and their interactions–in a way that respects the hypothesized prior distributions and the data we observed. The modern Bayesian approaches use a simulation method called Markov Chain Monte Carlo (MCMC; we’re not going into much detail here). In my mind–and this is a naive statement, not a statement of fact by a statistician–Bayesian approaches (which are implemented based on a method called ) can fill in the gaps where we cannot apply naive simulation like bootstrapping and permutation (and, really, surpass even the basic things we can do, at the cost of some intellectual and computational complexity).
A key advantage to Bayesian methods is that Bayes’ Theorem, regardless of the approach, estimates the probability of a hypothesis (operationalized as parameter values) given a prior hypothesis (plausible values of those parameters) and observed data. This is the opposite of what we typically do with significance testing, in which we estimate the probability of data given a particular (e.g., null hypothesis). With Bayesian methods, we end up getting a plausible range of the parameters we care about.
In a Bayesian approach, just like in our bootstrapping or permutation approach, we will make many “draws” from a distribution of possible results. The way we get to them is different (and I am not going to be covering it), but the end result is similar. We can then marginalize or split our estimated probabilities in order to do posthoc testing that gives us nice estimates of uncertainty. We see above that Timing seems to be important for hydrogen sulfide production, but its difficult to actually estimate the average effect of Timing (all at once vs split) because Timing is interacted with our other variables.
To get at this effect, we can average over all other effects by examining the generated simulation.
simulated_Split<-h2s_bayes_aov|># Here our "new data" is an experimental design with 90 rows with our other 2# factors as usual (Nutrient and Strain), but with all Timing set to be# "Split". The function will simulate outcomes for these combinations based# on what was calculated from our original model.posterior_predict(newdata =h2s_data|>mutate(Timing ="Split"))|>as_tibble()|>pivot_longer(everything(), values_to ="Split")simulated_Single<-h2s_bayes_aov|># Here our "new data" is just a version of our experimental factors where all# of the Timing is set to "Single".posterior_predict(newdata =h2s_data|>mutate(Timing ="Single"))|>as_tibble()|>pivot_longer(everything(), values_to ="Single")bind_cols(simulated_Split[, 2], simulated_Single[, 2])|>mutate(`Split - Single` =Split-Single)|>pivot_longer(everything())|>mutate(name =factor(name)|>fct_inorder())|>ggplot(aes(x =value))+geom_histogram(color ="black", fill ="grey", bins =30)+stat_slabinterval(point_interval =mode_hdi, width =0.95, color ="red")+facet_wrap(~name)+theme_classic()+labs(x =expression(estimated~H[2]*S~production), y =NULL)
We can see that the effect of nutrient addition (Timing) is really quite small. In addition, we can see the marked multimodality we get from the combination of our Strain variables in the data–the hump you see in the Split and Single estimates comes from the M2 yeast, which was so different. This implies that the Timing main effect probably isn’t important–we see that it spans 0 and is as likely to be negative as it is positive.
On the other hand, we are pretty sure that we should see some large main effects for our yeast Strain, even in the presence of such significant interactions. But we don’t just have a binary/two-level factor here–Strain has 3 levels. To compare these, we need to do a bit of basic algebra. We can of course directly compare 2/3 levels, but if we want to compare all 3 at once, we need to remember to average them appropriately to keep the scale controlled.
# When we have multiple levels we want to compare, the direct wrangling can be a# little lengthy. There are some shortcuts in `tidybayes` that I haven't# explored as much.simulated_EC1118<-h2s_bayes_aov|>posterior_predict(newdata =h2s_data|>mutate(Strain ="EC1118"))|>as_tibble()|>pivot_longer(everything(), values_to ="EC1118")simulated_ICVOK<-h2s_bayes_aov|>posterior_predict(newdata =h2s_data|>mutate(Strain ="ICV OKAY"))|>as_tibble()|>pivot_longer(everything(), values_to ="ICV OKAY")simulated_M2<-h2s_bayes_aov|>posterior_predict(newdata =h2s_data|>mutate(Strain ="M2"))|>as_tibble()|>pivot_longer(everything(), values_to ="M2")bind_cols(simulated_EC1118[, 2], simulated_ICVOK[, 2], simulated_M2[, 2])|># If we want to do contrasts between >2 treatments, we need to remember to# average appropriately.mutate(`M2 vs. other yeasts` =M2-(EC1118+`ICV OKAY`)/2, `ICV OKAY vs EC1118` =`ICV OKAY`-EC1118)|>pivot_longer(everything())|>mutate(name =factor(name)|>fct_inorder())|>ggplot(aes(x =value))+geom_vline(xintercept =0, color ="black", linetype ="dashed")+geom_histogram(color ="black", fill ="grey", bins =30)+stat_slabinterval(point_interval =mode_hdi, width =0.95, color ="red")+facet_wrap(~name)+theme_classic()+labs(x =expression(estimated~H[2]*S~production), y =NULL)
Here we can see that indeed we have what look like rather large effects of yeast Strain. The M2 yeast definitely produces more hydrogen sulfide, on average, than the other 2 yeast strains: we see a large, positive estimate for the contrast. On the other hand, we are less certain that ICV OKAY and EC1118 are actually that different. While the modal estimate for the effect is positive, a range of estimates that span positive and negative numbers are compatible with the data–we should feel less confident saying that there is a large, consistent difference in hydrogen sulfide production between ECV OKAY and EC1118. I didn’t bother to include the contrast between M2 and EC1118: it is clear given the shape of the data that it will be large and positive.
If we look back at our estimated effects from the complete model, we see that the 3-way interaction between strain, nutrient, and addition time seems to be quite large. We can use the same simulation-from-the-posterior approach to estimate this effect.
Returning to our example, we really might want to look at some derived parameter–perhaps the marginal difference between estimated parameters of our model for two levels of one of our factors.
# Let's say we want to examine the estimated predicted difference between two# ciders, made with the same method (M2 yeast, DAP nutrient) except for split# nutrient additionsimulated_results<-h2s_bayes_aov|>posterior_epred(newdata =tibble(Strain =c("M2", "M2"), Nutrient =c("Experimental", "Experimental"), Timing =c("Single", "Split")))|># Now with a little wrangling we can calculate the estimated marginal# difference in H2S productionas_tibble()|>rename(`Single @ M2, Experimental` =1, `Split @ M2, Experimental` =2)|>mutate(`Split - Single @ M2, Experimental` =`Split @ M2, Experimental`-`Single @ M2, Experimental`)simulated_results
simulated_results|>pivot_longer(everything())|>mutate(name =factor(name)|>fct_inorder())|>ggplot(aes(x =value))+geom_histogram(bins =30, fill ="grey", color ="black")+stat_slabinterval(point_interval =mode_hdi, width =0.95, color ="red")+facet_wrap(~name, nrow =3)+theme_classic()+labs(x =expression(mean~H[2]*S~production), y =NULL)
Here we can see that, indeed, the estimated effect of the nutrient addition timing for the M2 yeast and the Experimental nutrient is large and positive. This is probably driving the estimated positive effect of the Timing main effect–but it really is limited to this yeast. With a Bayesian MCMC simulation approach, we can examine these different possible effects carefully, without resoriting to indirect post-hoc testing.
Recommended reading for Bayesian methods
Bayesian methods are a bit of lift–you might wonder how I determined to use the type of posterior simulation I used above, for example. I recommend the following resources (perhaps in reverse order, in fact).
Moore, Amy N, Jacob Lahne, Elizabeth Burzynski-Chang, et al. 2020. “Impact of Yeast Nutrient Supplementation on Hydrogen Sulfide Production During Cider Fermentation.”World Brewing Congress.
