20 Ковариационный анализ
Это будет очень короткая глава, так как всё, что нужно знать о ковариационном анализе, мы уже знаем. Просто пока не знаем, что то, что мы знаем, это ковариационный анализ. 🙃
Вспомним данные, с которым работали в главах о простой и множественной линейной регрессии.
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──## ✓ ggplot2 3.3.2 ✓ purrr 0.3.4
## ✓ tibble 3.0.4 ✓ dplyr 1.0.2
## ✓ tidyr 1.1.2 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.0## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## ID = col_double(),
## Age = col_double(),
## FEV = col_double(),
## Height = col_double(),
## Sex = col_character(),
## Smoker = col_character()
## )## tibble [654 × 6] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ ID : num [1:654] 301 451 501 642 901 ...
## $ Age : num [1:654] 9 8 7 9 9 8 6 6 8 9 ...
## $ FEV : num [1:654] 1.71 1.72 1.72 1.56 1.9 ...
## $ Height: num [1:654] 57 67.5 54.5 53 57 61 58 56 58.5 60 ...
## $ Sex : chr [1:654] "Female" "Female" "Female" "Male" ...
## $ Smoker: chr [1:654] "Non" "Non" "Non" "Non" ...
## - attr(*, "spec")=
## .. cols(
## .. ID = col_double(),
## .. Age = col_double(),
## .. FEV = col_double(),
## .. Height = col_double(),
## .. Sex = col_character(),
## .. Smoker = col_character()
## .. )Описание данных из источника: «FEV (forced expiratory volume) is an index of pulmonary function that measures the volume of air expelled after one second of constant effort. The data contains determinations of FEV on 654 children ages 6-22 who were seen in the Childhood Respiratory Desease Study in 1980 in East Boston, Massachusetts. The data are part of a larger study to follow the change in pulmonary function over time in children.»
Файл содержит данные об объёме форсированного выдоха (ОФВ) (forced expiratory volume, FEV)1 у курящих и некурящих детей (6–22 лет) и информацию о поле, возрасте, росте и ID пациента.
Поставим задачу: мы хотим исследовать, влияет ли статус курения на объем формированного выдоха у детей. Ожидаем, что влияет. Однако мы также понимаем, что рост ребенка будет связан с объёмом его лёгких, а значит, и с нашей целевой переменной.
Построим соответствующую модель. В качестве целевой переменной возьмем ОФВ (FEV), а в качестве предикторов статус курения (Smoker) и рост (Height), а также их взаимодействие, так как оно может быть потенциально значимо:
##
## Call:
## lm(formula = FEV ~ Smoker * Height, data = ofv)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.74332 -0.26960 -0.00462 0.23909 2.12941
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.40165 1.11228 -6.654 6.04e-11 ***
## SmokerNon 2.03678 1.12846 1.805 0.0716 .
## Height 0.16191 0.01685 9.612 < 2e-16 ***
## SmokerNon:Height -0.03106 0.01713 -1.813 0.0703 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4302 on 650 degrees of freedom
## Multiple R-squared: 0.7549, Adjusted R-squared: 0.7538
## F-statistic: 667.3 on 3 and 650 DF, p-value: < 2.2e-16Модель сошлась, статистически значима, объясняет 75% дисперсии данных. Это очень хорошо.
Мы также говорили, что чем проще модель, тем лучше, и учились упрощать модели. Мы видим, что взаимодействие предикторов оказалось статистически незначимым. Посмотрим, влияет ли его удаление на информативность модели:
## Single term deletions
##
## Model:
## FEV ~ Smoker * Height
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 120.32 -1099.2
## Smoker:Height 1 0.60843 120.93 -1097.9 3.2868 0.0703 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Частный F-критерий показывает, что удаление взаимодействия из модели не повлияет на её информативность, следовательно, взаиможействие можно удалить — оно не значимо.
##
## Call:
## lm(formula = FEV ~ Smoker + Height, data = ofv)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.7505 -0.2660 -0.0041 0.2447 2.1207
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.421300 0.210104 -25.803 <2e-16 ***
## SmokerNon -0.006319 0.058686 -0.108 0.914
## Height 0.131883 0.003081 42.808 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.431 on 651 degrees of freedom
## Multiple R-squared: 0.7537, Adjusted R-squared: 0.7529
## F-statistic: 995.9 on 2 and 651 DF, p-value: < 2.2e-16Внимательно на неё посмотрим. И познакомимся. Ведь мы только что построили модель ковариационного анализа (analysis of covariance, ANCOVA).
20.1 Модель ковариационного анализа
Модель ковариационного анализа — это модель множественной линейной регрессии, в которой отсутствуют взаимодействия дискретных и непрерывных предикторов. Мы можем пойти с другой стороны и сказать, что модель ковариационного анализа — это модель дисперсионного анализа, в которую включён (включены) один или несколько непрерывных предикторов. Непрерывный предиктор называется ковариатой.
Ковариата — это переменная, которая [потенциально] связана с нашей целевой переменной, но её влияние не является целью нашего анализа. Например, в случае, рассмотренном выше, рост — это переменная, которая, несомненно, связана с ОФВ, но эта связь — не наша основная цель, так как мы исследуем влияние статуса курения.
Математическая модель ровным счётом та же самая, что мы видели в случае множественной линейной регрессии:
\[ y_i = b_0 + b_1 I_{\mathrm{Group2}} + b_2 x_2 + e = b_0 + b_1 x_1 + b_2 x_2 \]
Это самый простой вариант. Здесь \(x_1\) — категориальный предиктор (принимает значения \(0\) и \(1\)), а \(x_2\) — ковариата, или непрерывный предиктор. Как можно заметить, взаимодействие в модель не включено.
20.2 Влияние ковариаты
Зачем вообще учитывать ковариаты? Рассмотрим такую модель:
##
## Call:
## lm(formula = FEV ~ Smoker, data = ofv)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.7751 -0.6339 -0.1021 0.4804 3.2269
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.2769 0.1043 31.407 < 2e-16 ***
## SmokerNon -0.7107 0.1099 -6.464 1.99e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8412 on 652 degrees of freedom
## Multiple R-squared: 0.06023, Adjusted R-squared: 0.05879
## F-statistic: 41.79 on 1 and 652 DF, p-value: 1.993e-10## Anova Table (Type II tests)
##
## Response: FEV
## Sum Sq Df F value Pr(>F)
## Smoker 29.57 1 41.789 1.993e-10 ***
## Residuals 461.35 652
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Видим, что предиктор «статус курения» оказывается статистически значимым. Если мы на этом завершим наш анализ, то мы сделаем некорректный вывод: статус курения детей значимо влияет на ОФВ, при этом у группы некурящих ОФВ на 0.7 л меньше. Здесь нас наталкивает на мысль об ошибке направление влияния нашего предиктора — согласно здравому смыслу и биологии и медицины ОФВ должен быть ниже у курящей группы. Но мы можем получить результаты, адекватные текущей реальности и ничего не заподозрить.
Проблема здесь кроется в том, что, как мы уде сказали, ОФВ связан не только со статусом курения. Ковариата, как и любой предиктор, будет объяснять часть дисперсии данных, поэтому значимость других предикторов модели при её введении может изменяться.
Взглянем ещё раз на модель ANCOVA:
##
## Call:
## lm(formula = FEV ~ Smoker + Height, data = ofv)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.7505 -0.2660 -0.0041 0.2447 2.1207
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.421300 0.210104 -25.803 <2e-16 ***
## SmokerNon -0.006319 0.058686 -0.108 0.914
## Height 0.131883 0.003081 42.808 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.431 on 651 degrees of freedom
## Multiple R-squared: 0.7537, Adjusted R-squared: 0.7529
## F-statistic: 995.9 on 2 and 651 DF, p-value: < 2.2e-16Ковариата значима, а вот наш основной предиктор свою значимость утратил — причем практически абсолютно (p > 0.9). Что мы можем сказать, увидев такие результаты? Мы можем заключить, что разница в ОФВ между группами курящих и некурящих, которую мы видели в model3 объясняется не статусом курения, а возрастом респондентов, вошедших в эту группу. Теперь и вывод о направлении связи не кажется таким безумным — в группа курящих, вероятно, просто в среднем старше, а соответственно, и объем легких, следовательно и ОФВ, у них больше.
Проверьте данное предположение, сравнив рост и возраст в двух группах.
Таким образом, введя ковариату в анализ, мы обнаружили, что связи между статусом курения и одним из объективных показателей здоровья дыхательной системы нет. Как это интерпретировать? Ну, ça dépend. Одним из возможных вариантов может быть как раз таки возраст респондентов — «запас прочности» организма ещё не исчерпан.
20.3 Более сложные модели
В модель можно ввести несколько дискретных предикторов, их взаиможействия, а также несколько ковариат. Однако надо быть внимательным.
Пусть мы хотим немного усложнить нашу модель и ввести в неё одним из предикторов пол. Итого у нас будут следующие предикторы:
- основные:
- пол (
Sex) - статус курения (
Smoker) - их взаимодействие (
Sex:Smoker)
- пол (
- ковариаты:
- рост (
Height)
- рост (
Однако нам придется начать с полной модели, чтобы понять, можно ли вообще проводить здесь ковариационный анализ.
Ковариационный анализ можно проводить только в том случае, если отсутствует значимое взаиможействие между факторами и ковариатами!
Строим полную модель:
##
## Call:
## lm(formula = FEV ~ Sex * Smoker * Height, data = ofv)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.55289 -0.25070 0.00711 0.24854 2.03200
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.02132 1.90900 -0.011 0.991093
## SexMale -10.27579 2.59616 -3.958 8.39e-05 ***
## SmokerNon -4.37701 1.93495 -2.262 0.024024 *
## Height 0.04628 0.02956 1.566 0.117886
## SexMale:SmokerNon 8.96579 2.62577 3.415 0.000679 ***
## SexMale:Height 0.16002 0.03925 4.077 5.13e-05 ***
## SmokerNon:Height 0.06743 0.03002 2.246 0.025048 *
## SexMale:SmokerNon:Height -0.13649 0.03978 -3.431 0.000640 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4173 on 646 degrees of freedom
## Multiple R-squared: 0.7708, Adjusted R-squared: 0.7683
## F-statistic: 310.4 on 7 and 646 DF, p-value: < 2.2e-16Саммари нам как бы намекает, что это самая простая модель из возможных, но попытаемся её упростить. Рассчитаем частный F-критерий:
## Single term deletions
##
## Model:
## FEV ~ Sex * Smoker * Height
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 112.52 -1135.0
## Sex:Smoker:Height 1 2.0503 114.57 -1125.2 11.771 0.0006398 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Значимо взаиможействие самого высокого уровня, а значит его нельзя исключить из модели. Поэтому ковариационный анализ в этом случае невозможен. 😢
20.4 Тестирование значимости предикторов и диагностика модели
Так как ANCOVA — это всё ещё линейная регрессия, то и всё, что можно делать с регрессией, можно делать и с ней. Тестирование значимости предикторов происходит точно так же, как и две главы назад. Диагностика модели также аналогична. Более того, если возникает потребность, то можно воспользоваться и функцией predict(), чтобы предсказать значения на новых данных. Всё абсолютно то же самое.
Как видите, несмотря на то, что мы уже освоили достаточно много методов анализа данных, все они базируются на одной и той же модели — модели линейной регрессии — которая время от времени несколько модифицируется под те или иные типы данных. Спойлер, далее тоже ничего не измениться — мы будем в рамках всё той же линейной регрессии, только ход адаптации модели будет несколько более сложным. Но если сейчас вы чувствуете себя достаточно уверенно, то и дальше всё пойдёт ровно. 😉
Объём воздуха, выдыхаемого за первую секунду манёвра форсированного выдоха. Используется для расчёта индекса наличия/отсутствия ухудшения проходимости дыхательных путей.↩︎