# Chapter 10 Linear models with a single, categorical *X* (“t-tests” and “ANOVA”)

## 10.1 A linear model with a single, categorical *X* variable estimates the effects of the levels of *X* on the response

To introduce a linear model with a single, categorical \(X\) variable, I’ll use data from a set of experiments designed to measure the effect of the lipid 12,13-diHOME on brown adipose tissue (BAT) thermoregulation and the mechanism of this effect.

Download the source data files and move to a new folder named “The cold-induced lipokine 12,13-diHOME promotes fatty acid transport into brown adipose tissue”.

Cold temperature and the neurotransmitter/hormone norepinephrine are known to stimulate increased thermogenesis in BAT cells. In this project, the researchers probed the question “what is the pathway that mediates the effect of cold-exposure on BAT thermogenesis?”. In the “discovery” component of this project, the researchers measured plasma levels of 88 lipids with known signaling properties in humans exposed to one hour of both normal (20 °C) and cold temperature (14 °C) temperature. Of the 88 lipids, 12,13-diHOME had the largest response to the cold treatment. The researchers followed this up with experiments on mice.

### 10.1.1 Example 1 (fig3d) – two treatment levels (“groups”)

Let’s start with the experiment in Figure 3d, which was designed to measure the effect of 12,13-diHOME on plasma triglyceride level. If 12,13-diHOME stimulates BAT activity, then levels in the 12,13-diHOME mice should be less than levels in the control mice.

#### 10.1.1.1 Step 1 – Understand the experiment design

**Response variable**: \(\texttt{serum_tg}\), A measure of serum triglycerides (mg/dl). \(\texttt{serum_tg}\) is a continuous variable.

**Factor variable**: \(\texttt{treatment}\), with levels: “Vehicle”, “12,13-diHOME” (the control or “Vehicle” mice were injected with saline). Coded as a factor.

**Design**: single, categorical *X*

#### 10.1.1.2 Step 2 – import

Open the data and, if necessary, wrangle into an analyzable format. The script to import these data is in the section Hidden code below.

#### 10.1.1.3 Step 3 – inspect the data

The second step is to examine the data to

- get a sense of sample size and balance
- check for biologically implausible outliers that suggest measurement failure, or transcription error (from a notebook, not in a cell)
- assess outliers for outlier strategy or robust analysis
- assess reasonable distributions and models for analysis.

`qplot(x = treatment, y = serum_tg, data = fig_3d)`

There are no obviously implausible data points. A normal distribution is a good, reasonable start. This can be checked more thoroughly after fitting the model.

#### 10.1.1.4 Step 4 – fit the model

`<- lm(serum_tg ~ treatment, data = fig_3d) fig3d_m1 `

#### 10.1.1.5 Step 5 – check the model

```
set.seed(1)
qqPlot(fig3d_m1, id=FALSE)
```

The Q-Q plot indicates the distribution of residuals is well within that expected for a normal sample and there is no cause for concern with inference.

`spreadLevelPlot(fig3d_m1, id=FALSE)`

```
##
## Suggested power transformation: 0.8167264
```

The spread-location plot shows no conspicuous trend in how the spread changes with the conditonal mean. There is no cause for concern with inference.

#### 10.1.1.6 Step 6 – inference

##### 10.1.1.6.1 coefficient table

```
<- cbind(coef(summary(fig3d_m1)),
fig3d_m1_coef confint(fig3d_m1))
kable(fig3d_m1_coef, digits = c(1,2,1,4,1,1)) %>%
kable_styling()
```

Estimate | Std. Error | t value | Pr(>|t|) | 2.5 % | 97.5 % | |
---|---|---|---|---|---|---|

(Intercept) | 42.6 | 1.67 | 25.6 | 0.0000 | 38.9 | 46.3 |

treatment12,13-diHOME | -7.2 | 2.36 | -3.0 | 0.0125 | -12.4 | -1.9 |

##### 10.1.1.6.2 emmeans table

```
<- emmeans(fig3d_m1, specs = "treatment")
fig3d_m1_emm kable(fig3d_m1_emm, digits = c(1,1,2,0,1,1)) %>%
kable_styling()
```

treatment | emmean | SE | df | lower.CL | upper.CL |
---|---|---|---|---|---|

Vehicle | 42.6 | 1.67 | 10 | 38.9 | 46.3 |

12,13-diHOME | 35.5 | 1.67 | 10 | 31.7 | 39.2 |

##### 10.1.1.6.3 contrasts table

```
<- contrast(fig3d_m1_emm,
fig3d_m1_pairs method = "revpairwise") %>%
summary(infer = TRUE)
kable(fig3d_m1_pairs, digits = c(1,1,1,0,1,1,2,4)) %>%
kable_styling()
```

contrast | estimate | SE | df | lower.CL | upper.CL | t.ratio | p.value |
---|---|---|---|---|---|---|---|

(12,13-diHOME) - Vehicle | -7.2 | 2.4 | 10 | -12.4 | -1.9 | -3.04 | 0.0125 |

#### 10.1.1.7 Step 6 – plot the model

```
ggplot_the_model(
fig3d_m1,
fig3d_m1_emm,
fig3d_m1_pairs,legend_position = "none",
y_label = "Serum TG (µg/dL)",
effect_label = "Effects (µg/dL)",
palette = pal_okabe_ito_blue,
rel_heights = c(0.5,1)
)
```

#### 10.1.1.8 Step 7 – report the model results

- Different ways of reporting the results without making claims that are not evidenced by the statistical analysis

“The estimated average effect of 12,13-diHOME on serum TG is 7.17 µg/dL (95% CI: -12.4, -1.9, \(p = 0.012\)).”

“Mean serum TG in mice with 12,13-diHOME (35.5 µg/dL, 95% CI: 31.7, 39.2) was 7.17 µg/dL less (95% CI: -12.4, -1.9, p = 0.012) than mean serum TG in control mice (42.6 µg/dL, 95% CI: 38.9, 46.3).”

- Problematic reports

