- Линейные модели c непрерывными и дискретными предикторами
- Взаимодействие дискретных и непрерывных предикторов
Мы рассмотрим
Геометрическая интерпретация коэффициентов регрессии
Но ведь может быть еще и так
- Как все эти изменения закодировать в уравнении линейной модели?
Похожий интерсепт, похожий угол наклона
Если зависимости для обеих групп выглядят одинаково, то мы можем пользоваться общей линейной регрессией, не учитывать группу
\(Y = b_0 + b_1X\)
Разный интерсепт, похожий угол наклона
\(Y = b_0 + b_1X + b_2 I_{Group2}\)
\(I_{Group2} = 1\) если наблюдение из группы 2, и \(I_{Group2} = 0\) в остальных случаях
Уравнение для группы 1: \(Y = b_0 + b_1X + b_2\cdot0\)
Уравнение для группы 2: \(Y = b_0 + b_1X + b_2\cdot1 =(b_0 + b_2) + b_1X\)
Похожий интерсепт, разный угол наклона
\(Y = b_0 + b_1X + b_2X\cdot I_{Group2}\)
\(I_{Group2}\) если наблюдение из группы 2, и \(I_{Group2} = 0\) в остальных случаях
Уравнение для группы 1: \(Y = b_0 + b_1X + b_2X \cdot 0 = b_0 + b_1X\)
Уравнение для группы 2: \(Y = b_0 + b_1X + b_2X \cdot 1 = b_0 + (b_1 + b_2)X\)
Разный интерсепт и разный угол наклона
\(Y = b_0 + b_1X + b_2 \cdot I_{Group2} + b_3X \cdot I_{Group2}\)
\(I_{Group2} = 1\) если наблюдение из группы 2, и \(I_{Group2} = 0\) в остальных случаях
Уравнение для группы 1: \(Y = b_0 + b_1X + b_2 \cdot 0 + b_3X \cdot 0 = b_0 + b_1X\)
Уравнение для группы 2: \(Y = b_0 + b_1X + b_2 \cdot 1 + b_3X \cdot 1 = (b_0 + b_2) + (b_1 + b_3)X\)
Физический смысл взаимодействия факторов
Не называйте взаимодействие факторов “корреляцией” или “связью”!
Взаимодействие предикторов означает, что для разных значений одного из факторов зависимость от другого предиктора выглядит по-разному.
Т.е. в данном случае можно сказать: “Для группы 1 и группы 2 зависимость Y от Х выглядит по-разному”
Если у дискретного предиктора становится больше уровней…
\[y_i = b_0 + b_1x_i + \big(b_2I_{2i} + ... + b_k I_{ki}\big) + \big(b_{k+1}I_{2i}x_i + ... + b_{2k-1} I_{ki} x_i\big) + e_i\]
В случае, если у дискретного фактора \(k\) градаций, уравнение полной модели будет включать несколько групп коэффициентов:
- \(b_0\) — значение свободного члена для базового уровня дискретного фактора
- \(b_1\) — коэффициент угла наклона для базового уровня дискретного фактора
Поправки для свободного члена:
- \(b_2, ..., b_k\) — коррекция свободного члена для других уровней (всего их \(k - 1\))
Поправки для угла наклона:
- \(b_{k+1}, ..., b_{2k - 1}\) — коррекция коэффициентов угла наклона для других уровней (всего их \(k - 1\))
Модель со взаимодействием на примере данных о весе новорожденных
Вес новорожденных и курение
Как вес новорожденных зависит от возраста матери и того, курит ли она
age— возраст материlwt— вес матери до беременностиrace— раса (1-белые, 2-черные, 3-другие)smoke— курение во время беременности (1-да,2-нет)ptl— число предыдущих преждевременных родовht— гипертензияui— гипертонус маткиftv— число визитов к врачу в последний триместрbwt— вес новорожденного, г
wt <- read.table("data/birthwt.csv", header = TRUE, sep = ";")
Преобразуем значения перменных
wt$smoke[wt$smoke == 1] <- 'Smoker' wt$smoke[wt$smoke == 0] <- 'Non smoker' wt$smoke <- factor(wt$smoke)
Задание
- Исследуйте данные о весе новорожденных
- Постройте модель зависимости веса новорожденных от возраста матери, того курит она или нет и взаимодействия предикторов
- Проверьте условия применимости этой модели
- Упростите модель, если это возможно
- Напишите общее уравнение и отдельные уравнения модели
- Постройте график предсказаний модели
Решение
str(wt)
## 'data.frame': 189 obs. of 9 variables: ## $ age : int 19 33 20 21 18 21 22 17 29 26 ... ## $ lwt : int 182 155 105 108 107 124 118 103 123 113 ... ## $ race : int 2 3 1 1 1 3 1 3 1 1 ... ## $ smoke: Factor w/ 2 levels "Non smoker","Smoker": 1 1 2 2 2 1 1 1 2 2 ... ## $ ptl : int 0 0 0 0 0 0 0 0 0 0 ... ## $ ht : int 0 0 0 0 0 0 0 0 0 0 ... ## $ ui : int 1 0 0 1 1 0 0 0 0 0 ... ## $ ftv : int 0 3 1 2 0 0 1 1 1 0 ... ## $ bwt : int 2523 2551 2557 2594 2600 2622 2637 2637 2663 2665 ...
colSums(is.na(wt))
## age lwt race smoke ptl ht ui ftv bwt ## 0 0 0 0 0 0 0 0 0
Решение
table(wt$smoke, wt$age)
## ## 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 ## Non smoker 2 3 5 7 3 8 10 7 8 9 10 12 3 2 5 4 5 2 4 ## Smoker 1 0 2 5 7 8 8 5 5 4 3 3 5 1 4 3 2 3 2 ## ## 33 34 35 36 45 ## Non smoker 2 0 1 2 1 ## Smoker 1 1 1 0 0
table(wt$smoke, wt$race)
## ## 1 2 3 ## Non smoker 44 16 55 ## Smoker 52 10 12
table(wt$smoke, wt$ui)
## ## 0 1 ## Non smoker 100 15 ## Smoker 61 13
Решение
wt$race <- factor(wt$race) wt$smoke <- factor(wt$smoke) wt$ptl <- factor(wt$ptl) wt$ht <- factor(wt$ht) wt$ui <- factor(wt$ui) wt$ftv <- factor(wt$ftv)
Решение
gg_dot <- ggplot(wt, aes(y = 1:nrow(wt))) + geom_point() gg_dot + aes(x = age)
Решение
gg_dot + aes(x = bwt)
Решение
Проверка на колинеарность
wt_mod_1 <- lm(bwt ~ age+smoke, data = wt) library(car) vif(wt_mod_1)
## age smoke ## 1 1
- Колинеарности нет
Проверка на гомогенность углов наклона
Подберем полную модель.
wt_mod_2 <- lm(bwt ~ age*smoke, data = wt) drop1(wt_mod_2, test = "F")
## Single term deletions ## ## Model: ## bwt ~ age * smoke ## Df Sum of Sq RSS AIC F value Pr(>F) ## <none> 93062605 2485 ## age:smoke 1 2609683 95672288 2488 5.19 0.024 * ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
- От исключения взаимодействия модель становится значительно хуже. Оставляем
Решение
# Данные для графиков остатков
wt_mod_2_diag <- fortify(wt_mod_2)
wt_mod_2_diag <- data.frame(
wt_mod_2_diag,
wt[, c("lwt", "race", "smoke", "ptl", "ht", "ui", "ftv")])
График расстояния Кука
ggplot(wt_mod_2_diag, aes(x = 1:nrow(wt_mod_2_diag), y = .cooksd)) + geom_bar(stat = "identity")
Решение
График остатков от предсказанных значений
gg_resid <- ggplot(data = wt_mod_2_diag, aes(x = .fitted, y = .stdresid)) + geom_point() + geom_hline(yintercept = 0) + geom_smooth() gg_resid
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Решение
Графики остатков от предикторов в модели и не в модели
library(gridExtra)
grid.arrange(gg_resid + aes(x = age),
gg_resid + aes(x = lwt),
nrow = 1)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x' ## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Решение
gg_box <- ggplot(data = wt_mod_2_diag, aes(x = smoke, y = .stdresid)) +
geom_boxplot() + geom_hline(yintercept = 0)
grid.arrange(gg_box + aes(x = smoke),
gg_box + aes(x = ftv),
gg_box + aes(x = race),
gg_box + aes(x = ht),
gg_box + aes(x = ui),
nrow = 2)
Решение
Квантильный график остатков
qqPlot(wt_mod_2)
## [1] 131 132
Проверяем значимость коэффициентов
Подумайте, что означают эти коэффициенты
summary(wt_mod_2)
## ## Call: ## lm(formula = bwt ~ age * smoke, data = wt) ## ## Residuals: ## Min 1Q Median 3Q Max ## -2189.3 -458.5 51.5 527.3 1521.4 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 2406.1 292.2 8.23 3.2e-14 *** ## age 27.7 12.1 2.28 0.024 * ## smokeSmoker 798.2 484.3 1.65 0.101 ## age:smokeSmoker -46.6 20.4 -2.28 0.024 * ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 709 on 185 degrees of freedom ## Multiple R-squared: 0.0691, Adjusted R-squared: 0.054 ## F-statistic: 4.58 on 3 and 185 DF, p-value: 0.00407
Таблица дисперсионного анализа
Anova(wt_mod_2, type = 3)
## Anova Table (Type III tests) ## ## Response: bwt ## Sum Sq Df F value Pr(>F) ## (Intercept) 34110331 1 67.81 3.2e-14 *** ## age 2620946 1 5.21 0.024 * ## smoke 1366143 1 2.72 0.101 ## age:smoke 2609683 1 5.19 0.024 * ## Residuals 93062605 185 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Анатомируем модель
Результаты подбора модели
summary(wt_mod_2)
## ## Call: ## lm(formula = bwt ~ age * smoke, data = wt) ## ## Residuals: ## Min 1Q Median 3Q Max ## -2189.3 -458.5 51.5 527.3 1521.4 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 2406.1 292.2 8.23 3.2e-14 *** ## age 27.7 12.1 2.28 0.024 * ## smokeSmoker 798.2 484.3 1.65 0.101 ## age:smokeSmoker -46.6 20.4 -2.28 0.024 * ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 709 on 185 degrees of freedom ## Multiple R-squared: 0.0691, Adjusted R-squared: 0.054 ## F-statistic: 4.58 on 3 and 185 DF, p-value: 0.00407
Можем ли мы избавиться от взаимодействия?
drop1(wt_mod_2, test = 'F')
## Single term deletions ## ## Model: ## bwt ~ age * smoke ## Df Sum of Sq RSS AIC F value Pr(>F) ## <none> 93062605 2485 ## age:smoke 1 2609683 95672288 2488 5.19 0.024 * ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Если удалить взаимодействия, то модель значимо изменится.
В этой системе присутствует взаимодействие!
Модельная матрица при наличии взаимодействий дискретного и непрерывного предиктора
X <- model.matrix(wt_mod_2) head(X, 10)
## (Intercept) age smokeSmoker age:smokeSmoker ## 1 1 19 0 0 ## 2 1 33 0 0 ## 3 1 20 1 20 ## 4 1 21 1 21 ## 5 1 18 1 18 ## 6 1 21 0 0 ## 7 1 22 0 0 ## 8 1 17 0 0 ## 9 1 29 1 29 ## 10 1 26 1 26
В колонке (Intercept) — как всегда единицы.
В колонке age — возраст матери
Базовый уровень smoke = Non smoker, в колонке smokeSmoker стоят единицы, если мать курящая, если она не курит, то стоят нули.
Колонка age:smokeSmoker описывает взаимодействие: если мать курит (Smoker), то в этой колонке стоит ее возраст, если мать не курит (Non smoker), то стоит ноль.
Вектор коэффициентов
b <- coefficients(wt_mod_2) b
## (Intercept) age smokeSmoker age:smokeSmoker ## 2406.0580 27.7314 798.1749 -46.5719
Предсказанные значения в матричном виде
\[ \widehat{\textbf{y}} = \textbf{X} \textbf{b} \]
as.numeric(fitted(wt_mod_2)[1])
## [1] 2932.95
(X %*% b) [1]
## [1] 2932.95
Предсказанное значение для объекта, относящегося к базовому уровню
X[1,] # Первая строка в модельной матрице
## (Intercept) age smokeSmoker age:smokeSmoker ## 1 19 0 0
b
## (Intercept) age smokeSmoker age:smokeSmoker ## 2406.0580 27.7314 798.1749 -46.5719
Что происходит при матричном умножении
1 * 2406.0580 + 19.0 * 27.7314 + 0 * 798.1749 + 0 * -46.5719
## [1] 2932.95
Общий вид формулы предсказанного значения для базового уровня
\[ \widehat{bwt_i} = 2406.0580 + 27.7314 age_i \]
Предсказанное значение для объекта, не относящегося к базовому уровню
as.numeric(fitted(wt_mod_2)[3])
## [1] 2827.42
X[3,] # Третья строка в модельной матрице
## (Intercept) age smokeSmoker age:smokeSmoker ## 1 20 1 20
b # Коэффциенты
## (Intercept) age smokeSmoker age:smokeSmoker ## 2406.0580 27.7314 798.1749 -46.5719
Что происходит при матричном умножении
1 * 2406.0580 + 20 * 27.7314 + 1 * 798.1749 + 20 * -46.5719
## [1] 2827.42
Общий вид формулы предсказанного значения для объектов, не относящихся к базовому уровню
\[ \widehat{bwt}_i = 2406.0580 + 27.7314 age_i + 798.1749 + (-46.5719) age_i \]
\[ \widehat{bwt}_i = (2406.0580 + 798.1749) + (27.7314 - 46.5719) age_i = 3204.23 +(-18.8405)age_i \]
Визуализация модели
Сложные модели лучше по возможности изображать в виде графика
Данные для графика
# Диапазон возрастов разный для курящих и некурящих, поэтому... library(dplyr) new_data <- wt %>% group_by(smoke)%>% do(data.frame(age = seq(min(.$age), max(.$age), length.out = 100))) # Предсказанные значения Predictions <- predict(wt_mod_2, newdata = new_data, se.fit = TRUE) new_data$fit <- Predictions$fit # Стандартные ошибки new_data$se <- Predictions$se.fit # Критические значения t t_crit <- qt(0.975, df = nrow(wt) - length(coef(wt_mod_2))) # Доверительный интервал new_data$upr <- new_data$fit + t_crit * new_data$se new_data$lwr <- new_data$fit - t_crit * new_data$se
Рисуем график предсказаний
Pl_smoke <- ggplot(new_data, aes(x = age, y = fit)) + geom_ribbon(alpha = 0.2, aes(ymin = lwr, ymax = upr, group = smoke)) + geom_line(aes(colour = smoke)) Pl_smoke
На графике предсказаний можно показать исходные значения
Pl_final <- Pl_smoke +
geom_point(data = wt, aes(x = age, y = bwt, colour = smoke)) +
scale_colour_discrete('Курила ли мать', labels = c('Не курила', 'Курила')) +
labs(x = 'Возраст матери', y = 'Вес новорожденного')
Pl_final
О чем говорят результаты
- Наличие значимого взаимодействия говорит о том, что связь веса новорожденного и возраста матери выглядит по-разному у курящих и некурящих матерей.
- Мы не можем судить о влиянии каждого из предикторов по отдельности (main effects), так как они взаимодействуют.
- При визуализации модели становится видно, что у некурящих матерей есть тенденция с увеличением возраста рожать более крупных детей, у курящих матерей эта связь не выражена.
- У курящих женщин с возрастом более 25 лет есть тенденция рожать более легких детей.
Take home messages
- У дискретных и непрерывных факторов могут быть взаимодействия, поэтому анализ нужно начинать с более полной модели (чтобы можно было протестировать и исключить недостоверные взаимодействия)
Что почитать
Must read paper! Zuur, A.F. and Ieno, E.N., 2016. A protocol for conducting and presenting results of regression‐type analyses. Methods in Ecology and Evolution, 7(6), pp.636-645.
- Кабаков Р.И. R в действии. Анализ и визуализация данных на языке R. М.: ДМК Пресс, 2014
- Zuur, A., Ieno, E.N. and Smith, G.M., 2007. Analyzing ecological data. Springer Science & Business Media.
- Quinn G.P., Keough M.J. 2002. Experimental design and data analysis for biologists
Logan M. 2010. Biostatistical Design and Analysis Using R. A Practical Guide