16 Ковариационный анализ
Вообще-то мы это уже сделали, просто не заметили.
Начнём в этот раз с данных1.
## # A tibble: 654 × 6
## ID Age FEV Height Sex Smoker
## <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 301 9 1.71 57 Female Non
## 2 451 8 1.72 67.5 Female Non
## 3 501 7 1.72 54.5 Female Non
## 4 642 9 1.56 53 Male Non
## 5 901 9 1.90 57 Male Non
## 6 1701 8 2.34 61 Female Non
## 7 1752 6 1.92 58 Female Non
## 8 1753 6 1.42 56 Female Non
## 9 1901 8 1.99 58.5 Female Non
## 10 1951 9 1.94 60 Female Non
## # … with 644 more rowsУ нас есть данные об объёме форсированного выдоха (ОФВ) (forced expiratory volume, FEV)2 у курящих и некурящих детей (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).
16.1 Модель ковариационного анализа
Модель ковариационного анализа — это модель множественной линейной регрессии, в которой отсутствуют взаимодействия дискретных и непрерывных предикторов. Мы можем пойти с другой стороны и сказать, что модель ковариационного анализа — это модель дисперсионного анализа, в которую включён (включены) один или несколько непрерывных предикторов. Непрерывный предиктор называется ковариатой.
Ковариата — это переменная, которая [потенциально] связана с нашей целевой переменной, но её влияние не является целью нашего анализа. Например, в случае, рассмотренном выше, рост — это переменная, которая, несомненно, связана с ОФВ (наша целевая переменная), но эта связь — не наша основная цель, так как мы исследуем влияние статуса курения.
Математическая модель ровным счётом та же самая, что мы видели в случае множественной линейной регрессии:
\[ \hat y_i = b_0 + b_1 I_1 + b_2 x_2 \]
Это самый простой вариант. Здесь \(I_1\) — переменная-индикатор, кодирующая категориальный предиктор (принимает значения \(0\) и \(1\)), а \(x_2\) — ковариата, или непрерывный предиктор. Как можно заметить, взаимодействие в модель не включено.
16.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Видим, что предиктор Smoker (статус курения) оказывается статистически значимым. Если мы на этом завершим наш анализ, то мы сделаем некорректный вывод: статус курения детей значимо влияет на ОФВ, при этом у группы некурящих ОФВ на 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). Что мы можем сказать, увидев такие результаты? Мы можем заключить, что разница в ОФВ между группами курящих и некурящих, которую мы видели в предыдущей модели, объясняется не статусом курения, а возрастом респондентов, вошедших в эту группу. Теперь и вывод о направлении связи не кажется таким безумным — в группу курящих, вероятно, просто в среднем старше, а соответственно, и объем легких — и ОФВ — у них больше.
Таким образом, введя ковариату в анализ, мы обнаружили, что связи между статусом курения и одним из объективных показателей здоровья дыхательной системы нет. Как это интерпретировать? Ну, это зависит. Одним из возможных вариантов может быть как раз таки возраст респондентов — «запас прочности» организма ещё не исчерпан.
Самое главное про ковариационный анализ:
Ковариационный анализ можно проводить только в том случае, если отсутствует значимое взаиможействие между факторами и ковариатами!
Если взаимодействие значимо, мы не можем исключить его из модели — придется интерпретировать результаты со взаимодействием.
16.3 Тестирование значимости предикторов и диагностика модели
Так как ANCOVA — это всё ещё линейная регрессия, то и всё, что можно делать с регрессией, можно делать и с ней. Тестирование значимости предикторов происходит точно так же, как и в предыдущих главах. Диагностика модели также аналогична.
Как видите, несмотря на то, что мы уже освоили достаточно много методов анализа данных, все они базируются на одной и той же модели — модели линейной регрессии — которая время от времени несколько модифицируется под те или иные типы данных. Спойлер, далее тоже ничего не измениться — мы будем в рамках всё той же линейной регрессии, только ход адаптации модели будет несколько более сложным. Но если сейчас вы чувствуете себя достаточно уверенно с регрессией, то дальше всё пойдёт ровно.
Описание данных из источника: «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.»↩︎
Объём воздуха, выдыхаемого за первую секунду манёвра форсированного выдоха. Используется для расчёта индекса наличия/отсутствия ухудшения проходимости дыхательных путей.↩︎