“12,13-diHOME significantly reduced serum TG (\(p = 0.012\))” (1. In English, the current use of “significant” means “large” or “important”. Only knowledge of the physiological consequences of TG reduction over the range of the CI can be use as evidence of importance. 2. The statistics do not provide evidence that it is the serum TG that caused this reduction. While 7.17 µg/dL is the difference, some or all of the effect could be due to, in general, sampling, execution errors such as non-random treatment assignment, clustered caging, clustered application of assays, or non-blinded measurement of the response.

### 10.1.2 Understanding the analysis with two treatment levels

The variable \(\texttt{treatment}\) in the Figure 3d mouse experiment, is a single, categorical \(X\) variable. In a linear model, categorical variables are called **factors**. \(\texttt{treatment}\) can take two different values, “Vehicle” and “12,13-diHOME”. The different values in a factor are the **factor levels** (or just “levels”). “Levels” is a strange usage of this word; a less formal name for levels is “groups”. In a **Nominal** categorical factor, the levels have no units and are unordered, even if the variable is based on a numeric measurement. For example, I might design an experiment in which mice are randomly assigned to one of three treatments: one hour at 14 °C, one hour at 18 °C, or one hour at 26 °C. If I model this treatment as a nominal categorical factor, then I simply have three levels. While I would certainly choose to arrange these levels in a meaningful way in a plot, for the analysis itself, these levels have no units and there is no order. **Ordinal** categorical factors have levels that are ordered but there is no information on relative distance. The treatment at 18 °C is not more similar to 14 °C than to 26 °C. Nominal categorical factors is the default in R and how all factors are analyzed in this text.

#### 10.1.2.1 Linear models are regression models

The linear model fit to the serum TG data is

\[\begin{align} serum\_tg &= treatment + \varepsilon\\ \varepsilon &\sim N(0, \sigma^2) \tag{10.1} \end{align}\]

This notation is potentially confusing because the variable \(\texttt{treatment}\) is a factor containing the words “Vehicle” and “12,13-diHOME” and not numbers. The linear model in (10.1) can be specified using notation for a regression model using

\[\begin{align} serum\_tg &= \beta_0 + \beta_1 treatment_{12,13-diHOME} + \varepsilon\\ \varepsilon &\sim N(0, \sigma^2) \tag{10.2} \end{align}\]

Model (10.2) is a regression model where \(treatment_{12,13-diHOME}\) is not the variable \(\texttt{treatment}\), containing the words “Vehicle” or “12,13-diHOME” but a numeric variable that indicates membership in the group “12,13-diHOME”. This variable contains the number 1 if the mouse belongs to “12,13-diHOME” and the number 0 if the mouse doesn’t belong to “12,13-diHOME”. \(treatment_{12,13-diHOME}\) is known as an **indicator variable** because it indicates group membership. There are several ways of coding indicator variables and the way described here is called **dummy** or treatment coding. Dummy-coded indicator variables are sometimes called **dummy variables**.

The `lm`

function creates indicator variables under the table, in something called the **model matrix**.

```
<- model.matrix(~ treatment, data = fig_3d)
X <- nrow(X)
N 1:N,] X[
```

```
## (Intercept) treatment12,13-diHOME
## 1 1 0
## 2 1 0
## 3 1 0
## 4 1 0
## 5 1 0
## 6 1 0
## 7 1 1
## 8 1 1
## 9 1 1
## 10 1 1
## 11 1 1
## 12 1 1
```

The columns of the model matrix are the names of the **model terms** in the fit model. R names dummy variables by combining the names of the factor and the name of the level within the factor. So the \(X\) variable that R creates in the model matrix for the fit linear model in model (10.2) is \(treatment12,13-diHOME\). You can see these names as terms in the coefficient table of the fit model.

There are alternatives to dummy coding for creating indicator variables. Dummy coding is the default in R and it makes sense when thinking about experimental data with an obvious control level. I also like the interpretation of a “interaction effect” using Dummy coding. The classical coding for ANOVA is deviation effect coding, which creates coefficients that are deviations from the grand mean. In contrast to R, Deviation coding is the default in many statistical software packages including SAS, SPSS, and JMP. The method of coding can make a difference in an ANOVA table. Watch out for this – I’ve found several published papers where the researchers used the default dummy coding but interpreted the ANOVA table as if they had used deviation coding. This is both getting ahead of ourselves and somewhat moot, because I don’t advocate reporting ANOVA tables.

Recall from stats 101 that the slope of \(X\) in the model \(Y = b_0 + b_1 X\) is \(b_1 = \frac{\textrm{COV}(X,Y)}{\textrm{VAR}(X)}\). This can be generalized using the equation

\[\begin{equation} \mathbf{b} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y} \end{equation}\]

where \(\mathbf{X}\) is the model matrix containing a column for an intercept, columns for all indicator variables, and columns for all numeric covariates. \(\mathbf{b}\) is a vector containing the model coefficients, including the intercept in the first element. The first part of the RHS (\(\mathbf{X}^\top \mathbf{X}\)) is a matrix of the “sums of squares and cross-products” of the columns of \(\mathbf{X}\). Dividing each element of this matrix by \(N-1\) gives us the covariance matrix of the \(\mathbf{X}\), which contains the variances of the \(X\) columns along the diagonal, so this component has the role of the denominator in the stats 101 equation. Matrix algebra doesn’t do division, so the inverse of this matrix is multiplied by the second part. The second part or the RHS (\(\mathbf{X}^\top \mathbf{y}\)) is a vector containing the cross-products of each column of \(\mathbf{X}\) with \(\mathbf{y}\). Dividing each element of this vector by \(N-1\) gives us the covariances of each \(X\) with \(y\), so this component has the role of the numerator in the stats 101 equation.

Self-learning. `lm`

fits the model `y ~ X`

where X is the model matrix. Fit the model using the standard formula and the model using the model matrix. The coefficient table should be the same.

`m1 <- lm(serum_tg ~ treatment, data = fig_3d)`

`X <- model.matrix(~ treatment, data = fig_3d)`

`m2 <- lm(serum_tg ~ X, data = fig_3d)`

`coef(summary(m1))`

`coef(summary(m2))`

#### 10.1.2.2 The Estimates in the coefficient table are estimates of the parameters of the linear model fit to the data.

Estimate | Std. Error | t value | Pr(>|t|) | 2.5 % | 97.5 % | |
---|---|---|---|---|---|---|

(Intercept) | 42.6 | 1.67 | 25.6 | 0.0000 | 38.9 | 46.3 |

treatment12,13-diHOME | -7.2 | 2.36 | -3.0 | 0.0125 | -12.4 | -1.9 |

The row names of the coefficient table are the column names of the model matrix. These are the model terms. There are two terms (two rows) because there are two parameters in the regression model (10.2). The values in the column \(\texttt{Estimate}\) in the coefficient table are the estimates of the regression parameters \(\beta_0\) and \(\beta_1\). These estimates are the coefficients of the fit model, \(b_0\) and \(b_1\).

#### 10.1.2.3 The coefficients of a linear model using dummy coding have a useful interpretation

coefficient | parameter | model term | interpretation |
---|---|---|---|

\(b_0\) | \(\beta_0\) | (Intercept) | \(\overline{Vehicle}\) |

\(b_1\) | \(\beta_1\) | treatment12,13-diHOME | \(\overline{12,13\;diHOME} - \overline{Vehicle}\) |

It is important to understand the interpretation of the coefficients of the fit linear model (10.1) (Table 10.1).

- The coefficient \(b_0\) is the is the conditional mean of the response for the reference level, which is “Vehicle”. Remember that a conditional mean is the mean of a group that all have the same value for one or more \(X\) variables.
- The coefficient \(b_1\) is the difference between the conditional means of the “12,13-diHOME” level and the reference (“Vehicle”) level:

\[\begin{equation} \mathrm{E}[serum\_tg|treatment = \texttt{"12,13-diHOME"}] - \mathrm{E}[serum\_tg|treatment = \texttt{"Vehicle"}] \end{equation}\]

Because there are no additional covariates in model, this difference is equal to the difference between the sample means \(\bar{Y}_{12,13-diHOME} - \bar{Y}_{Vehicle}\). The *direction* of this difference is important – it is non-reference level minus the reference level.

The estimate \(b_1\) is the **effect** that we are interested in. Specifically, it is the measured effect of 12,13-diHOME on serum TG. When we inject 12,13-diHOME, we find the mean serum TG decreased by -7.2 µg/dL relative to the mean serum TG in the mice that were injected with saline. Importantly, the reference level is not a property of an experiment but is set by whomever is analyzing the data. Since the non-reference estimates are differences in means, it often makes sense to set the “control” treatment level as the reference level.

Many beginners mistakenly memorize the coefficient \(b_1\) to equal the mean of the non-reference group (“12,13-diHOME”). Don’t do this. In a regression model, only \(b_0\) is a mean. The coefficient \(b_1\) in model (10.2) is a difference in means.

A geometric interpretation of the coefficients is illustrated in Figure 10.1. \(b_0\) is the conditional mean of the reference level (“Vehicle”) and is an estimate of \(\beta_0\), the true, conditional mean of the population. \(b_1\) is the difference in the conditional means of the first non-reference level (“12,13-diHOME”) and the reference level (“Vehicle”) and is an estimate of \(\beta_1\), the true difference in the conditional means of the population with and without the treatment 12,13-diHOME.

tl;dr. What is a population? In the experimental biology examples in this text, we might consider the population as a very idealized, infinitely large set of mice, or fish, or fruit flies, or communities from which our sample is a reasonably representative subset. For the experiments in the 12,13-diHOME study, the population might be conceived of as the hypothetical, infinitely large set of 12-week-old, male, C57BL/6J mice, raised in the mouse facility at Joslin Diabetes Center. An even more abstract way to way to think about what the population could be is the infinitely large set of values that could generated by the linear model.

