10 - The linear model
Introduction
Last week we got a taste of statistical inference, with a focus on using computational approaches to estimate stability/confidence for our parameters (bootstrapping) and for testing whether null hypotheses are plausible (permutation). This week we’re going to review (because you learned these in statistics, right?) the concept of the linear model–very broadly–which encompasses both simple and multiple regression as well as one-way and multi-way ANOVA.
Today we’re going to start by talking about the general problem of predicting one continuous (interval or ratio) variable from another. This is the classical recipe for what is usually called linear regression–we’ll talk about the lm() function in R, how to interpret results, and look into a bootstrapping approach for estimating stability of our estimates and a permutation approach that will help us see whether our results are significant. We will also talk a tiny bit about the history of the method and why I will be trying to call it the linear model instead of “regression”.
We will spend some time looking at how to plot results from this linear model using a combination of geom_point()/geom_jitter() and geom_smooth(), as well as faceting. We’ll spend relatively little time on the assumptions we should really examine, but I’ll try to mention them.
Then we will spend a little while talking about a particularly vexing kind of variable type you may already have encountered in R: factor-type data. These data are key in both the analysis and the plotting of ANOVA.
Then we’ll look at a very closely related problem: how to predict a continuous (interval or ratio) variable from a discrete (interval or ordinal) predictor. This is the classic recipe for ANOVA (analysis of variance), and we will look at the aov() function (a thin wrapper for the lm()) as well as pointing towards some other, more powerful packages (including the afex package). We’ll also revisit our permutation approach for ANOVA. And, again, we’ll note that this is really more properly just called linear modeling, and why we might want to move away from “ANOVA”.
We will also spend some time on plotting results from this kind of linear model, with a special focus on interaction (mean-effect) plots. These plots give powerful visualizations of both the main effects (treatment effects) and interaction effects (differential effects of one treatment based on levels of another treatment).
We will not be looking in any detail at using multiple predictors (categorical or continuous) in the linear model this week; that material will be introduced next week, when we will think more broadly about different methods for exploring how complex, interacting independent variables can be investigated in our models.
Datasets for today
Berry data (continuous predictor)
The first set of data we’ll be working with today is from research I conducted with Driscoll’s, Inc. in 2019. It is a subset of data from a Central Location Test on berry liking, in which we asked subjects to rate their liking on different sets of berries (blackberries, blueberries, strawberries, raspberries) using different kinds of rating scales. Some of the results of the study were published in a paper by Yeung et al. (2021). If you’re really interested in the details you can find them there, as well as some examples of applied bootstrapping and permutation tests!
Today, however, we’re just interested in using the data to explore some applications of the linear model. Let’s take a look at our data:
| Name | berry_data |
| Number of rows | 3350 |
| Number of columns | 8 |
| _______________________ | |
| Column type frequency: | |
| character | 3 |
| numeric | 5 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| age | 270 | 0.92 | 5 | 5 | 0 | 3 | 0 |
| gender | 270 | 0.92 | 4 | 6 | 0 | 3 | 0 |
| berry | 0 | 1.00 | 9 | 10 | 0 | 4 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| subject | 0 | 1.00 | 1075965.50 | 136631.18 | 1031946 | 1032642 | 1033289 | 1034014 | 1522576 | ▇▁▁▁▁ |
| sample | 0 | 1.00 | 3.40 | 1.66 | 1 | 2 | 3 | 5 | 6 | ▇▃▃▃▃ |
| 9pt_overall | 905 | 0.73 | 5.68 | 2.10 | 1 | 4 | 6 | 7 | 9 | ▂▅▂▇▅ |
| us_overall | 790 | 0.76 | 8.49 | 3.81 | 0 | 6 | 9 | 11 | 15 | ▃▅▇▇▅ |
| lms_overall | 848 | 0.75 | 20.24 | 43.74 | -100 | -9 | 26 | 56 | 100 | ▁▃▆▇▃ |
This data is semi-tidy: each row has 3 observations on it: the rating of the berry by the particular subject in terms of the 9-pt hedonic scale (9pt_overall), the Labeled Magnitude Scale (lms_overall) and the Unstructured Line Scale (us_overall). The first of these 3 scales produces ordinal data that is often treated as interval-level, and the latter two are reputed to produce interval-level data (again, see Yeung et al. (2021) for details on some of these assumptions and their validity in this dataset).
You’ll notice we have some missing data in our dataset. We can take a closer look at this using the naniar package (which is great for dealing with the common situation of missing data):
berry_data |>
vis_miss(cluster = TRUE)
It looks like we have some big chunks of missing data (because of the way the data were collected: subjects did not have to attend every session). This is not a situation made-up for your benefit, just a characteristic of these data!
We are going to use these data to investigate the relationships among the scales, as well as predicting liking outcomes from categorical variables like age-category and gender.
Salmonella data (categorical predictor)
Today we’re also going to work with a subset of data from Michael Wesolowski’s 2017 thesis work on the effect of ethanol vapor on Salmonell spp. on produce. Michael ran a number of experiments, but we are going to consider just the experiment on the effect of exposure time on Salmonella CFUs on tomatoes. This is the first column of data in the Thesis Data Mastersheet.xls data file which you can download from a few places on Canvas, in particular the Files/Week 10 directory in the course setup package.
tomato_salmonella <-
read_excel("data/Week 10/salmonella-survival.xls",
range = "A6:I78")
skimr::skim(tomato_salmonella)| Name | tomato_salmonella |
| Number of rows | 72 |
| Number of columns | 9 |
| _______________________ | |
| Column type frequency: | |
| character | 3 |
| numeric | 6 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| Time (sec) | 68 | 0.06 | 1 | 2 | 0 | 4 | 0 |
| Replicate | 0 | 1.00 | 1 | 1 | 0 | 3 | 0 |
| Count | 0 | 1.00 | 1 | 5 | 0 | 30 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| Trial # | 0 | 1 | 2.00 | 0.82 | 1.00 | 1.00 | 2.00 | 3.00 | 3.00 | ▇▁▇▁▇ |
| Plate # | 0 | 1 | 1.50 | 0.50 | 1.00 | 1.00 | 1.50 | 2.00 | 2.00 | ▇▁▁▁▇ |
| Additional Dilution | 0 | 1 | 0.08 | 0.04 | 0.00 | 0.08 | 0.10 | 0.10 | 0.10 | ▂▁▁▁▇ |
| Total Dilution | 0 | 1 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | ▂▁▁▁▇ |
| Total CFU/mL | 0 | 1 | 1132885.96 | 2181060.41 | 909.92 | 909.92 | 909.92 | 450864.42 | 7279344.86 | ▇▁▁▁▁ |
| Log CFU/mL | 0 | 1 | 4.17 | 1.55 | 2.96 | 2.96 | 2.96 | 5.19 | 6.86 | ▇▁▂▁▃ |
Our Time (sec) variable–which is our key independent variable–seems to be missing a lot of data. If we looked at the nicely formatted Excel workbook that this comes from, this is just a formatting artifact: Michael didn’t want to repeat the 0, 5, etc exposure time for each row. We can fix this by using fill() to just fill in all empty columns with the last, non-empty value. This should be used with care (you don’t want to make fake data) but in this case it’s easier than editing the Excel file.
tomato_salmonella <-
tomato_salmonella |>
fill(`Time (sec)`)
tomato_salmonella# A tibble: 72 × 9
`Time (sec)` `Trial #` Replicate `Plate #` Count `Additional Dilution`
<chr> <dbl> <chr> <dbl> <chr> <dbl>
1 0 1 A 1 352 0.01
2 0 1 A 2 51 0.001
3 0 1 B 1 52 0.001
4 0 1 B 2 75 0.001
5 0 1 C 1 53 0.001
6 0 1 C 2 55 0.001
7 0 2 A 1 62 0.001
8 0 2 A 2 66 0.001
9 0 2 B 1 251.5 0.01
10 0 2 B 2 267.5 0.01
# ℹ 62 more rows
# ℹ 3 more variables: `Total Dilution` <dbl>, `Total CFU/mL` <dbl>,
# `Log CFU/mL` <dbl>The key independent variable in this experiment (one part of the thesis) is the time for which the tomatoes were exposed to ethanol vapor–Time (sec)–and the dependent (outcome variable) is the total number of colony forming units/mL, which we will express on a log scale–Log CFU/ml. The question Michael had was whether the treatment decreased the number of colony-forming units (CFUs) on the surface of the tomatoes. We can get an intuition for the answer to that question by looking at the mean and SD for each time exposure:
# A tibble: 4 × 3
`Time (sec)` mean sd
<chr> <dbl> <dbl>
1 0 6.60 0.220
2 10 3.25 0.671
3 15 3.32 0.686
4 5 3.53 0.831It certainly looks like any time exposure to ethanol vapor reduces the number of CFUs on the tomatoes. It isn’t so clear that longer exposures lead to further decreases, though.
Review of measurement levels
Typically, we talk about four “levels of measurement” in regards to data that dictate what kinds of analysis we should conduct (although we can do whatever we want, with sometimes disastrous results):
The borders between these categories is not as hard and fast as we’d wish. For example, the famous 9-pt hedonic scale (well, famous in my discipline!) is really an ordinal response scale:
We can have different kinds of variables in all different roles in our datasets. For example, in berry_data we have both nominal and ordinal predictors (independent or experimental variables):
- Berry type (nominal)
- Gender (nominal)
- Age (ordinal)
- Berry sample ID (nominal, even though it is listed as an integer)
- Subject ID (nominal, even though given as numeric)
We also have ordinal and interval-level outcomes (dependent or measured variables):
- 9-pt scale (ordinal, treated as interval)
- Labeled Affective Magnitude Scale (interval, although often claimed to be ratio)
- Visual Analog Scale (interval)
Technically, the choice of what statistical model we can fit to our data depends on the types of data we have! Frequently, we will violate the rules, but knowing that they exist will let us understand what risks we are taking.
With that, let’s talk about one of the most frequent situations we encounter in research: predicting the value of one interval/ratio-level variable from another. This is known as linear modeling (aka linear regression).
Linear modeling
Let’s begin by taking a quick look at the correlations among the various forms of measured liking for our berry_data:
9pt_overall us_overall lms_overall
9pt_overall 1.0000 0.2551 0.2573
us_overall 0.2551 1.0000 0.2687
lms_overall 0.2573 0.2687 1.0000# Then, we might plot them
berry_data |>
select(contains("overall")) |>
ggplot() +
geom_jitter(aes(x = .panel_x, y = .panel_y), alpha = 0.1) +
geom_smooth(aes(x = .panel_x, y = .panel_y), method = "lm", color = "red") +
geom_autodensity() +
facet_matrix(vars(1:3), layer.diag = 3) +
theme_classic()
We can see that while there are positive relationships among our variables, they are not as strongly correlated as we might assume for instruments that putatively measure the same underlying construct (overall liking). This is going to be our justification for exploring modeling what otherwise seems like it would be a kind of silly example.
We already learned that correlation (usually represented with the symbols \(r\) or \(\rho\)) represents the strength of linear association between two numeric variables. Correlations fall into the range \([-1, 1]\):
- \(-1\) indicates perfect anticorrelation: when \(x\) goes up, \(y\) goes down
- \(0\) indicates no relationship between the variables
- \(1\) indicates perfect correlation: when \(x\) goes up, \(y\) goes up
However, while correlation can give us an idea of the strength of a relationship, it doesn’t quantify the relationship. Informally (as I do most things in this class, because I’m not a statistician!), this is exactly what linear modeling does: in its simple form, linear regression quantifies the relationship between \(x\) and \(y\) whose strength is represented by the correlation \(r\).
Bivariate (2-variable) relationships imply a third variable
Before we begin to explore these models, a quick aside: these kinds of bivariate relationships always imply a third, implicit variable: the pairing variable, which is often give as \(i\), as in the index for pairs of variables: \((x_1, y_1), (x_2, y_2), ..., (x_i, y_i)\). This seems like something trivial to point out, but it is actually important and can be kind of puzzle.
For example, this is the problem with apparent (false) historical correlations:
The real, underlying relationship is between shark attacks and time, and between ice-cream consumption and time, but timepoint is being used as the pairing variable and so it creates a false correlation between the two other variables.
We often assume that our underlying pairing variable, \(i\), is a meaningless index, but it can often have consequences for our data. It is always worth at least explicitly stating how our data is structured to help us understand these (and other) kinds of errors we might encounter.
Hold up, why do I keep saying “linear model”?
The term “regression” from comes Sir Francis Galton. Galton was a late 19th-century brilliant scientist and pioneering statistician. He was also super committed to eugenics and what we could fairly call “scientific racism”. He sought to explain perceived differences among racial, ethnic, and class groups on a genetic basis. The story is long and ugly, and many of his scientific innovations–which are important and influential–were developed to support terrible goals and ideals.
The term “regression”–which comes from Galton’s paper Regression towards mediocrity in hereditary stature–clearly has an ideological bent built into it. The idea that “regression” modeling embodies is that, by understanding the mean of a population, we can make good linear predictions about it. This is why the common term we hear is “regression towards the mean”.
The more I learn about Galton and his view towards minoritized groups and others he thought were inferior to him, the weirder it feels that “regression”, a term that is freighted with meaning, has become the standard term we use for a method that could be more neutrally called the “linear model”, since it seeks to explain variation in \(y\) based on a linear model using \(x\). Knowing about this stuff–what the assumptions and biases that have shaped how we do science–is pretty important, and it is sometimes frustrating and unsettling.
So here I will often say “linear models”/“linear modeling” more and “regression” less, although I am sure I will flip between the two.
Simple linear models
With that all in mind, let’s state that the overall goal of linear regression is to obtain an estimate of the effect of one (usually continuous or interval/ratio) variable, \(x\), on the value of another (continuous) variable, \(y\). In order to estimate this relationship, we need some set of paired \(\{(x_i, y_i)\}\) so that we can estimate the relationship. We typically write this relationship as:
\[\hat{y}_i = \beta*x_i + \alpha\]
In plain English, we are saying that from each observed \(x_i\), we generate an estimate for \(y_i\), which we write as \(\hat{y}_i\), meaning it is an estimate that contains error, by multiplying \(x_i\) by some quantity \(\beta\) and adding a constant quantity \(\alpha\). Typically, we calculate both \(\alpha\) and \(\beta\) through Ordinary Least Squares, which gives us an exact numerical solution in order to minimize the quantity \(\sum{(y_i-\hat{y}_i)^2}\), which is the difference between our estimated and observed values.
I will not be deriving this method at all in this class.
I had to memorize the derivation (based in linear algebra) once, and I have mostly forgotten it. It is not too tricky, and it is interesting for the right sort of person, but it isn’t super pertinent to our everyday job of analyzing data as scientists. Instead, we’re going to take a look at how to run linear models in R, how to interpret and plot the output, and some basic considerations.
The command for the linear model in R (which includes both regression and ANOVA) is lm(). With our data, we can write the following:
# note that this is not a pipe-friendly function, need to use `data = ./_`
# depending on pipe
berry_lm <- lm(`9pt_overall` ~ lms_overall, data = berry_data)
berry_lm
Call:
lm(formula = `9pt_overall` ~ lms_overall, data = berry_data)
Coefficients:
(Intercept) lms_overall
5.34898 0.01254 class(berry_lm)[1] "lm"typeof(berry_lm)[1] "list"Notice that when we print berry_lm we just get a simple couple lines to the console: a list of the “function call” which tells us the model we fit, and estimates for the coefficients. In this case, we have only an (Intercept) coefficient, which is our \(\alpha\), and an lms_overall coefficient, which is our \(\beta\).
This doesn’t provide us with any information on whether this model is a good model. In order to get that information, we are going to do two things: examine a summary table of the output and make some visual plots to explore the model.
# We can get a pretty informative summary of our lm objects by using summary()
summary(berry_lm)
Call:
lm(formula = `9pt_overall` ~ lms_overall, data = berry_data)
Residuals:
Min 1Q Median 3Q Max
-5.6032 -1.5622 0.3626 1.6510 4.8927
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.348976 0.051718 103.43 <2e-16 ***
lms_overall 0.012542 0.001082 11.59 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.041 on 1893 degrees of freedom
(1455 observations deleted due to missingness)
Multiple R-squared: 0.06623, Adjusted R-squared: 0.06574
F-statistic: 134.3 on 1 and 1893 DF, p-value: < 2.2e-16This gives us the same information and more! We can see, for example, that a t-test run on the lms_overall \(\beta\) coefficient is apparently “significant” (more on this later), and we also can learn that despite this our model is pretty awful: the \(R^2\), which estimates for us how much variability in \(y\) is explained by our model and in simple regression is literally the squared correlation coefficient (\(r^2\)), is very low! Even though our model is a significant fit, we really aren’t explaining much at all about the variability of the results.
Why is this? Let’s take a look at our actual data:
berry_data |>
ggplot(aes(x = lms_overall, y = `9pt_overall`)) +
geom_jitter() +
# geom_smooth adds a fitted model estimate to the plot; note the
# "method = lm" call, where lm matches our `lm()` function
geom_smooth(method = 'lm') +
theme_classic()
We can see that, while there is indeed an underlying positive trend, it is quite weak: we have a lot of low values of 9pt_overall for high values of lms_overall, and vice versa.
We can use extractor functions to get some more useful information from our model to help us understand what’s going on. We use coefficients() to get the coefficients of the model (\(\alpha\) and \(\beta\) in this case) if we want to use them somewhere else (for example in bootstrapping, see below).
coefficients(berry_lm)(Intercept) lms_overall
5.34897610 0.01254223 We can extract predicted values for the model using the predict() function. We can also use this function to predict values for new \(x\) values which we want to estimate our expected \(y\) values for by providing them in the newdata = argument (as a vector). We can also get fitted (predicted) \(\hat{y}\) values using the fitted() function.
# A tibble: 1,895 × 2
observation predicted_9pt
<chr> <dbl>
1 7 5.85
2 8 5.49
3 9 5.24
4 10 5.49
5 11 4.92
6 12 4.11
7 13 6.28
8 14 6.28
9 15 6.05
10 16 6.05
# ℹ 1,885 more rowsWe can use this to take a look at our model’s effectiveness visually. We’re going to
- Select the observed
9pt_overallvalues from our data - Get the predicted observations from our model
- Join them into a single data table
- Plot them against each other
If the resulting plot is close to a diagonal line, we’ll know that our model is performing well, because ideally we want the error (\(y_i-\hat{y}_i\)) to be very small!
predicted_y |>
# we create a new observation variable to filter our left_join
left_join(berry_data |> mutate(observation = as.character(row_number()))) |>
select(observation, `9pt_overall`, predicted_9pt) |>
ggplot(aes(x = `9pt_overall`, y = predicted_9pt)) +
geom_jitter() +
geom_smooth(method = "lm", color = "red") +
theme_classic() +
labs(y = "Predicted 9-pt score", x = "Observed 9-pt score")
Unfortunately, it is clear that our model really doesn’t do much to predict the variation in our outcome variable. Equivalently, there is only a weak relationship between our predicted \(\hat{y}_i\) values and our actual \(y_i\) values.
This is just one of the interesting things we found out in this study!
Bonus: what’s the underlying (implied) variable(s) in this model?
After spending a little bit of time on thinking about implicit links caused by the indexing variable, we immediately skimmed over what \(i\) means in this study Looking back at our skim(), what variable(s) links the observations of 9pt_overall and lms_overall?
# Give your answers here!What does something like “statistical significance” actually mean?
Before we dive into our example, let’s pause for a minute. When we used summary(berry_data), we saw that the effect of lms_overall on 9pt_overall was “significant”. This sounds like it’s important. But is it?
This is a bit of a loaded question, and one that I have gotten less comfortable with (in a technical sense!) since I began to teach this material.
Briefly, you’ve learned about null-hypothesis significance testing (NHST) and associated terms like the “null hypothesis”, “\(p\)-values”, and of course “significance” in your stats class. These are terms that come from a particular school of statistical theory called “frequentist” statistics. Without derailing this class too much, I want to point out that these are increasingly no longer recognized as the best tools for the job, but they are very common tools in scientific data analysis, so we’ll be relating much of what we do to these tools.
Technically, when we speak about a result being statistically significant, we mean something like
The probability of observing results as extreme or more extreme than the ones observed, given that all of our assumptions are correct, is less than a pre-selected chance of a false positive (typically called \(\alpha\); the chance of a Type I error).
This means that we “reject the null hypothesis”, because it seems unlikely given these results. Two important corollaries that are often misunderstood come from this type of reasoning. The following statements are about what NHST does not tell us:
- A significant result does not prove an alternative hypothesis.
- A non-significant result does not prove the null hypothesis.
We will return to this topic when we discuss critical thinking with data in a few weeks. But I want to just make sure I define what we’re doing here so that you can understand what we are doing–and what we’re not doing.
A resampling approach to estimating quality of regression
Given how bad our model explanation is, and how little effect lms_overall seems to have (we could read our model as saying that our estimated 9pt_overall score increases by 0.013 points for every 1 point increase in lms_overall), it is curious that we think this effect is important at all.
First, we might want to ask about the stability of our estimate. We can use this via bootstrapping, just as we did for means and variances in the last class.
# As a warm up, we need to remember we need a place to store our data
bootstrapped_betas <- numeric(1000)
# We may as well make ourselves a progress bar so we know the computer isn't
# frozen
pb <- progress::progress_bar$new(total = 1000)
# Then we run a loop for some large number of times (i) that will get us plausible, sample-based estimates of our statistic of interest: the
for(i in 1:1000){
# Resample our data with replacement
bootstrap_data <- berry_data[sample(1:nrow(berry_data), replace = TRUE), ]
# run the same model on the new sample
bootstrap_model <- lm(`9pt_overall` ~ lms_overall, data = bootstrap_data)
# Get the lms_overall (beta) coefficient
bootstrapped_betas[i] <- coefficients(bootstrap_model)[2]
# Update our progress bar
pb$tick()
}
# We'll grab these quantile values to make a nice confidence interval
beta_limits <- quantile(bootstrapped_betas, c(0.05, 0.95))
# Now we can plot the distribution of these beta coefficients and compare them to our observed coefficient
bootstrapped_betas |>
as_tibble() |>
ggplot(aes(x = value)) +
geom_histogram(bins = 20, fill = "grey", color = "black") +
# This is a hacky way to draw a single line segment
annotate(geom = "segment",
x = beta_limits[1], xend = beta_limits[2],
y = 0, yend = 0,
color = "red", linewidth = 2) +
# And here we plot our observed beta coefficient
geom_point(aes(x = coefficients(berry_lm)[2], y = 0),
color = "red", fill = "white",
size = 4,
shape = 21) +
theme_classic() +
labs(x = expression(beta), y = "count")
# The limits of the confidence interval are
beta_limits 5% 95%
0.01057674 0.01427106 There is an almost 40% difference between our bootstrapped confidence-interval’s upper and lower bounds, which is not super tight. But we see that, for example, 0 is not a plausible value from any of our resamples. The order of magnitude of our estimate doesn’t change very much acrosss simulations–we’re reasonably stable in absolute estimates, although relatively we change a fair amount.
We could also use a permutation approach to consider the significance of our estimate. When we break the implied pairing between our two variables (disrupt the \(i\) variable), what do we learn? We will only run 1000 of the almost infinitely large possible reorderings of the \(i\) variable:
# Don't forget to store your results somewhere
permuted_betas <- numeric(1000)
# And it's also always good practice to make sure we're not stuck just because
# we don't see progress
pb <- progress::progress_bar$new(total = 1000)
for(i in 1:1000){
# Scramble our data (not just resample)
permuted_data <- tibble(`9pt_overall` = sample(berry_data$`9pt_overall`, replace = FALSE),
lms_overall = berry_data$lms_overall)
# run the same model on the new sample
bootstrap_model <- lm(`9pt_overall` ~ lms_overall, data = permuted_data)
# Get the lms_overall (beta) coefficient
permuted_betas[i] <- coefficients(bootstrap_model)[2]
pb$tick()
}
permuted_betas |>
as_tibble() |>
ggplot(aes(x = value)) +
geom_histogram(bins = 20, fill = "grey", color = "black") +
theme_classic() +
labs(x = expression(permuted~beta~coefficient), y = "count") +
annotate(geom = "point", x = coefficients(berry_lm)[2], y = 0,
color = "red", size = 4, shape = 21, fill = "white")
Happily, it does appear that the relationship, while weak, is indeed stronger than we’d expect to observe if there were no relationship between our subjects’ reported liking on 2 scales! That’s a relief.
So, can you explain the difference between what we are doing for the bootstrap and for the permutation? What different questions are we answering?
# Let's write some possible responses!Checking your model
We’ve spent a little bit of time talking about how to interpret a model and how to investigate whether it is telling us anything useful. But an important facet of using any statistical model (that I am barely going to touch on) is investigating the fit of the model. The lm() function tries to help you out with this: by using plot() on an lm object, you can see a set of plots that give you information about model fit.
plot(berry_lm)
{width=c(“50%”, “50%”, “50%” “50%”)}
{width=c(“50%”, “50%”, “50%” “50%”)}
{width=c(“50%”, “50%”, “50%” “50%”)}
{width=c(“50%”, “50%”, “50%” “50%”)}
Each of these plots gives you information about how well the model fits. In this set, I tend to look most at the first two plots: the fitted values against the residuals and the QQ-plot of theoretical quantiles vs standardized residuals.
- The first plot tells us whether there is a systematic pattern in (mis)-fitting our data. The question is whether the model fits worse (bigger absolute residuals) at some set of fitted values than at others. It does appear that there is a systematic relationship between fitted values and residuals (\(y_i - \hat{y_i}\)), which tells us that this model is not unbiased: it will tend to predict more extreme values at either end of the range of our predicted values.
- The second plot tells us something about the normality of our data. It orders the residuals from lowest to highest, and then assigns them to quantiles of the normal distribution. Then, it plots those assigned quantiles against the actual normalized residuals. In a well-fitting model, the residuals should cling pretty close to a straight line–in our model, we see the same systematic deviation at both low and high levels (in this case represented as quantiles). This is an indicator of non-normality in our data, which in turn means our estimates are probably bad.
The last 2 plots give more technical information about the model fit (and look wild for this model), and I will not go into them here.
Beyond these simple plots, there are many important different ways to assess a model. But we will not go into them here, because we start to get into statistical assumptions, and I think you should learn about these methods from someone who is more expert in them. My goal is to equip you with the tools to manipulate data and execute these analyses, and interpretation is best left to the domain experts!
A note about causality
We are all familiar with the phrase “correlation does not imply causation”. But why do we say this so frequently?
Without getting too into the weeds, we know a few things that should make us cautious!
- Correlation is symmetrical:
cor(x, y) == cor(y, x). So we can’t say anything about what is causing what! - Correlation is a binary relationship: we only look at
cor(x, y), so if \(z\) is causing both \(x\) and \(y\), we will see a strong correlation but won’t know why. - Human brains seem to be simultaneous pattern-recognition and hypothesis-generation machines. So when we see associations (like correlations) we immediately make up plausible theories for them. We need to be cautious about our intuitions in the presence of these kinds of data explorations!
A coding digression: factor
OK. We’ve done a bit of inference. Let’s take a minute before we move on from linear models with continuous variables to those with categorical variables–typically called “ANOVA” models.
So far we have talked about using one continuous variable to explain/predict variation in another continuous variable. Often, however, we will have categorical variables whose effect on some outcome variable \(y\) we want to explore. In this case, we are going to end up working with a king of R data type we haven’t really explore yet: factor.
I always recommend reading in the weekly assignment, but I am going to call it out here as well to remind you to check these out: the R4DS chapter on factors (Wickham and Grolemund 2017) is probably the best resource for this class. The links at the “Learning More” section of that chapter are great resources for those of you interested in the programming history of R
The short story behind the existence of factors is that they are memory- and computationally-efficient ways to deal with limited sets of (possible non-numeric) labels, which to us are strings like "male", "female", "red“, etc. In early R history it was not really necessary to deal with strings much, and so it was easier to largely treat them as labels for underlying integers. This is something we all do frequently in our own data-management without thinking about: think about how in an Excel lab notebook we fill the cells in a column like”Income level” by writing a value like “2”, which represents something like “$20-$30,000/year”. This is the idea behind how R treats factor data.
Factors also have the advantage that they are arbitrarily ordered, so that we can tell R to ignore things like alphabetical or numerical order.
A factor data type is a set of (small) integers with an attribute tag of levels that correspond to those labels. Let’s examine these in the context of the salmonella data.
time_steps <- tomato_salmonella$`Time (sec)`
unique(time_steps) # we don't want to print all 3000+ duplicates[1] "0" "5" "10" "15"class(time_steps)[1] "character"typeof(time_steps)[1] "character"Note that because of the way in which the data were imported, these are characters (and not integers). We could use as.integer() to convert them, but we don’t necessarily want to do math with these. We want to use them as discrete labels: in this experiment, there were no 4-second, 13-second, or 23-second intervals. We only have 4 categories. This is what a factor does.
In the language of factor, the label values are called levels.
Factor w/ 4 levels "0","5","10","15": 1 1 1 1 1 1 1 1 1 1 ...typeof(time_factor)[1] "integer"class(time_factor)[1] "factor"attributes(time_factor)$levels
[1] "0" "5" "10" "15"
$class
[1] "factor"We can get the levels back out of a factor by using levels()
levels(time_factor)[1] "0" "5" "10" "15"When we try to add a non-included level to a factor, R automatically coerces it to an NA.
Thus, a factor is R’s way of creating a map between (computationally expensive) strings and inexpensive, easy to manipulate integers. Many modern (tidyverse-type) functions no longer use factors as extensively, but a lot of base R functions do use them (in particular the aov() and lm() functions we are using a lot today). They either will request factors or convert strings to factors in the background silently.
Another place that factors can be quite useful is in ggplot2. If we have an integer label (like sample in berry_data, which represents a sample #), by default it will be treated as numeric:
berry_data |>
ggplot(aes(x = lms_overall, y = `9pt_overall`, color = sample)) +
geom_jitter() +
theme_classic()Warning: Removed 1455 rows containing missing values or values outside the scale range
(`geom_point()`).
Notice the gradient assigned to the color mapping for sample.
If we tell R that sample is a factor (by using the as.factor() coercion function), ggplot() will recognize this and give us a much more useful (discrete) color mapping:
berry_data |>
mutate(sample = as.factor(sample)) |>
ggplot(aes(x = lms_overall, y = `9pt_overall`, color = sample)) +
geom_jitter() +
theme_classic()Warning: Removed 1455 rows containing missing values or values outside the scale range
(`geom_point()`).
Now that you’ve been introduced to factor data, you may have questions about, for example, how to add or subtract levels to a factor, or how to reorder those levels.
The forcats package, part of the tidyverse, has a lot of useful functions for doing these tasks, all of which start with fct_*(). This makes it really easy to search for the one we want with tab-completion. How do we do this?
Please consult the R4DS chapter on factors (Wickham and Grolemund 2017) to learn more–it is easy but slightly fiddly. And now, if you encounter weird effects in your treatment of string data, you’ll know why!
Linear models with categorical predictors (Analysis of variance/ANOVA)
Now that we’ve introduced a new type of data in R, let’s justify it’s use! Often, we want to explain the variation in a (continuous) \(y\) variable by the effect of a categorical \(x\) variable–in base R like lm() these are usually treated (or coerced to) factor data. Happily, unless we want to make changes to levels, we can let R handle this behind the scenes.
When we have some \(x\) variable that can only take on a set of discrete values, this is a categorical predictor, and what has historically been called linear regression becomes instead what has historically been called Analysis of Variance: “ANOVA”.
Just like regression, ANOVA has a somewhat fraught history. Sir Ronald Fisher, the inventor of ANOVA (and a host of other major statistical methods and theories) was another eugenicist and promoter of what we would now consider scientific racism. Some have argued that we should reconsider terms like “variance” because of their ties with ideas of “deviation from the norm”–another obviously loaded phrase–but I wonder if this is more effort than its worth. “Regression” seems, somehow, more obviously contentious than “variance”, so I will probably keep trying to reduce my use of “regression”, while not worrying as much about “ANOVA”. YMMV!
In ANOVA, we are interested in the effect of category membership \(\{x_1, x_2, ..., x_n\}\) on our predicted outcome variable \(y\). Instead of our classic linear model equation (recall \(y = a + bx\)), we tend to write ANOVA something like the following:
\[\hat{y}_i = \mu + \tau_k\]
We are familiar with the LHS of this equation–our estimate for a particular, observed \(y_i\) is indicated as \(\hat{y}_i\). The RHS looks a bit different: it is read as “the grand mean” (\(\mu\)) “plus the group mean of the level of \(x_i\)” (which is written here as \(\tau_k\), because we often talk about these as treatment means).
This bit of mathiness can be rephrased intuitively as “our estimate of a particular \(y\)-value is the sum of the average of \(y\) plus the average effect of the treatment variable \(x\)”. Even though this sounds like a different problem, it’s actually our same equation from linear regression reframed for a situation in which \(x\) can only take on a delimited set of values. We have an example(s) of this in our tomato_salmonella data that will give us some context in this question.
To put it another way, ANOVA relies on the simple observation that, if some grouping variable (“treatment”) is important, there will be less variation within groups than there is between different groups. This intution is shown here (with the assumption of normality that is important for linear-model ANOVA statistics):
One-way ANOVA
In its simplest form, ANOVA is used to relate a single outcome \(y\) and a single predictor \(x\). We still have the implicit \(i\) variable that links our pairs of observations \((x_i, y_i)\), but now \(x_i \in \{x_1, x_2, ..., x_n\}\)–it does not have an unlimited number of values it can take on, as it does in standard linear regression. Thus, our \(\beta\) from linear regression is replaced by a set of \(\tau_k\) values that represent the average effect for each possible group \(x_i\) can be in. But we could as easily write \(\beta_k\)–the meaning is the same. It is just convention.
In tomato_salmonella, Time (sec) is a categorical variable that represents the time the tomato was exposed to ethanol vapor: it can be one of four values.
unique(tomato_salmonella$`Time (sec)`)[1] "0" "5" "10" "15"Let’s model Log CFU/ml (log colony-forming units per ml) on Time (sec). This model asks: does the exposure time predict (reduce, one would assume) the number of logCFUs of salmonella?
Call:
lm(formula = `Log CFU/mL` ~ time_factor, data = tomato_salmonella)
Coefficients:
(Intercept) time_factor5 time_factor10 time_factor15
6.605 -3.078 -3.357 -3.289 summary(tomato_aov)
Call:
lm(formula = `Log CFU/mL` ~ time_factor, data = tomato_salmonella)
Residuals:
Min 1Q Median 3Q Max
-0.5676 -0.3565 -0.2888 0.1736 1.6088
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.6050 0.1518 43.50 <2e-16 ***
time_factor5 -3.0784 0.2147 -14.34 <2e-16 ***
time_factor10 -3.3571 0.2147 -15.63 <2e-16 ***
time_factor15 -3.2894 0.2147 -15.32 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.6442 on 68 degrees of freedom
Multiple R-squared: 0.8348, Adjusted R-squared: 0.8275
F-statistic: 114.5 on 3 and 68 DF, p-value: < 2.2e-16Notice we can just use our standard lm() function, even though our \(x\) variable, time_factor, is categorical. The model object and summary() for tomato_aov has coefficients for 3/4 exposure times–this is because lm() needs to have an intercept value (the \(a\) in \(y = a + bx\)), which in ANOVA is the “reference category”. Since 0 the first level (we defined it that way with the levels = statement in the factor() function), the coefficients in the model are the values for a change from 0 (the average CFU count for 0 is the (Intercept) term). So 15 seconds of exposure sees a drop of the most average CFUs, whereas 5 sees the least. We also see our typical t-test estimates for significance for each level, but these are a bit harder to interpret for this model because these are actually each the same \(x\) variable.
Therefore, we might want to use the anova() function, which is a version of summary() specific for ANOVA.
anova(tomato_aov)Analysis of Variance Table
Response: Log CFU/mL
Df Sum Sq Mean Sq F value Pr(>F)
time_factor 3 142.623 47.541 114.55 < 2.2e-16 ***
Residuals 68 28.222 0.415
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1This tells us that the \(x\) variable time_factor is overall “significant”, but doesn’t tell us about specific factor levels. We can get to this from the very start by using the aov() function, which is a version of lm() for ANOVA:
tomato_aov_2 <-
tomato_salmonella |>
aov(`Log CFU/mL` ~ time_factor, data = _)
tomato_aov_2Call:
aov(formula = `Log CFU/mL` ~ time_factor, data = tomato_salmonella)
Terms:
time_factor Residuals
Sum of Squares 142.62330 28.22201
Deg. of Freedom 3 68
Residual standard error: 0.6442279
Estimated effects may be unbalancedsummary(tomato_aov_2) Df Sum Sq Mean Sq F value Pr(>F)
time_factor 3 142.62 47.54 114.5 <2e-16 ***
Residuals 68 28.22 0.42
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The typical way to report this kind of model is “there was a significant effect of berry type on overall liking as measured by an unstructured line scale”. But can we do better?
Resampling approaches to categorical data
Just like with our continuous linear model, we have a couple questions we’d like to answer about our proposed model. First, we’d like to know whether the observed effect on treatments–the different variances that give us a large \(F\)-ratio for the ANOVA table above–are really that large. Second, we’d like to know something about our individual treatment levels: for example. is the effect of 15 seconds really higher than that of 20 seconds? That’s a bit weird!
Permutation for whether an effect is real
We will answer the first question with a permutation test: we will break the relationship between the groups and the samples, and see if we can find as large an \(F\)-ratio as we observed in our actual data.
# first, just to show how this works, let's do it once:
permutation_example <- tomato_salmonella |>
select(`Time (sec)`, `Log CFU/mL`) |>
# unlike bootstrapping we use replace = FALSE to just permute our sample
bind_cols(sample(tomato_salmonella$`Time (sec)`, replace = FALSE)) |>
rename(actual_time = 1, logCFU = 2, permuted_time = 3)
permutation_example# A tibble: 72 × 3
actual_time logCFU permuted_time
<chr> <dbl> <chr>
1 0 6.51 10
2 0 6.67 5
3 0 6.68 5
4 0 6.83 5
5 0 6.68 5
6 0 6.70 15
7 0 6.75 15
8 0 6.78 5
9 0 6.36 5
10 0 6.39 15
# ℹ 62 more rowsWe can then compare our observed means to our permuted means:
# A tibble: 4 × 2
actual_time mean
<chr> <dbl>
1 0 6.60
2 10 3.25
3 15 3.32
4 5 3.53# A tibble: 4 × 2
permuted_time mean
<chr> <dbl>
1 0 3.57
2 10 3.71
3 15 4.39
4 5 5.02It is immediately apparent that the difference between the groups is less extreme in or permuted example. This makes sense: some observations that actually came from the high logCFU 0-second treatment have now been assigned to other times, while some from the low logCFU treatments have been assigned to 0 seconds. But just a single comparison isn’t enough to prove anything. What we will do is repeat this experiment thousands of times:
permutation_reps <- 1000 # set this to 100 or so if your computer is slow for the demo
# Here's one way we can store our results
permuted_F_ratio <- tibble()
# It's always nice to create a progress bar
pb <- progress::progress_bar$new(total = 1000)
for(i in 1:1000){
permuted_data <-
# Make a new table of permuted treatments by combining the scrambled times
# with the (unscrambled) logCFUs. You could scramble these too, but there's
# no real need. Why?
tibble(time = sample(tomato_salmonella$`Time (sec)`, replace = FALSE),
logCFUs = tomato_salmonella$`Log CFU/mL`)
perm_aov <-
permuted_data |>
aov(logCFUs ~ time, data = _) |>
# We use the `broom` package, which turns (linear) model outputs into tidy
# data frames to get our F statistic
broom::tidy()
# don't forget to store the results!
permuted_F_ratio <-
bind_rows(permuted_F_ratio,
perm_aov)
# Increment the progress bar
pb$tick()
}Now we can ask the question we care about of these simulated data: how many times do we observe a ratio of variances that is larger than that in our real data? How “real” is our effect?
# A tibble: 1 × 1
mean_permuted_F
<dbl>
1 1.04permuted_F_ratio |>
drop_na() |>
ggplot(aes(x = statistic)) +
geom_histogram(color = "black", fill = "grey") +
geom_point(data = tomato_aov_2 |> broom::tidy(),
mapping = aes(y = 0),
color = "red") +
theme_classic() +
labs(x = expression(italic(F)*-statistic), y = NULL)
As is evident, we never observe a value of the \(F\)-statistic that is as extreme as the one we observe in our actual data. This is pretty strong evidence that there is indeed a relationship between the treatment (time exposure) and the outcome of interest (Salmonella logCFUs).
Posthoc testing
In ANOVA, we typically are interested in whether there are significant differences among our treatment levels: sure, we have estimates for \(\tau_{0s}\), \(\tau_{5s}\), etc, but are they actually different from each other?
To answer these questions, we often use what are called post-hoc tests, which have names like “Fisher’s Least Significant Difference” (Fisher’s LSD) or “Tukey’s Honestly Significant Difference” (Tukey’s HSD). I think the easiest way to calculate these are from the (I no longer recommend this as a first choice, see the section on this package for how to use it.)agricolae package.
Currently, I use the emmeans (Estimated Marginal Means) package for post-hoc comparisons. emmeans does a good job of allowing multiple comparisons even in the presence of interactions, more complex models, etc. When combined with the multcomp package, it can also give us easy outputs of “compact letter displays”, which we’ll examine below.
time_factor emmean SE df lower.CL upper.CL
0 6.60 0.152 68 6.30 6.91
5 3.53 0.152 68 3.22 3.83
10 3.25 0.152 68 2.94 3.55
15 3.32 0.152 68 3.01 3.62
Confidence level used: 0.95 # We can do all pairwise comparisons using `pairs()`;
# default adjustment is Tukey but others are available
tomato_emm |>
pairs() contrast estimate SE df t.ratio p.value
time_factor0 - time_factor5 3.0784 0.215 68 14.335 <0.0001
time_factor0 - time_factor10 3.3571 0.215 68 15.633 <0.0001
time_factor0 - time_factor15 3.2894 0.215 68 15.318 <0.0001
time_factor5 - time_factor10 0.2788 0.215 68 1.298 0.5673
time_factor5 - time_factor15 0.2111 0.215 68 0.983 0.7597
time_factor10 - time_factor15 -0.0677 0.215 68 -0.315 0.9891
P value adjustment: tukey method for comparing a family of 4 estimates # the multcomp::cld() function gives us the "compact letter display" we want
tomato_tukey <-
tomato_emm |>
cld(Letters = letters)
tomato_tukey time_factor emmean SE df lower.CL upper.CL .group
10 3.25 0.152 68 2.94 3.55 a
15 3.32 0.152 68 3.01 3.62 a
5 3.53 0.152 68 3.22 3.83 a
0 6.60 0.152 68 6.30 6.91 b
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 4 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. The key output of these kinds of tests are the “letters”–the groups labeled with the same letter are close enough to each other based on the test’s criteria to be considered not significantly different.
All of these methods have various responses to the problem of familywise error: when we run multiple tests on the same set of data, we are increasing our chance of finding a relationship that is just coincidental. This is often referred to as “p-hacking” or “significance mining”. We will return to this issue in Week 14, but this XKCD actually very aptly sums up the problem:
I will say this: do not use Fisher’s LSD. It is frequently used in published papers, and it is very very prone to this kind of error. While there is no perfect approach, Tukey’s HSD (and other methods) are less likely to give this kind of false positive.
As an aside, we can use cld() to get us nice tables for output to Word, PDF, whatever if we wrangle it a bit. First make sure everything is correctly saved in your output. Then:
| time_factor | emmean | SE | df | lower.CL | upper.CL | .group |
|---|---|---|---|---|---|---|
| 0 | 6.604979 | 0.151846 | 68 | 6.301975 | 6.907983 | b |
| 5 | 3.526595 | 0.151846 | 68 | 3.223591 | 3.829599 | a |
| 10 | 3.247841 | 0.151846 | 68 | 2.944837 | 3.550845 | a |
| 15 | 3.315539 | 0.151846 | 68 | 3.012535 | 3.618543 | a |
tomato_cld_tibble |>
ggplot(aes(x = time_factor, y = emmean)) +
geom_point() +
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
size = 1/4, width = 1/10) +
geom_text(aes(y = upper.CL + 0.1, label = .group))Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
Or even better, you can write this to a CSV to import into Word, whatever you want:
tomato_cld_tibble |>
write_csv("your/directory/here/output.csv")Alternative post-hoc approach: {agricolae} package {agricolae}
What I used to focus on:
Study: tomato_aov_2 ~ "time_factor"
LSD t Test for Log CFU/mL
Mean Square Error: 0.4150296
time_factor, means and individual ( 95 %) CI
Log.CFU.mL std r se LCL UCL Min Max
0 6.604979 0.2202668 18 0.151846 6.301975 6.907983 6.200798 6.862092
10 3.247841 0.6712713 18 0.151846 2.944837 3.550845 2.959002 4.856629
15 3.315539 0.6864083 18 0.151846 3.012535 3.618543 2.959002 4.592471
5 3.526595 0.8305656 18 0.151846 3.223591 3.829599 2.959002 4.816335
Q25 Q50 Q75
0 6.377689 6.679142 6.778546
10 2.959002 2.959002 2.959002
15 2.959002 2.959002 2.959002
5 2.959002 2.959002 4.460705
Alpha: 0.05 ; DF Error: 68
Critical Value of t: 1.995469
least Significant Difference: 0.4285123
Treatments with the same letter are not significantly different.
Log CFU/mL groups
0 6.604979 a
5 3.526595 b
15 3.315539 b
10 3.247841 bHSD.test(y = tomato_aov_2, trt = "time_factor", console = TRUE)
Study: tomato_aov_2 ~ "time_factor"
HSD Test for Log CFU/mL
Mean Square Error: 0.4150296
time_factor, means
Log.CFU.mL std r se Min Max Q25 Q50
0 6.604979 0.2202668 18 0.151846 6.200798 6.862092 6.377689 6.679142
10 3.247841 0.6712713 18 0.151846 2.959002 4.856629 2.959002 2.959002
15 3.315539 0.6864083 18 0.151846 2.959002 4.592471 2.959002 2.959002
5 3.526595 0.8305656 18 0.151846 2.959002 4.816335 2.959002 2.959002
Q75
0 6.778546
10 2.959002
15 2.959002
5 4.460705
Alpha: 0.05 ; DF Error: 68
Critical Value of Studentized Range: 3.724638
Minimun Significant Difference: 0.5655713
Treatments with the same letter are not significantly different.
Log CFU/mL groups
0 6.604979 a
5 3.526595 b
15 3.315539 b
10 3.247841 bBootstrapping for post-hoc comparisons
A better approach than either of these, I think, is a bootstrapping approach.
# Find some place to keep the results
tomato_boots <- tibble()
# Be nice, keep track of your progress
pb <- progress::progress_bar$new(total = 1000)
for(i in 1:1000){
boot_results <-
tomato_salmonella |>
# This is the key step where we get a sample of the same size WITH
# REPLACEMENT
group_by(time_factor) |>
slice_sample(prop = 1, replace = TRUE) |>
summarize(mean_CFUs = mean(`Log CFU/mL`)) |>
mutate(boot_rep = i)
# Don't forget to store your results
tomato_boots <-
bind_rows(
tomato_boots,
boot_results
)
# And increment the progress bar
pb$tick
}
tomato_boots# A tibble: 4,000 × 3
time_factor mean_CFUs boot_rep
<fct> <dbl> <int>
1 0 6.53 1
2 5 3.50 1
3 10 3.14 1
4 15 3.50 1
5 0 6.58 2
6 5 3.51 2
7 10 3.27 2
8 15 3.31 2
9 0 6.55 3
10 5 3.39 3
# ℹ 3,990 more rowsThis gets us 1000 “new” experiments, with the means for each treatment. We can then ask ourselves about the stability and confidence intervals around our estimated means, to see which are indeed different from each other. This is a “non-parametric” HSD-type approach.
tomato_averages <-
tomato_salmonella |>
group_by(time_factor) |>
summarize(mean_CFUs = mean(`Log CFU/mL`))
tomato_boots |>
ggplot(aes(x = time_factor, y = mean_CFUs)) +
geom_boxplot() +
geom_point(data = tomato_averages, color = "red") +
theme_classic() +
labs(x = "Ethanol Exposure", y = "Average logCFUs")
We can see that the effect of any ethanol treatment is very consisten across resamplings, and that the 5 second exposure has a consistently higher logCFU load than the 10 and 15 second exposures, but that the other 2 are not consistently different.
We might consider plotting these as densities, which would give us a look at overlap a little more easily.
tomato_boots |>
ggplot(aes(x = mean_CFUs, fill = time_factor)) +
geom_density(alpha = 2/3) +
theme_bw() +
labs(fill = "Ethanol Exposure", x = "Average logCFUs", y = NULL) +
scale_fill_brewer(palette = "Accent")
It is clear that there is a strong effect from any level of ethanol exposure, but that the different levels are not strongly differentiated.
Plotting ANOVA results
Because ANOVA is based on categorical data, it doesn’t look very good to plot that data in scatterplots. Generally, boxplots or barplots will be the most effective visual representation of data from ANOVA. We go back to the berry data but use some of the categories (berry types) to demonstrate some other approaches.
The following visualizations were adapted to some degree from Philipp Masur’s blog:
berry_data |>
ggplot(aes(x = berry, y = us_overall)) +
geom_boxplot(aes(fill = berry), notch = TRUE) +
geom_jitter(alpha = 0.1) +
theme_classic()Warning: Removed 790 rows containing non-finite outside the scale range
(`stat_boxplot()`).Warning: Removed 790 rows containing missing values or values outside the scale range
(`geom_point()`).
berry_data |>
ggplot(aes(x = berry, y = us_overall)) +
stat_summary(aes(fill = berry), fun.y = mean, geom = "bar") +
stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = 0.25) +
theme_classic()Warning: The `fun.y` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
ℹ Please use the `fun` argument instead.Warning: Removed 790 rows containing non-finite outside the scale range
(`stat_summary()`).
Removed 790 rows containing non-finite outside the scale range
(`stat_summary()`).
Another classic plot for ANOVA that will become even more powerful when we begin to examine multiway analyses is the interaction plot, which is really a dot + line plot representing the same data as the bar plot:
berry_data |>
# We will make our own summary table rather than do transformations in ggplot
group_by(berry) |>
summarize(mean_liking = mean(us_overall, na.rm = TRUE),
se_liking = sd(us_overall, na.rm = TRUE)/sqrt(n())) |>
ggplot(aes(x = berry, y = mean_liking, color = berry, group = berry)) +
geom_line(aes(group = 1), color = "black") +
geom_point(aes(y = mean_liking)) +
geom_errorbar(aes(ymin = mean_liking - 2 * se_liking, ymax = mean_liking + 2 * se_liking), width = 0.25) +
theme_classic()
Incorporating multiple predictors
As a tease of our material for next week, we can always include multiple predictors in linear regression and ANOVA. This lets us simultaneously estimate the effect of multiple predictors on some outcome. So, just for fun, let’s go back to ANOVA for the berries and estimate the effect of berry type and subject age on over all berry liking.
Df Sum Sq Mean Sq F value Pr(>F)
berry 3 849 283.10 19.970 8.78e-13 ***
age 2 188 93.77 6.615 0.00137 **
berry:age 6 51 8.52 0.601 0.72949
Residuals 2388 33852 14.18
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
950 observations deleted due to missingnessWe can see that there is an overall effect of berry that remains even when the significant effect of age is also in the model. There doesn’t appear to be an interactive effect of berry and age, as shown by the third line for berry:age. We’ll talk about this in more detail tomorrow, but this plot should give an idea of what that means:
berry_data |>
ggplot(aes(x = age, y = us_overall, color = berry)) +
geom_boxplot(notch = TRUE) +
theme_classic()Warning: Removed 790 rows containing non-finite outside the scale range
(`stat_boxplot()`).
Often, when model multiple regression and ANOVA, we are not only interested in predicting \(y\), but in determining which different \(x\)s matter. This is the topic we will explore next week!
Reading
This week, you should read:
- Practical Statistics for Data Scientists, Chapter 3 (especially on ANOVA), Chapter 4 (on Simple Regression)
- R for Data Science, Chapters 22-24
Session info
R version 4.5.3 (2026-03-11)
Platform: aarch64-apple-darwin20
Running under: macOS Tahoe 26.5.2
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: America/New_York
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] agricolae_1.3-7 multcomp_1.4-30 TH.data_1.1-5 MASS_7.3-65
[5] survival_3.8-6 mvtnorm_1.3-7 emmeans_2.0.3 readxl_1.4.5
[9] ggforce_0.5.0 naniar_1.1.0 skimr_2.2.2 lubridate_1.9.5
[13] forcats_1.0.1 stringr_1.6.0 dplyr_1.2.1 purrr_1.2.2
[17] readr_2.2.0 tidyr_1.3.2 tibble_3.3.1 ggplot2_4.0.2
[21] tidyverse_2.0.0
loaded via a namespace (and not attached):
[1] tidyselect_1.2.1 farver_2.1.2 S7_0.2.1
[4] fastmap_1.2.0 tweenr_2.0.3 rpart_4.1.24
[7] digest_0.6.39 timechange_0.4.0 estimability_1.5.1
[10] lifecycle_1.0.5 cluster_2.1.8.2 multcompView_0.1-11
[13] magrittr_2.0.5 compiler_4.5.3 Hmisc_5.2-6
[16] rlang_1.2.0 progress_1.2.3 tools_4.5.3
[19] utf8_1.2.6 yaml_2.3.12 data.table_1.18.2.1
[22] knitr_1.51 prettyunits_1.2.0 labeling_0.4.3
[25] htmlwidgets_1.6.4 bit_4.6.0 repr_1.1.7
[28] RColorBrewer_1.1-3 foreign_0.8-91 withr_3.0.2
[31] nnet_7.3-20 grid_4.5.3 polyclip_1.10-7
[34] AlgDesign_1.2.1.2 colorspace_2.1-2 xtable_1.8-8
[37] scales_1.4.0 cli_3.6.6 rmarkdown_2.31
[40] crayon_1.5.3 generics_0.1.4 otel_0.2.0
[43] rstudioapi_0.18.0 tzdb_0.5.0 splines_4.5.3
[46] parallel_4.5.3 cellranger_1.1.0 base64enc_0.1-6
[49] vctrs_0.7.3 Matrix_1.7-4 sandwich_3.1-1
[52] jsonlite_2.0.0 hms_1.1.4 bit64_4.6.0-1
[55] visdat_0.6.0 htmlTable_2.5.0 Formula_1.2-5
[58] rematch_2.0.0 glue_1.8.1 codetools_0.2-20
[61] stringi_1.8.7 gtable_0.3.6 pillar_1.11.1
[64] htmltools_0.5.9 R6_2.6.1 vroom_1.7.1
[67] evaluate_1.0.5 lattice_0.22-9 backports_1.5.1
[70] broom_1.0.12 checkmate_2.3.4 gridExtra_2.3
[73] coda_0.19-4.1 nlme_3.1-168 mgcv_1.9-4
[76] xfun_0.57 zoo_1.8-15 pkgconfig_2.0.3