# Analysis of Covariance

This section of notes will explore an analysis procedure classically known as Analysis of Covariance (ANCOVA). This methodology is aimed to answer the following question, Is there a difference in groups, after controlling for other covariates. Here, covariates refer to another attribute in the data that is thought to influence or explain variation in the outcome, therefore, this method is meant to remove, control, or partition the variance associated with the covariate (or multiple covariates) to evaluate the conditional or adjusted means.

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.0      ✔ purrr   0.3.5
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ tidyr   1.2.1      ✔ stringr 1.4.1
## ✔ readr   2.1.3      ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::lag()    masks stats::lag()
library(ggformula)
## Loading required package: ggstance
##
## Attaching package: 'ggstance'
##
## The following objects are masked from 'package:ggplot2':
##
##     geom_errorbarh, GeomErrorbarh
##
##
## Attaching package: 'scales'
##
## The following object is masked from 'package:purrr':
##
##
##
##     col_factor
##
##
## New to ggformula?  Try the tutorials:
##  learnr::run_tutorial("introduction", package = "ggformula")
##  learnr::run_tutorial("refining", package = "ggformula")
library(broom)

theme_set(theme_bw(base_size = 18))

mutate(gestational_mean = gestational_days - mean(gestational_days),
maternal_weight_mean = maternal_pregnancy_weight - mean(maternal_pregnancy_weight),
maternal_age_mean = maternal_age - mean(maternal_age),
maternal_height_mean = maternal_height - mean(maternal_height))
## Rows: 1174 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (5): birth_weight, gestational_days, maternal_age, maternal_height, mate...
## lgl (1): maternal_smoker
##
## ℹ Use spec() to retrieve the full column specification for this data.
## ℹ Specify the column types or set show_col_types = FALSE to quiet this message.
head(baby)
## # A tibble: 6 × 10
##   birth_weight gestati…¹ mater…² mater…³ mater…⁴ mater…⁵ gesta…⁶ mater…⁷ mater…⁸
##          <dbl>     <dbl>   <dbl>   <dbl>   <dbl> <lgl>     <dbl>   <dbl>   <dbl>
## 1          120       284      27      62     100 FALSE     4.90   -28.5   -0.228
## 2          113       282      33      64     135 FALSE     2.90     6.52   5.77
## 3          128       279      28      64     115 TRUE     -0.101  -13.5    0.772
## 4          108       282      23      67     125 TRUE      2.90    -3.48  -4.23
## 5          136       286      25      62      93 FALSE     6.90   -35.5   -2.23
## 6          138       244      33      62     178 FALSE   -35.1     49.5    5.77
## # … with 1 more variable: maternal_height_mean <dbl>, and abbreviated variable
## #   names ¹​gestational_days, ²​maternal_age, ³​maternal_height,
## #   ⁴​maternal_pregnancy_weight, ⁵​maternal_smoker, ⁶​gestational_mean,
## #   ⁷​maternal_weight_mean, ⁸​maternal_age_mean

Let’s first explore the model that only includes gestational days in the model, that is, this is a simple linear regression model. Note: I mean centered gestational days for this model.

$weight = \beta_{0} + \beta_{1} days + \epsilon$

baby_lm <- lm(birth_weight ~ gestational_mean, data = baby)

glance(baby_lm)
## # A tibble: 1 × 12
##   r.squared adj.r.squa…¹ sigma stati…²  p.value    df logLik   AIC   BIC devia…³
##       <dbl>        <dbl> <dbl>   <dbl>    <dbl> <dbl>  <dbl> <dbl> <dbl>   <dbl>
## 1     0.166        0.165  16.7    233. 3.40e-48     1 -4973. 9953. 9968. 328608.
## # … with 2 more variables: df.residual <int>, nobs <int>, and abbreviated
## #   variable names ¹​adj.r.squared, ²​statistic, ³​deviance
tidy(baby_lm) |>
mutate_if(is.double, round, 3)
## # A tibble: 2 × 5
##   term             estimate std.error statistic p.value
##   <chr>               <dbl>     <dbl>     <dbl>   <dbl>
## 1 (Intercept)       119.        0.489     244.        0
## 2 gestational_mean    0.467     0.031      15.3       0

Let’s next explore another simple linear regression model that only includes whether a mother smoked or not. This model is:

$weight = \beta_{0} + \beta_{2} smoker + \epsilon$

baby_smoker <- lm(birth_weight ~ maternal_smoker, data = baby)

glance(baby_smoker)
## # A tibble: 1 × 12
##   r.squared adj.r.sq…¹ sigma stati…²  p.value    df logLik    AIC    BIC devia…³
##       <dbl>      <dbl> <dbl>   <dbl>    <dbl> <dbl>  <dbl>  <dbl>  <dbl>   <dbl>
## 1    0.0609     0.0601  17.8    76.0 9.46e-18     1 -5043. 10092. 10107. 370056.
## # … with 2 more variables: df.residual <int>, nobs <int>, and abbreviated
## #   variable names ¹​adj.r.squared, ²​statistic, ³​deviance
tidy(baby_smoker) |>
mutate_if(is.double, round, 3)
## # A tibble: 2 × 5
##   term                estimate std.error statistic p.value
##   <chr>                  <dbl>     <dbl>     <dbl>   <dbl>
## 1 (Intercept)           123.       0.665    185.         0
## 2 maternal_smokerTRUE    -9.27     1.06      -8.72       0