#### 10.1.2.4 Better know the coefficient table

```
%>%
fig3d_m1_coef kable(digits = c(1,2,1,4,1,1)) %>%
kable_styling()
```

Estimate | Std. Error | t value | Pr(>|t|) | 2.5 % | 97.5 % | |
---|---|---|---|---|---|---|

(Intercept) | 42.6 | 1.67 | 25.6 | 0.0000 | 38.9 | 46.3 |

treatment12,13-diHOME | -7.2 | 2.36 | -3.0 | 0.0125 | -12.4 | -1.9 |

- The \(\texttt{(Intercept)}\) row contains the statistics for \(b_0\) (the estimate of \(\beta_0\)). Remember that \(b_0\) is the conditional mean of the reference treatment (“Vehicle”).
- The \(\texttt{treatment12,13-diHOME}\) row contains the statistics for \(b_1\) (the estimate of \(\beta_1\)). Remember that \(b_1\) is the difference in conditional means of the groups “12,13-diHOME” and “Vehicle”.
- The column \(\texttt{Estimate}\) contains the model coefficients, which are estimates of the parameters.
- The column \(\texttt{Std. Error}\) contains the model SEs of the coefficients. The SE of \(\texttt{(Intercept)}\) is a standard error of a mean (SEM). The SE of \(\texttt{treatment12,13-diHOME}\) is a standard error of a difference (SED).
- The column \(\texttt{t value}\) contains the test statistic of the coefficients. This value is the ratio \(\frac{Estimate}{SE}\). For this model, we are only interested in the test statistic for \(b_1\). Effectively, we will never be interested in the test statistic for \(b_0\) because the mean of a group will never be zero.
- The column \(\texttt{Pr(>|t|)}\) contains the
*p*-values for the test statistic of the coefficients. For this model, and all models in this text, we are only interested in the*p*-value for the non-intercept coefficients. - The columns \(\texttt{2.5 %}\) and \(\texttt{97.5 %}\) contain the lower and upper limits of the 95% confidence interval of the estimate.

#### 10.1.2.5 The emmeans table is a table of modeled means and inferential statistics

treatment | emmean | SE | df | lower.CL | upper.CL |
---|---|---|---|---|---|

Vehicle | 42.6 | 1.67 | 10 | 38.9 | 46.3 |

12,13-diHOME | 35.5 | 1.67 | 10 | 31.7 | 39.2 |

The table of marginal means for the model fit to the fig3d serum TG data (Model (10.1)) is given in Table @ref(tab:fig3d_m1_emm_kable). The table of marginal means gives the **modeled* mean, standard error and confidence interval for all specified groups. There is no test-statistic with a *p*-value because there is no significance test. In this text, I’ll refer to this table as the “emmeans table”, since it is the output from the `emmeans`

function (“em” is the abbreviation for “estimated marginal”). I’ll use “modeled means” to refer to the means themselves as these are the estimate of means from a fit linear model.

A **marginal mean** is the mean over a set of conditional means. For example, if a treatment factor has three levels, the conditional means are the means for each level and the marginal mean is the mean of the three means. Or, if the conditional means are the expected values given a continous covariate, the marginal mean is the expected value at the mean of covariate. The specified emmeans table of the fig3d data is not too exciting because it simply contains the conditional means – the values are not marginalized over any \(X\). Because the emmeans table contains different sorts of means (conditional, marginal, adjusted), this text will generally refer to the means in this table as “modeled means”.

treatment | emmean | SE | df | lower.CL | upper.CL |
---|---|---|---|---|---|

Vehicle | 42.62004 | 1.66723 | 10 | 38.90523 | 46.33485 |

12,13-diHOME | 35.45233 | 1.66723 | 10 | 31.73752 | 39.16714 |

treatment | mean | SE | lower.CL | upper.CL |
---|---|---|---|---|

Vehicle | 42.62004 | 1.77325 | 38.66899 | 46.57109 |

12,13-diHOME | 35.45233 | 1.55398 | 31.98984 | 38.91482 |

- The means in the emmeans table are modeled means. Here, and for many linear models, these will be equal to the sampled means. This will not be the case in models with one or more continuous covariates, or if marginal means are explicitly specified.
- Unlike the modeled means, the modeled standard errors and confidence intervals will, effectively, never equal sample standard errors and confidence intervals. In many models, plots using sample statistics can lead to very deceiving inference of differences between groups. This is why this text advocates plotting the model – using modeled means and confidence intervals in plots.

It is exceptionally important to understand the difference between the means, SEs, and CIs in the emmeans table and the statistics of the same name in a summary table of the data.

- The statistics in a summary table are
**sampled**means, SEs, and CIs – these statistics are computed for the group using*only the data in the group*. - To understand modeled SEs and CIs, recall that the standard error of a sample mean is \(\frac{s}{\sqrt{n}}\), where \(s\) is the sample standard deviation and \(n\) is the sample size in the group. The computation of the SE in the emmeans table uses the same equation, except the numerator is not the sample standard deviation of the group but the model standard deviation \(\hat{\sigma}\). As with the sample SE, the denominator for the modeled SE is the sample size \(n\) for the group. Since the numerator of the modeled SE is the same for all groups, the modeled SE will be the same in all groups that have the same sample size, as seen in the marginal means table for the model fit to the Figure 3d data. This is not true for sampled SEs, since sampled standard deviations will always differ.

It may seem odd to use a common standard deviation in the computation of the modeled SEs. It is not. Remember that an assumption of the linear model is homogeneity of variances – that all residuals \(e_i\) are drawn from the same distribution (\(N(0, \sigma^2)\)) (a “single hat”) regardless of group. The model standard deviation \(\hat{\sigma}\) is the estimate of the square root of the variance of this distribution. Given this interpretation, it is useful to think of each sample standard deviation as an estimate of \(\sigma\) (the linear model assumes that all differences among the sample standard deviations are due entirely to sampling). The model standard deviation is a more precise estimate of \(\sigma\) since it is computed from a larger sample (all \(N\) residuals).

The model standard deviation is called the “pooled” standard deviation in the ANOVA literature and is computed as a sample-size weighted average of the sample standard deviations.

The modeled standard error of the mean uses the estimate of \(\sigma\) from the fit model. This estimate is

\[\begin{equation} \hat{\sigma} = \sqrt{\frac{\sum{(y_i - \hat{y}_i)^2}}{df}} \end{equation}\]

Create a code chunk that computes this. Recall that \((y_i - \hat{y}_i)\) is the set of residuals from the model, which can be extracted using `residuals(fit)`

where “fit” is the fit model object. \(df\) is the model degrees of freedom, which is \(N-k\), where \(N\) is the total sample size and \(k\) is the number of parameters that are fit. This makes sense – for the sample variance there is one parameter that is fit, the mean of the group. In model fig3d_m1, there are two parameters that are fit, the intercept and the coefficient of treatment12,13-diHOME.

#### 10.1.2.6 Estimates of the effects are in the contrasts table

contrast | estimate | SE | df | lower.CL | upper.CL | t.ratio | p.value |
---|---|---|---|---|---|---|---|

(12,13-diHOME) - Vehicle | -7.2 | 2.36 | 10 | -12.4 | -1.9 | -3 | 0.012 |

- This table is important for reporting treatment effects and CIs and for plotting the model. A contrast is a difference in means. With only two treatment levels, the table of contrasts doesn’t give any more information than the coefficient table – the single contrast is the coefficient \(b_1\) in the coefficient table. Nevertheless, I advocate computing this table to stay consistent and because the script (or function) to plot the model uses this table and not the coefficient table.
- The value in the column \(\texttt{Estimate}\) is the mean of the non-reference group (“12,13-diHOME”) minus the mean of the reference group (“Vehicle”).
- The value in the “SE” column is the standard error of the difference (SED), specifically the difference in the estimate column. This SE is computed using the model standard deviation \(\sigma\).
- The values in the “lower.CL” and “upper.CL” columns are the bounds of the 95% confidence interval of the estimate. Remember (from Chapter
**??**) to think of this interval as containing potential values of the true parameter (the true difference in means between the two groups) that are reasonably compatible with the data. Don’t think of the interval as having 95% probability of containing the true effect. Remember that a confidence interval applies to the procedure and not a parameter – 95% of the CIs from hypothetical, replicate experiments that meet all the assumptions used to compute the CI will include the true effect. - The columns “t.ratio” and “p.value” contains the
*t*and*p*values of the significance (not hypothesis!) test of the estimate. The*t*-statistic is the ratio of the estimate to the SE of the estimate (use the console to confirm this given the values in the table). It is a signal (the estimate) to noise (SE of the estimate) ratio. The*p*-value is the probability of sampling from normal distribution with the observed standard deviation, randomly assigning the sampled values to either “Vehicle” or “12,13-diHOME”, fitting the linear model, and observing a*t*-value as or more extreme than the observed*t*. A very small*p*-value is consistent with the experiment “not sampling from distributions with the same mean” – meaning that adding a treatment affects the mean of the distribution. This is the logic used to infer a treatment effect. Unfortunately, it is also consistent with the experiment not approximating other conditions of the model, including non-random assignment, non-independence, non-normal conditional responses, and variance heterogeneity. It is up to the rigorous researcher to be sure that these other model conditions are approximated or “good enough” to use the*p*-value to infer a treatment effect on the mean.

#### 10.1.2.7 *t* and *p* from the contrasts table – when there are only two levels in \(X\) – are the same as *t* and *p* from a *t*-test

Compare

coefficient table:

```
<- lm(serum_tg ~ treatment, data = fig_3d)
m1 coef(summary(m1)) %>%
kable() %>%
kable_styling()
```

Estimate | Std. Error | t value | Pr(>|t|) | |
---|---|---|---|---|

(Intercept) | 42.620042 | 1.667226 | 25.563447 | 0.000000 |

treatment12,13-diHOME | -7.167711 | 2.357813 | -3.039982 | 0.012463 |

contrast table:

```
emmeans(m1, specs = "treatment") %>%
contrast(method = "revpairwise") %>%
kable() %>%
kable_styling()
```

contrast | estimate | SE | df | t.ratio | p.value |
---|---|---|---|---|---|

(12,13-diHOME) - Vehicle | -7.167711 | 2.357813 | 10 | -3.039982 | 0.012463 |

t-test:

```
<- t.test(fig_3d[treatment == "12,13-diHOME", serum_tg],
m2 == "Vehicle", serum_tg],
fig_3d[treatment var.equal = TRUE)
glance(m2) %>% # glance is from the broom package
kable() %>%
kable_styling()
```

estimate | estimate1 | estimate2 | statistic | p.value | parameter | conf.low | conf.high | method | alternative |
---|---|---|---|---|---|---|---|---|---|

-7.167711 | 35.45233 | 42.62004 | -3.039982 | 0.012463 | 10 | -12.42125 | -1.914175 | Two Sample t-test | two.sided |

Notes

- The default
`t.test`

in R is the Welch t-test for heterogenous variance. To compute the Student t-test, use`var.equal = TRUE`

. - The “statistic” in the t-test output contains the
*t*-value of the*t*-test. It is precisely the same as the*t*-statistic in the coefficient table and the contrast table. - The
*p*-values in all three tables are precisely the same.

The *t* and *p* values for the *t*-test are the same as those for the linear model, because the *t*-test is a specific case of the linear model. Reasons to abandon classic *t*-tests and learn the linear modeling strategy include

- A linear modeling strategy encourages researchers to think about the effect and uncertainty in the effect and not just a
*p*-value. - The linear model is nearly infinitely flexible and expandible while the
*t*-test is limited to a few variations.

There is rarely a reason to ever use the `t.test()`

function. Throw the function away. Ignore web pages that teach you to use it. The *t*-test is easy to learn, which encourages its overuse. If your only tool is a t-test, every problem looks like a comparison between two-means.

### 10.1.3 Example 2 – three treatment levels (“groups”)

#### 10.1.3.1 Understand the experiment design

The data come from the experiment reported in Figure 2a of the 12,13-diHOME article described above. This experiment was designed to probe the hypothesis that 12,13-diHOME is a mediator of known stimulators of increased BAT activity (exposure to cold temperature and sympathetic nervous system activation). Mice were assigned to control (30 °C), one-hour exposure to 4 °C, or 30 minute norepinephrine (NE) treatment level (NE is the neurotransmitter of the sympathetic neurons targeting peripheral tissues).

**response variable**: \(\texttt{diHOME}\), the serum concentration of 12,13-diHOME. a continuous variable.

**treatment variable**: \(\texttt{treatment}\), with levels: “Control”, “1 hour cold”, “30 min NE”. Coded as a factor.

- “Control” – the negative control. We expect diHOME to be low relative to the two treated conditions.
- “1 hour cold” – focal treatment 1. If diHOME is a mediator, then response should be relatively high.
- “30 min NE” – focal treatment 2. If diHOME is a mediator, then response should be relatively high.

**planned comparisons**

- “1 hour cold” - “Control” – If diHOME is a mediator of cold, then difference should be positive.
- “30 min NE” - “Control” – If diHOME is a mediator of NE, then difference should be positive.

**design**: single, categorical *X*

#### 10.1.3.2 fit the model

`<- lm(diHOME ~ treatment, data = fig2a) fig2a_m1 `

#### 10.1.3.3 check the model

```
set.seed(1)
qqPlot(fig2a_m1)
```

`## [1] 3 18`

The Q-Q plot indicates potential issues at the extreme quantiles, what is called “heavy tails”. The two values are the extreme values in the “30 min NE” group. This could be the result of a small sample from a response with a larger variance.

`spreadLevelPlot(fig2a_m1, id=FALSE)`

```
##
## Suggested power transformation: 0.4430799
```

The combination of the raw residuals and the spread-level plot suggests heterogeneity but low confidence in anything given the small sample size.

#### 10.1.3.4 Inference from the model

##### 10.1.3.4.1 coefficient table

```
<- cbind(coef(summary(fig2a_m1)),
fig2a_m1_coef confint(fig2a_m1))
%>%
fig2a_m1_coef kable(digits = c(1,2,1,4,1,1)) %>%
kable_styling()
```

Estimate | Std. Error | t value | Pr(>|t|) | 2.5 % | 97.5 % | |
---|---|---|---|---|---|---|

(Intercept) | 12.0 | 3.08 | 3.9 | 0.0016 | 5.4 | 18.6 |

treatment1 hour cold | 7.1 | 4.57 | 1.6 | 0.1405 | -2.7 | 16.9 |

treatment30 min NE | 14.8 | 4.36 | 3.4 | 0.0044 | 5.4 | 24.1 |

##### 10.1.3.4.2 emmeans table

```
<- emmeans(fig2a_m1, specs = "treatment")
fig2a_m1_emm
%>%
fig2a_m1_emm kable(digits = c(1,1,2,0,1,1)) %>%
kable_styling()
```

treatment | emmean | SE | df | lower.CL | upper.CL |
---|---|---|---|---|---|

Control | 12.0 | 3.08 | 14 | 5.4 | 18.6 |

1 hour cold | 19.2 | 3.38 | 14 | 11.9 | 26.4 |

30 min NE | 26.8 | 3.08 | 14 | 20.2 | 33.4 |

##### 10.1.3.4.3 contrasts table