The idea behind ANCOVA is that we know that gestational days is an important predictor in understanding variation in the birth weight of the baby. This may not be a predictor that is of most interest to our research question of knowing whether or not there are differences in birth weight between those that smoked vs did not smoke when the baby was in the womb. There may be some sampling bias occurring, that is, those that smoked may have fewer gestational days compared to those that did not smoke. This fact may drive the differences in the mean difference of birth weight. As a results, we would want to adjust for the effect of gestational days prior to evaluating the mean difference of smoker status.

Below is a descriptive figure trying to show the relationship between birth weight and gestational days by smoker status.

gf_point(birth_weight ~ gestational_days, data = baby, color = ~maternal_smoker, size = 4) %>%
gf_smooth(method = 'lm', size = 1)
## Warning: Using size aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use linewidth instead. Below, let’s fit the traditional ANCOVA model, that contains a single continuous covariate (gestational days) and the categorical predictor. The model would be:

$weight = \beta_{0} + \beta_{1} days + \beta_{2} smoker + \epsilon$

baby_ancova <- lm(birth_weight ~ gestational_mean + maternal_smoker, data = baby)

glance(baby_ancova)
## # A tibble: 1 × 12
##   r.squared adj.r.squa…¹ sigma stati…²  p.value    df logLik   AIC   BIC devia…³
##       <dbl>        <dbl> <dbl>   <dbl>    <dbl> <dbl>  <dbl> <dbl> <dbl>   <dbl>
## 1     0.216        0.214  16.2    161. 1.71e-62     2 -4937. 9883. 9903. 309075.
## # … with 2 more variables: df.residual <int>, nobs <int>, and abbreviated
## #   variable names ¹​adj.r.squared, ²​statistic, ³​deviance
tidy(baby_ancova) |>
mutate_if(is.double, round, 3)
## # A tibble: 3 × 5
##   term                estimate std.error statistic p.value
##   <chr>                  <dbl>     <dbl>     <dbl>   <dbl>
## 1 (Intercept)          123.        0.608    202.         0
## 2 gestational_mean       0.451     0.03      15.2        0
## 3 maternal_smokerTRUE   -8.37      0.973     -8.60       0

The model above can be partioned to better show the specific terms in the model.

$weight_{smoker} = (\beta_{0} + \beta_{2}) + \beta_{1} days + \epsilon$

$weight_{non-smoker} = \beta_{0} + \beta_{1} days + \epsilon$

More explicitly, the term, $$\beta_{2}$$ above reflects a mean level change on the intercept after controlling for the effects of gestational days.

Let’s compare the means from this ANCOVA model compared to the unconditional means.

uncond_mean <- df_stats(birth_weight ~ maternal_smoker, data = baby, mean) %>%
mutate(maternal_smoker = as.character(maternal_smoker))

cond_mean <- data.frame(maternal_smoker = c('FALSE', 'TRUE'), cond_mean = c(122.74, 122.74 - 8.37))

combine_means <- full_join(uncond_mean, cond_mean)
## Joining, by = "maternal_smoker"
combine_means
##       response maternal_smoker     mean cond_mean
## 1 birth_weight           FALSE 123.0853    122.74
## 2 birth_weight            TRUE 113.8192    114.37

Note, the conditional mean used the fact that the gestational days were mean centered, therefore, a value of 0 was included in the above equations.

## Interaction between continuous and categorical predictor

One assumption in classical ANCOVA is to assume that the effect of the adjusting covariate is the same for each group. Although, this is a classical assumption, this is an empirical question that we can explore given the data. This can be explored by the introduction of an interaction.

This new model would be:

$weight = \beta_{0} + \beta_{1} days + \beta_{2} smoker + \beta_{3} days:smoker + \epsilon$

where the $$:$$ in the equation above represents the interaction effect between gestational days and smoker status. An interaction is testing if the slope for gestational days is the same for those that smoke vs those that did not smoke.

More explicitly, we can partition the model above into the following two regression equations.

$weight_{smoker} = (\beta_{0} + \beta_{2}) + (\beta_{1} + \beta_{3}) days + \epsilon$

$weight_{non-smoker} = \beta_{0} + \beta_{1} days + \epsilon$

baby_lm_int <- lm(birth_weight ~ gestational_mean * maternal_smoker, data = baby)
# baby_lm_int <- lm(birth_weight ~ gestational_mean + maternal_smoker + gestational_mean:maternal_smoker,
#    data = baby)