```
<- contrast(fig2a_m1_emm,
fig2a_m1_planned method = "trt.vs.ctrl",
adjust = "none",
level = 0.95) %>%
summary(infer = TRUE)
%>%
fig2a_m1_planned kable(digits = c(1,1,1,0,1,1,2,4)) %>%
kable_styling()
```

contrast | estimate | SE | df | lower.CL | upper.CL | t.ratio | p.value |
---|---|---|---|---|---|---|---|

1 hour cold - Control | 7.1 | 4.6 | 14 | -2.7 | 16.9 | 1.56 | 0.1405 |

30 min NE - Control | 14.8 | 4.4 | 14 | 5.4 | 24.1 | 3.40 | 0.0044 |

#### 10.1.3.5 plot the model

```
ggplot_the_model(
fig2a_m1,
fig2a_m1_emm,
fig2a_m1_planned,legend_position = "none",
y_label = "12,13-diHOME (pmol/mL)",
effect_label = "Effects (pmol/mL)",
palette = pal_okabe_ito_blue,
rel_heights = c(0.5,1)
)
```

#### 10.1.3.6 Report the model results

Compared to control mice (12.0 pmol/mL, 95% CI: 5.4, 18.6), mean serum 12,13-diHOME in mice exposed to one-hour cold (19.2 pmol/mL, 95% CI: 11.9, 26.4) was 7.1 pmol/mL higher (95% CI: -2.7, 16.9, \(p = 0.14\)) while mean Serum 12,13-diHOME in mice exposed to 30 minutes NE (26.8 pmol/mL, 95% CI: 20.2, 33.3) was 14.8 pmol/mL higher (95% CI: 5.4, 24.1 \(p = 0.004\)).

### 10.1.4 Understanding the analysis with three (or more) treatment levels

#### 10.1.4.1 Better know the coefficient table

The fit regression model for the data in Figure 2a is

\[\begin{equation} diHOME_i = b_0 + b_1 treatment_{1\_hour\_cold,i} + b_2 treatment_{30\_min\_NE,i} + e_i \tag{10.3} \end{equation}\]

The coefficients of the model are in the \(\texttt{Estimate}\) column of the coefficient table.

Estimate | Std. Error | t value | Pr(>|t|) | 2.5 % | 97.5 % | |
---|---|---|---|---|---|---|

(Intercept) | 12.0 | 3.1 | 3.9 | 0 | 5.4 | 18.6 |

treatment1 hour cold | 7.1 | 4.6 | 1.6 | 0 | -2.7 | 16.9 |

treatment30 min NE | 14.8 | 4.4 | 3.4 | 0 | 5.4 | 24.1 |

- The \(\texttt{(Intercept)}\) row contains the statistics for \(b_0\) (the estimate of \(\beta_0\)). Here, \(b_0\) is the mean of the reference group, which is “Control”.
- The \(\texttt{treatment1 hour cold}\) row contains the statistics for \(b_1\) (the estimate of \(\beta_1\)). Here, \(b_1\) is the difference \(\mathrm{E}[diHOME|treatment = \texttt{"1 hour cold"}] - \mathrm{E}[diHOME|treatment = \texttt{"Control"}]\). This difference in conditional means is equal to the difference in the sample means of the two groups
*for this model*because there are no additional covariates in the model. - The \(\texttt{treatment30 min NE}\) row contains the statistics for \(b_2\) (the estimate of \(\beta_2\)). Here, \(b_2\) is the difference \(\mathrm{E}[diHOME|treatment = \texttt{"30 min NE"}] - \mathrm{E}[diHOME|treatment = \texttt{"Control"}]\). Do not make the mistake in thinking that the value in \(\texttt{Estimate}\) is the mean of the “30 min NE” group.
- The number of non-intercept coefficients generalizes to any number of levels of the factor variable. If there are \(k\) levels of the factor, there are \(k-1\) indicator variables, each with its own coefficient (\(b_1\) through \(b_{k-1}\)) estimating the effect of that treatment level relative to the control (if using dummy coding).

Again – Do not make the mistake in thinking that the values in \(\texttt{Estimate}\) for the \(\texttt{treatment1 hour cold}\) and \(\texttt{treatment30 min NE}\) rows are the *means* of the “1 hour cold” and “30 min NE” groups. These coefficients are differences in means. And, to emphasize further understanding of these coefficients, both \(b_1\) and \(b_2\) are “slopes”. Don’t visualize this as a single line from the control mean through both non-control means. as two lines. Slopes is plural – there are two regression lines. \(b_1\) is the slope of the line from the control mean to the “1 hour cold” mean. \(b_2\) is the slope of the line from the control mean to the “30 min NE” mean. The numerator of each slope is the difference between that group’s mean and the control mean. The denominator of each slope is 1 (because each has the value 1 when the row is assigned to that group).

Two understand the names of the model terms, it’s useful to recall the order of the factor levels of \(\texttt{treatment}\), which is

`levels(fig2a$treatment) `

`## [1] "Control" "1 hour cold" "30 min NE"`

Given this ordering, the `lm`

function creates a regression model with an intercept column for the “Control” group (because it is first in the list), an indicator variable for the “1 hour cold” group called `treatment1 hour cold`

, and an indicator variable for the “30 min NE” group called `treatment30 min NE`

. We can see these model names by peeking at the model matrix of the fit model

```
<- model.matrix(fig2a_m1)
fig2a_m1_X
head(fig2a_m1_X)
```

```
## (Intercept) treatment1 hour cold treatment30 min NE
## 1 1 0 0
## 2 1 1 0
## 3 1 0 1
## 4 1 0 0
## 5 1 1 0
## 6 1 0 1
```

The column \(\texttt{treatment1 hour cold}\) is a dummy-coded indicator variable containing the number 1, if the individual is in the “1 hour cold” group, or the number 0, otherwise. The column \(\texttt{treatment30 min NE}\) is a dummy-coded indicator variable containing the number 1, if the individual is in the “30 min NE” group, or the number 0, otherwise.

The model coefficients, parameters, model term, and interpretation are summarized in the following table.

coefficient | parameter | model term | interpretation |
---|---|---|---|

\(b_0\) | \(\beta_0\) | (Intercept) | \(\overline{Control}\) |

\(b_1\) | \(\beta_1\) | treatment1 hour cold | \(\overline{1\;hour\;cold} - \overline{Control}\) |

\(b_2\) | \(\beta_2\) | treatment30 min NE | \(\overline{30\;min\;NE} - \overline{Control}\) |

#### 10.1.4.2 The emmeans table

treatment | emmean | SE | df | lower.CL | upper.CL |
---|---|---|---|---|---|

Control | 12.0 | 3.08 | 14 | 5.4 | 18.6 |

1 hour cold | 19.2 | 3.38 | 14 | 11.9 | 26.4 |

30 min NE | 26.8 | 3.08 | 14 | 20.2 | 33.4 |

This table is important for reporting means and CIs and for plotting the model. As in example 1, the values in the column “emmean” are the sample means of each group (what you would compute if you simply computed the mean for that group). Again, this is true for this model, but is not generally true.

And again, as in example 1, the SE for each mean is *not* the sample SE but the modeled SE – it is computed from the residuals of the model and is an estimate of \(\sigma\). *These are the SEs that you should report* because it is these SEs that are used to compute the *p*-value and CI that you report, that is, they tell the same “story”. The SE for the “1 hour cold” group is a bit higher because the sample size \(n\) for this group is smaller by 1.

#### 10.1.4.3 The contrasts table

contrast | estimate | SE | df | lower.CL | upper.CL | t.ratio | p.value |
---|---|---|---|---|---|---|---|

1 hour cold - Control | 7.1 | 4.57 | 14 | -2.7 | 16.9 | 1.6 | 0.1405 |

30 min NE - Control | 14.8 | 4.36 | 14 | 5.4 | 24.1 | 3.4 | 0.0044 |

- This table is important for reporting treatment effects and CIs and for plotting the model. A contrast is a difference in means.
- The contrast table here has no more information than is in the coefficient table, but that is not generally true for models with treatment factors with more than two groups. In “Working in R” below, I show how to compute a contrast table with all pairwise comparisons (contrasts between all possible pairings of the groups)
- The column \(\texttt{Contrast}\) contains the names of the contrasts. Note that the name gives the direction of the difference.
- The values in the column \(\texttt{estimate}\) are the contrasts. These are the differences in the conditional means of the groups identified in the \(\texttt{Contrast}\) column. These are the
**effects**that we are interested in. - The value in the “SE” column is the standard error of the difference (SED) of each contrast. This SE is computed using the model standard deviation \(\sigma\).
- The values in the “lower.CL” and “upper.CL” columns are the bounds of the 95% confidence interval of the estimate. Remember (from Chapter
**??**) to think of this interval as containing potential values of the true parameter (the true difference in means between the two groups) that are reasonably compatible with the data. Don’t think of the interval as having 95% probability of containing the true effect. Remember that a confidence interval applies to the procedure and not a parameter – 95% of the CIs from hypothetical, replicate experiments that meet all the assumptions used to compute the CI will include the true effect. - The columns “t.ratio” and “p.value” contains the
*t*and*p*values of the significance (not hypothesis!) test of the estimate. The*t*-statistic is the ratio of the estimate to the SE of the estimate (use the console to confirm this given the values in the table). It is a signal (the estimate) to noise (SE of the estimate) ratio. The*p*-value is the probability of sampling from normal distribution with the observed standard deviation, randomly assigning the sampled values to the three groups (using the original sample sizes for each), fitting the linear model, and observing a*t*-value as or more extreme than the observed*t*. A very small*p*-value is consistent with the experiment “not sampling from distributions with the same mean” – meaning that adding a treatment affects the mean of the distribution. This is the logic used to infer a treatment effect. Unfortunately, it is also consistent with the experiment not approximating other conditions of the model, including non-random assignment, non-independence, non-normal conditional responses, and variance heterogeneity. It is up to the rigorous researcher to be sure that these other model conditions are approximated or “good enough” to use the*p*-value to infer a treatment effect on the mean.

#### 10.1.4.4 *t* and *p* from the contrasts table – when there are more than two levels in \(X\) – are not the same as those from pairwise *t*-tests among pairs of groups

```
<- contrast(fig2a_m1_emm,
fig2a_m1_pairs method = "revpairwise",
adjust = "none") %>%
summary(infer = TRUE)
```

The chunk above computes a contrast table that includes comparisons of all pairs of groups in the factor \(\texttt{treatment}\) (this adds a 3rd comparison to the contrast table of planned comparisons above). The *t*-tests for the contrasts are derived from a single fit linear model.

In contrast to the analysis in the chunk above, researchers commonly fit separate *t*-tests for each pair of treatment levels.

```
# classic t-test
<- t.test(fig2a[treatment == "1 hour cold", diHOME],
test1 == "Control", diHOME],
fig2a[treatment var.equal = TRUE)
<- t.test(fig2a[treatment == "30 min NE", diHOME],
test2 == "Control", diHOME],
fig2a[treatment var.equal = TRUE)
<- t.test(fig2a[treatment == "30 min NE", diHOME],
test3 == "1 hour cold", diHOME],
fig2a[treatment var.equal = TRUE)
```

Notes

- Again, the default t.test in R is the Welch t-test for heterogenous variance. To compute the Student t-test, use
`var.equal = TRUE`

- To see the full
`t.test`

output, type “test1” into the console.

Compare the *t* and *p* values from the three independent tests with the *t* and *p*-values from the single linear model.

contrast | t (lm) | p (lm) | t (t-test) | p (t-test) |
---|---|---|---|---|

1 hour cold - Control | 1.562324 | 0.1405278 | 2.415122 | 0.0389207 |

30 min NE - 1 hour cold | 1.674696 | 0.1161773 | 1.380666 | 0.2007006 |

30 min NE - Control | 3.395015 | 0.0043559 | 3.238158 | 0.0088971 |

The *t* and *p*-values computed from three separate tests differ from the *t* and *p*-values computed from the single linear model shown in the contrasts table above. The values differ because the SE in the denominators used to compute the \(t\)-values differ. The linear model uses the same value of \(\sigma\) to compute the SED (the denominator of *t*) for all three t-tests in the contrast table. Each separate t-test uses a different value of \(\sigma\) to compute the SED. Which is correct? Neither – they simply make different assumptions about the data generating model.

- Most importantly,
**never**do both methods, look at the*p*-values, and then convince yourself that the method with the*p*-values that match your hypothesis is the correct method. Human brains are very, very good at doing this. This is called p-hacking. When you p-hack, the interpretation of the*p*-value is effectively meaningless. P-hacking leads to irreproducible science. - In general, using the linear model is a better practice than the separate
*t*-tests. The reason is the homogeneity of variance assumption. If we assume homogeneity of variances, then we can think of the sample standard deviation of all three groups as an estimate of \(\sigma\). In the linear model, we use three groups to estimate \(\sigma\) but in each separate*t*-test, we use only two groups. Consequently, our estimate of \(\sigma\) in the linear model is more precise than that in the*t*-tests. While the difference can be large with any individual data set (it’s pretty big with the fig2a data), the long-run advantage of using the linear model instead of separate*t*-tests is pretty small, especially with only three groups (the precision increases with more groups). - We can drop the homogeneity of variance assumption with either the linear model or the three, separate
*t*-tests. This is outlined below in “Heterogeneity of variance”. In this case, the*t*and*p*-values for the three comparisons are the same. Still, the linear model (that models heterogenity) is better practice than the separate*t*-tests because the linear model is much more flexible and expandable.

#### 10.1.4.5 The contrasts table – when there are more than two levels in \(X\) – has multiple p-values. How to handle this “multiple testing” is highly controversial

**Multiple testing** is the practice of adjusting *p*-values and confidence intervals to account for the expected increase in the frequency of Type I error in a batch, or *family*, of tests. Multiple testing is a concept that exists because of Neyman-Pearson hypothesis testing strategy. If multiple tests are used to answer the same question then, these are in the same family. Issues surrounding multiple testing are fleshed out in more detail in Chapter xxx “Best Practices”. Computing adjusted values is covered below in the “Working in R” section.

## 10.2 Working in R

### 10.2.1 Fit the model

`<- lm(diHOME ~ treatment, data = fig2a) m1 `

- As described in Fitting the linear model,
`diHOME ~ treatment`

is the model formula. All functions that will be used to fit models in this text will use a model formula. Many base R and package functions that are meant to be either “easy” or follow an ANOVA strategy use lists of dependent and independent variables instead of a model formula. We won’t use these because they are not consistent with a linear modeling way of thinking about data analysis. - \(\texttt{treatment}\) was specifically coded as a factor variable in the import and wrangle chunk and R will automatically create the correct indicator variables. If categorical variables on the RHS of the formula have not be converted to factors by the user, then R will treat character variables as factors and create the indicator variables. If the categorical variable is numeric (for example, the variable
`time`

might have levels 1, 2, 3 instead of “1 hour”, “2 hours”, “three hours”), R will treat the variable as numeric and not create the indicator variables. A user could use the formula`y ~ factor(time)`

to force R to create the indicator variable. I prefer to explicitly create a version of the variable that is recoded as factor using something like`dt[, time_fac := factor(time, levels = c("1", "2", "3"))]`

.

### 10.2.2 Controlling the output in tables using the coefficient table as an example

```
<- cbind(coef(summary(m1)),
m1_coef confint(m1))
m1_coef
```