glance(baby_lm_int)
## # A tibble: 1 × 12
##   r.squared adj.r.squa…¹ sigma stati…²  p.value    df logLik   AIC   BIC devia…³
##       <dbl>        <dbl> <dbl>   <dbl>    <dbl> <dbl>  <dbl> <dbl> <dbl>   <dbl>
## 1     0.225        0.223  16.2    113. 2.41e-64     3 -4930. 9871. 9896. 305427.
## # … with 2 more variables: df.residual <int>, nobs <int>, and abbreviated
## #   variable names ¹​adj.r.squared, ²​statistic, ³​deviance
tidy(baby_lm_int) |>
mutate_if(is.double, round, 3)
## # A tibble: 4 × 5
##   term                                 estimate std.error statistic p.value
##   <chr>                                   <dbl>     <dbl>     <dbl>   <dbl>
## 1 (Intercept)                           123.        0.605    203.         0
## 2 gestational_mean                        0.37      0.037     10.1        0
## 3 maternal_smokerTRUE                    -8.26      0.969     -8.52       0
## 4 gestational_mean:maternal_smokerTRUE    0.231     0.062      3.74       0

What do the 4 terms represent now?

1. $$\beta_{0}$$ is the intercept for the reference group (non-smokers) when all other attributes in the model are equal to 0.
2. $$\beta_{1}$$ is the slope for the reference group (non-smokers), controlling for the smoker status.
3. $$\beta_{2}$$ is the average change for those that smoke, controlling for differences in gestational days.
4. $$\beta_{3}$$ is the average slope difference for those that smoke compared to those that did not smoke.

Often, it is easier to combine the terms to be compared directly. For example, given the example above, we could combine $$\beta_{0}$$ and $$\beta_{2}$$ to get the estimated intercept for those that smoke. Then we could do something similar for the slope terms for those that smoke.

### Hypotheses being tested

The hypotheses being tested here for the $$\beta_{2}$$ and $$\beta_{3}$$ are as follows:

$H_{0}: \beta_{2} = 0$

and

$H_{0}: \beta_{3} = 0$

These can be a little bit untenable, therefore stating them in words can help.

• $$H_{0}: \beta_{2} = 0$$, is really testing, is there a mean difference between the reference and focal groups, controlling for the other terms in the model. This is inherently a conditional or adjusted means question.
• $$H_{0}: \beta_{3} = 0$$, is really testing, is there a slope difference comparing the focal and reference groups, controlling for other terms in the model.

## Including more than 1 covariate to adjust/condition

It is also possible to adjust for more than one covariate. For example, maybe we thought adjusting for the weight, height, and age of the mother would also be a good idea. This can be done to remove any variation associated with that attribute prior to making exploring the mean difference in the attributes of interest.

This model would now be:

$weight = \beta_{0} + \beta_{1} days + \beta_{4} mother\_weight + \beta_{5} mother\_height + \beta_{6} mother\_age + \beta_{2} smoker + \beta_{3} days:smoker + \epsilon$

These can also be partitioned into two separate equations.

$weight_{smoker} = (\beta_{0} + \beta_{2}) + (\beta_{1} + \beta_{3}) days + \beta_{4} mother\_weight + \beta_{5} mother\_height + \beta_{6} mother\_age + \epsilon$

$weight_{non-smoker} = \beta_{0} + \beta_{1} days + \beta_{4} mother\_weight + \beta_{5} mother\_height + \beta_{6} mother\_age + \epsilon$

baby_lm_int_w <- lm(birth_weight ~ gestational_mean + maternal_pregnancy_weight + maternal_height + maternal_age + maternal_smoker + gestational_mean:maternal_smoker,
data = baby)

glance(baby_lm_int_w)
## # A tibble: 1 × 12
##   r.squared adj.r.squa…¹ sigma stati…²  p.value    df logLik   AIC   BIC devia…³
##       <dbl>        <dbl> <dbl>   <dbl>    <dbl> <dbl>  <dbl> <dbl> <dbl>   <dbl>
## 1     0.260        0.256  15.8    68.3 6.14e-73     6 -4903. 9822. 9863. 291631.
## # … with 2 more variables: df.residual <int>, nobs <int>, and abbreviated
## #   variable names ¹​adj.r.squared, ²​statistic, ³​deviance
tidy(baby_lm_int_w) |>
mutate_if(is.double, round, 3)
## # A tibble: 7 × 5
##   term                                 estimate std.error statistic p.value
##   <chr>                                   <dbl>     <dbl>     <dbl>   <dbl>
## 1 (Intercept)                            43.5      12.3       3.54    0
## 2 gestational_mean                        0.365     0.036    10.1     0
## 3 maternal_pregnancy_weight               0.052     0.025     2.06    0.039
## 4 maternal_height                         1.10      0.204     5.41    0
## 5 maternal_age                            0.07      0.081     0.873   0.383
## 6 maternal_smokerTRUE                    -8.20      0.952    -8.62    0
## 7 gestational_mean:maternal_smokerTRUE    0.208     0.061     3.44    0.001

How could the appropriate conditional means be computed here? Also, what happened to the intercept in this model, why is it so small?

Previous