```
## Estimate Std. Error t value Pr(>|t|) 2.5 %
## (Intercept) 12.023075 3.081337 3.901902 0.001595771 5.414264
## treatment1 hour cold 7.140386 4.570362 1.562324 0.140527829 -2.662066
## treatment30 min NE 14.794354 4.357669 3.395015 0.004355868 5.448083
## 97.5 %
## (Intercept) 18.63189
## treatment1 hour cold 16.94284
## treatment30 min NE 24.14063
```

The standard print of the `m1_coef`

object is too long for the display and often wraps when printed on the display, which can make it hard to read. I use `knitr::kable`

to print the table with fewer decimal places and `kableExtra::kable_styling`

to make it a little prettier.

```
# the row names are not part of the m1_coef object
# so there is no digit designation for this column
%>% # pipe the m1_coef object to kable
m1_coef kable(digits = c(2,3,3,5,2,2)) %>%
kable_styling()
```

Estimate | Std. Error | t value | Pr(>|t|) | 2.5 % | 97.5 % | |
---|---|---|---|---|---|---|

(Intercept) | 12.02 | 3.081 | 3.902 | 0.00160 | 5.41 | 18.63 |

treatment1 hour cold | 7.14 | 4.570 | 1.562 | 0.14053 | -2.66 | 16.94 |

treatment30 min NE | 14.79 | 4.358 | 3.395 | 0.00436 | 5.45 | 24.14 |

`# explore other styles in the kableExtra package`

### 10.2.3 Using the emmeans function

```
<- emmeans(m1, specs = "treatment")
m1_emm m1_emm
```

```
## treatment emmean SE df lower.CL upper.CL
## Control 12.0 3.08 14 5.41 18.6
## 1 hour cold 19.2 3.38 14 11.92 26.4
## 30 min NE 26.8 3.08 14 20.21 33.4
##
## Confidence level used: 0.95
```

Notes

- Note that printing the
`emmeans`

object displays useful information. Here, this information includes the confidence level used. If the object is printed using`kable`

(as in the “Inference” and “Understanding” sections above), only the table is printed and the additional information is lost. `emmeans`

computes the modeled means of all combinations of the levels of the factor variables specified in`specs`

.- If there are two factor variables in the model, and both are passed to
`specs`

, then the modeled means of all combinations of the levels of the two variables are computed. If only one factor variable is passed, then the marginal means (averaged over all levels of the missing factor) are computed. This will become more clear in the chapter “Models for two (or more) categorical X variables”. - If there are continuous covariates in the model, the modeled means are computed at the average values of these covariates. These covariates do not need to be passed to the
`specs`

argument. - You can pass numeric and integer covariates to
`specs`

to control the value of the covariates used compute the modeled means. This is outlined in Adding covariates to a linear model

### 10.2.4 Using the contrast function

```
<- contrast(m1_emm,
m1_planned method = "trt.vs.ctrl",
adjust = "none",
level = 0.95) %>%
summary(infer = TRUE)
m1_planned
```

```
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## 1 hour cold - Control 7.14 4.57 14 -2.66 16.9 1.562 0.1405
## 30 min NE - Control 14.79 4.36 14 5.45 24.1 3.395 0.0044
##
## Confidence level used: 0.95
```

Notes

- Note that printing the
`contrast`

object displays useful information, including the confidence level used and the method of adjustment for multiple tests. If the object is printed using`kable() %>% kable_styling()`

(as in the “Inference” and “Understanding” sections above), only the table is printed and the additional information is lost. - The
`method`

argument is used to control the set of contrasts that are computed. See below. - The
`adjust`

argument controls if and how to adjust for multiple tests. Each`method`

has a default adjustment method. See below. - The
`level`

argument controls the percentile boundaries of the confidence interval. The default is 0.95. Including this argument with this value makes this level transparent.

#### 10.2.4.1 method argument

The `method =`

argument is used to control the set of contrasts that are computed. Type `help("contrast-methods")`

into the console to see the list of available methods. Also, read the comparisons and contrasts vignette for more on `emmeans::contrast()`

`method = "trt.vs.ctrl"`

computes all non-reference minus reference contrasts. This method was used in the “Inference” section because it gives the two contrasts of the planned comparisons identified in the “understand the experimental design” step.- The default adjustment for multiple tests is “dunnettx”, which is Dunnett’s test. “none” was specified in the inference section because both comparisons were planned.

```
<- contrast(m1_emm,
m1_tukey method = "revpairwise",
adjust = "tukey",
level = 0.95) %>%
summary(infer = TRUE)
m1_tukey
```

```
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## 1 hour cold - Control 7.14 4.57 14 -4.82 19.1 1.562 0.2936
## 30 min NE - Control 14.79 4.36 14 3.39 26.2 3.395 0.0114
## 30 min NE - 1 hour cold 7.65 4.57 14 -4.31 19.6 1.675 0.2490
##
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 3 estimates
## P value adjustment: tukey method for comparing a family of 3 estimates
```

- Both
`method = "pairwise"`

and`method = "revpairwise"`

compute all pairwise comparisons. I prefer “revpairwise” because the contrasts that include the reference are in the direction non-reference minus reference. - The default adjustment for multiple tests is “tukey”, which is the Tukey’s HSD method.

#### 10.2.4.2 Adjustment for multiple tests

See the Section 15.8.3 in the chapter “Linear models with two categorical \(X\) – Factorial linear models” for a more thorough explanation of the *p*-value adjustment for multiple tests that arise if an experiment has more than two treatment groups. Here, I just outline choices available in the `emmeans::contrast`

function.

- “none” – no adjustment
- “dunnettx” – Dunnett’s test is a method used when comparing all treatments to a single control.
- “tukey” – Tukey’s HSD method is a method used to compare all pairwise comparisons.
- “bonferroni” – Bonferroni is a general purpose method to compare any set of multiple tests. The test is conservative. A better method is “holm”
- “holm” – Holm-Bonferroni is a general purpose method like the Bonferroni but is more powerful.
- “fdr” – controls the false discovery rate not the Type I error rate for a family of tests. One might use this in an exploratory experiment.
- “mvt” – based on the multivariate
*t*distribution and using covariance structure of the variables.

#### 10.2.4.3 Planned comparisons using custom contrasts

With any but purely exploratory experiment, we have certain contrasts that test the focal hypotheses. We care little to none about other contrasts. We can limit the contrasts computed using `emmeans::contrast()`

by passing custom contrasts to `method`

.

- create a set of vectors with as many elements as there are rows in the emmeans table. Name each vector using the (or a) group name for the row in the emmeans table. For the vector created for row
*j*, set the value of the*j*th element to 1 and set all other values to zero. - create a list of contrasts, with each set in the list the difference of two of the named vectors from step 1.

```
# m1_emm # print in console to get row numbers
# set the mean as the row number from the emmeans table
<- c(1,0,0) # control is in row 1
cn <- c(0,1,0) # cold in in row 2
cold <- c(0,0,1) # ne is in row 3
ne
# contrasts are the difference in the vectors created above
# the focal contrasts are in the understand the experimental
# design section
# 1. (cold - cn)
# 2. (ne - cn)
<- list(
planned_contrasts "Cold - Cn" = c(cold - cn),
"NE - Cn" = c(ne - cn)
)<- contrast(m1_emm,
m1_planned method = planned_contrasts,
adjust = "none"
%>%
) summary(infer = TRUE)
%>%
m1_planned kable(digits = c(1,1,2,0,1,1,1,4)) %>%
kable_styling()
```

contrast | estimate | SE | df | lower.CL | upper.CL | t.ratio | p.value |
---|---|---|---|---|---|---|---|

Cold - Cn | 7.1 | 4.57 | 14 | -2.7 | 16.9 | 1.6 | 0.1405 |

NE - Cn | 14.8 | 4.36 | 14 | 5.4 | 24.1 | 3.4 | 0.0044 |

### 10.2.5 How to generate ANOVA tables

ANOVA is Analysis of Variance. Researchers frequently use the term “ANOVA” as the name for an analysis of an experiment with single-factor with more than two groups. However, ANOVA is a general method of inference for complex experimental designs. ANOVA models and regression models are different ways of expressing the same underlying linear model.

Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
---|---|---|---|---|---|

treatment | 2 | 656.9 | 328.4 | 5.765 | 0.015 |

Residuals | 14 | 797.5 | 57.0 |

Notes

- Do not confuse the statistics in an ANOVA table with those in a coefficient table.
- The ANOVA table 10.5 has two rows. The first row contains the statistics for the factor \(\texttt{treatment}\). The statistics address the null hypothesis “there is no effect of treatment – all groups have the same population mean”. This is a
*p*value about the factor as a whole and not the individual comparisons of different groups within the factor. - The second row contains the statistics for the error – the residuals of the model.
- \(\texttt{F value}\) is the test statistic. It is a ratio of variances, which is why this analysis is called Analysis of Variance. The numerator variance is the \(\texttt{Mean Sq}\) of \(\texttt{treatment}\). The numerator variance is the value in $ of \(\texttt{treatment}\). The denominator variance is the value in $ of \(\texttt{Residuals}\).
- \(\texttt{Mean Sq}\) contains the
**mean square**for each term (row) in the table. The mean square is a variance. Remember that the numerator of a variance is a sum of squared differences between observed and mean values. The numerator of \(\texttt{Mean Sq}\) is the value in \(\texttt{Sum Sq}\) from the same row. And, remember that the denominator of a variance is a degree of freedom. The denominator of \(\texttt{Mean Sq}\) is the value in \(\texttt{DF}\) from the same row.

An ANOVA table for a single factor with more than two groups has a single *p*-value for the treatment term. The single *p*-value is the probability of sampling a value of *F* as large or larger than the observed *F* under the null (no true effects of either treatment and all specifications of the generating model are true). There is not much we can do with this number - we want to estimate the effect sizes and their uncertainty and we don’t get this from an ANOVA table. Many textbooks, websites, and colleagues suggest to 1) fit the ANOVA, 2) check the *F*, and, if \(F < 0.05\), 3) do “tests after an ANOVA”. These tests after an ANOVA are the planned comparisons and post-hoc tests described above using the linear model. In classical ANOVA, the initial computation of the cell means (means of treatment combinations) and sums of squares was a logical first step to the decomposition of these sums of squares to compute the contrasts. With modern linear models using regression, the ANOVA first step is unnecessary and not recommended.

If your PI, manager, thesis committee, or journal editor insists that you do ANOVA, and you cannot convince them otherwise, here is how to generate an ANOVA table in R. Note that even though we are generating that table, the computation of the contrast table and all inference from that is not part of the ANOVA.

#### 10.2.5.1 The afex aov_4 function

The package afex was developed to make it much easier for researchers to generate ANOVA tables that look like those from other statistics packages including SAS, SPSS, JMP, and Graphpad Prism.

```
# .I is a data.table function that returns the row number
:= paste("mouse", .I)]
fig2a[, fake_id
<- aov_4(diHOME ~ treatment + (1|fake_id),
m1_aov4 data = fig2a)
anova(m1_aov4)
```

```
## Anova Table (Type 3 tests)
##
## Response: diHOME
## num Df den Df MSE F ges Pr(>F)
## treatment 2 14 56.968 5.7651 0.45163 0.01491 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

Notes

- The afex package has three function names for generating the same ANOVA table and statistics – here I’m using
`aov_4`

because this functions uses a linear model formula argument (specifically, that used in the lme4 package), which is consistent with the rest of this text. - The formula includes the addition of a random factor (
`(1|id)`

) even though there really is no random factor in this model. See Section 12.1 below for a brief explanation of a random factor. The random factor (the factor variable “id” created in the line before the fit model line) identifies the individual mouse from which the response variable was measured. Because the response was only measured once on each individual mouse, “id” is not really a random factor but the addition of this in the model formula is necessary for the`aov_4`

function to work. - It is easy to get an ANOVA table that you don’t want in R. If you want an ANOVA table that matches one from Graphpad Prism or JMP or similar software, the best practice is using the ANOVA functions from the afex package.
- What do I mean by “an ANOVA table that you don’t want”? With factorial ANOVA with unbalanced data, there are three ways to compute the sums of squares for the different terms of the ANOVA table. SAS termed these Type I, II, and III sums of squares and these names have stuck. Following SAS, almost all statistics packages use Type III as the default (or only) method for computing ANOVA tables. R uses Type I as the default. There are very good arguments for using Type II. This distinction is moot for single factor ANOVA or multi-factor ANOVA for balanced designs but is not moot for unbalanced multi-factor ANOVA or any ANOVA with covariates. If you want an ANOVA table from R to match what would be generated by Graphpad Prism or JMP (Type III), then the afex package is the best practice.

#### 10.2.5.2 The car Anova function

The car package has the extremely useful `Anova`

function although using it is a bit like doing brain surgery having only watched a youtube video.

```
<- list(treatment = contr.sum)
type3 <- lm(diHOME ~ treatment,
m1_type3 data = fig2a,
contrasts = type3)
Anova(m1_type3, type="3")
```

```
## Anova Table (Type III tests)
##
## Response: diHOME
## Sum Sq Df F value Pr(>F)
## (Intercept) 6308.4 1 110.7355 4.95e-08 ***
## treatment 656.9 2 5.7651 0.01491 *
## Residuals 797.5 14
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

Notes

`car::Anova`

has arguments for reporting the Type III sum of squares. Again, this is not relevant to a single factor ANOVA with no covariates but to avoid catastrophic code in the future, its good to know about best practices now, so I’m pre-peating what is written in Section 15.11.4.2.- Background: The default model matrix in the
`lm`

function uses dummy (or treatment) coding. For a Type 3 SS ANOVA (the kind that matches that in Graphpad Prism or JMP), we need to tell`lm`

to use sum (or deviation) coding. - The best practice method for changing the contrasts in the model matrix is using the
`contrasts`

argument within the`lm`

function, as in the code above to fit`m1_type3`

. This is the safest practice because this sets the contrasts only for this specific fit. - The coefficients of
`m1_type3`

will be different from`m1`

. The intercept will be the grand mean and the coefficients of the non-reference levels (the effects) will be their deviations from the grand mean. I don’t find this definition of “effects” very useful for most experiments in biology. - The contrasts (differences in the means among pairs of groups) in the contrast table will be the same, regardless of the contrast coding.

Danger!. Many online sites suggest this bit of code before a Type III ANOVA using car::Anova:

`options(contrasts = c("contr.sum", "contr.poly")`

If you’re reading this book, you almost certainly don’t want to do this because this code resets how R computes coefficients of linear models and SS of ANOVA tables. This will effect all future analysis until the contrasts are set to something else or a new R session is started.

**base R aov and anova**

```
<- aov(diHOME ~ treatment, data = fig2a)
m1_aov summary(m1_aov)
```

```
## Df Sum Sq Mean Sq F value Pr(>F)
## treatment 2 656.9 328.4 5.765 0.0149 *
## Residuals 14 797.5 57.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 1 observation deleted due to missingness
```

```
# same as m1 in the Example 2 section
<- lm(diHOME ~ treatment, data = fig2a)
m1 anova(m1)
```

```
## Analysis of Variance Table
##
## Response: diHOME
## Df Sum Sq Mean Sq F value Pr(>F)
## treatment 2 656.85 328.43 5.7651 0.01491 *
## Residuals 14 797.55 56.97
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

Notes

- Many introduction to statistics textbooks and websites use the base R
`aov`

function. I don’t find this function useful given the afex package functions. - The base R
`anova`

is useful if you know what you are doing with it.