- Техника подгонки множественных регрессионных моделей
- Проверка условий применимости множественных регрессионных моделей
Вы сможете
- Подобрать множественную линейную модель
- Протестировать ее статистическую значимость и валидность
Множественная регрессия
Пример: Опухоль простаты
Можно ли судить о концентрации простат-специфичного антигена базируясь на знаниях клинических параметров? (Данные Stamey et al., 1989)
Исследовано 97 пациентов, перенесших простатотомию.
Зависимая перменная
lpsa— логарифм концентрации простат-специфичного антигена
Предикторы
lcavol— логарифм объема опухолиlweight— логарифм весаage— возраст пациентаlbph— логарифм степени доброкачественной гиперплазииsvi— поражение семенных пузырьковlcp— логарифм меры поражения капсулыgleason— оценка по шкале Глисонаpgg45— доля оценок 4 и 5 по шкале Глисонаlpsa— логарифм концентрации простат-специфичного антигена
Читаем данные
library(readxl)
prost <- read_excel("data/Prostate.xlsx")
Проверяем, все ли правильно открылось
str(prost)
## Classes 'tbl_df', 'tbl' and 'data.frame': 97 obs. of 9 variables: ## $ lcavol : num -0.58 -0.994 -0.511 -1.204 0.751 ... ## $ lweight: num 2.77 3.32 2.69 3.28 3.43 ... ## $ age : num 50 58 74 58 62 50 64 58 47 63 ... ## $ lbph : num -1.39 -1.39 -1.39 -1.39 -1.39 ... ## $ svi : num 0 0 0 0 0 0 0 0 0 0 ... ## $ lcp : num -1.39 -1.39 -1.39 -1.39 -1.39 ... ## $ gleason: num 6 6 7 6 6 6 6 6 6 6 ... ## $ pgg45 : num 0 0 20 0 0 0 0 0 0 0 ... ## $ lpsa : num -0.431 -0.163 -0.163 -0.163 0.372 ...
Есть ли пропущенные значения?
colSums(is.na(prost))
## lcavol lweight age lbph svi lcp gleason pgg45 ## 0 0 0 0 0 0 0 0 ## lpsa ## 0
Можно ли ответить на вопрос таким методом?
round(cor(prost), 2)
## lcavol lweight age lbph svi lcp gleason pgg45 lpsa ## lcavol 1.00 0.19 0.22 0.03 0.54 0.68 0.43 0.43 0.73 ## lweight 0.19 1.00 0.31 0.43 0.11 0.10 0.00 0.05 0.35 ## age 0.22 0.31 1.00 0.35 0.12 0.13 0.27 0.28 0.17 ## lbph 0.03 0.43 0.35 1.00 -0.09 -0.01 0.08 0.08 0.18 ## svi 0.54 0.11 0.12 -0.09 1.00 0.67 0.32 0.46 0.57 ## lcp 0.68 0.10 0.13 -0.01 0.67 1.00 0.51 0.63 0.55 ## gleason 0.43 0.00 0.27 0.08 0.32 0.51 1.00 0.75 0.37 ## pgg45 0.43 0.05 0.28 0.08 0.46 0.63 0.75 1.00 0.42 ## lpsa 0.73 0.35 0.17 0.18 0.57 0.55 0.37 0.42 1.00
##Так нельзя!
- Обычная корреляция не учитывает, что взаимосвязь между переменными может находиться под контролем других переменых и их взаимодействий.
- Множественные тесты. При тестировании значимости множества коэффициентов корреляции нужно вводить поправку для уровня значимости. Лучше было бы учесть все в одном анализе.
Нам предстоит построить множественную регрессионную модель
\[y_i = \beta_0 + \beta_1x_{1\ i} + \beta_2x_{2\ i} + \beta_3x_{3\ i} + ... + \beta_{p - 1}x_{p - 1\ i} + \varepsilon_i\]
- \(y_i\) - значение зависимой переменной для \(i\)-того наблюдения
- \(\beta_0\) - свободный член (intercept). Значение \(y\) при \(x_1=x_2=x_3=....=x_{p-1}=0\)
- \(\beta_1\) - частный угловой коэффициент для зависимости \(y\) от \(x_1\). Показывает насколько единиц изменяется \(Y\) при изменении \(x_1\) на одну единицу и при условии, что все остальные предикторы не изменяются.
\(\beta_2\), \(\beta_3\), …., \(\beta_p\) - аналогично - \(\varepsilon_i\) - варьирование \(y_i\), не объясняемое данной моделью
Геометрическая интерпретация множественной линейной модели
Для случая с одним предиктором \(y_i = \beta_0 + \beta_1x_i + \varepsilon_i\) — линия регрессии
Геометрическая интерпретация множественной линейной модели
Для случая с двумя предикторами \(y_i = \beta_0 + \beta_1x_{i1} + \beta_2x_{i2} + \varepsilon_i\) — плоскость в трехмерном пространстве
Геометрическая интерпретация множественной линейной модели
Для случая с большим количеством предикторов
\[y_i = \beta_0 + \beta_1x_{1\ i} + \beta_2x_{2\ i} + \beta_3x_{3\ i} + ... + \beta_{p - 1}x_{p - 1\ i} + \varepsilon_i\]
Плоскость в n-мерном пространстве, оси которого образованы значениями предикторов
Исследование данных (Data Exploration)
Задание: Постройте dotplot для всех переменных, оцените присутствие отскакивающих значений
Решение
Явные проблемы — есть сильные корреляции между некоторым предикторами
Задание
- Напишите код, который позволит рассчитать параметры линейной модели, описывающей зависимость lpsa от всех остальных величин (lcavol, lweight, age, lbph, svi, lcp, gleason, pgg45)
Решение
mod1 <- lm(lpsa ~ lcavol + lweight + age + lbph + svi + lcp + gleason + pgg45, data = prost) mod1 <- lm(lpsa ~ . , data = prost) # то же самое summary(mod1)
## ## Call: ## lm(formula = lpsa ~ ., data = prost) ## ## Residuals: ## Min 1Q Median 3Q Max ## -1.733 -0.371 -0.017 0.414 1.638 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 0.66940 1.29638 0.52 0.6069 ## lcavol 0.58702 0.08792 6.68 0.0000000021 *** ## lweight 0.45446 0.17001 2.67 0.0090 ** ## age -0.01964 0.01117 -1.76 0.0823 . ## lbph 0.10705 0.05845 1.83 0.0704 . ## svi 0.76616 0.24431 3.14 0.0023 ** ## lcp -0.10547 0.09101 -1.16 0.2496 ## gleason 0.04514 0.15746 0.29 0.7751 ## pgg45 0.00453 0.00442 1.02 0.3089 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.708 on 88 degrees of freedom ## Multiple R-squared: 0.655, Adjusted R-squared: 0.623 ## F-statistic: 20.9 on 8 and 88 DF, p-value: <2e-16
Проверка валидности модели
Вспомним условия применимости линейных моделей
- Линейная связь между зависимой переменной (\(Y\)) и предикторами (\(X\))
- Независимость значений \(Y\) друг от друга
- Нормальное распределение \(Y\) для каждого уровня значений \(X\)
- Гомогенность дисперсий \(Y\) для каждого уровня значений \(X\)
- Отсутствие коллинеарности предикторов (для множественной регрессии)
Задание
- Проверьте условия применимости модели
Решение
1) График расстояния Кука
- Выбросов нет
mod1_diag <- fortify(mod1) ggplot(mod1_diag, aes(x = 1:nrow(mod1_diag), y = .cooksd)) + geom_bar(stat = "identity")
Решение
2) График остатков от предсказанных значений
- Выбросов нет
- Гетерогенность дисперсии не выявляется
- Есть намек на нелинейность связей
gg_resid <- ggplot(data = mod1_diag, aes(x = .fitted, y = .stdresid)) + geom_point() + geom_hline(yintercept = 0) gg_resid
Решение
3) Графики остатков от предикторов в модели и нет
Решение
3) Код для графиков остатков от предикторов в модели и нет
res_1 <- gg_resid + aes(x = lcavol)
res_2 <- gg_resid + aes(x = lweight)
res_3 <- gg_resid + aes(x = age)
res_4 <- gg_resid + aes(x = lbph)
res_5 <- gg_resid + aes(x = svi)
res_6 <- gg_resid + aes(x = lcp)
res_7 <- gg_resid + aes(x = gleason)
res_8 <- gg_resid + aes(x = pgg45)
library(gridExtra)
grid.arrange(res_1, res_2, res_3, res_4,
res_5, res_6, res_7, res_8, nrow = 2)
Решение
4) Квантильный график остатков
- Отклонения от нормального распределения остатков не выявляются
qqPlot(mod1)
## [1] 39 69
Мультиколлинеарность
Мультиколлинеарность
Мультиколлинеарность — наличие линейной зависимости между независимыми переменными (факторами) регрессионной модели.
При наличии мультиколлинеарности оценки параметров получаются неточными, а значит сложно будет дать интерпретацию влияния предикторов на отклик.
Косвенные признаки мультиколинеарности:
- Большие ошибки оценок параметров
- Большинство параметров модели недостоверно отличаются от нуля, но F критерий говорит, что вся модель значима
Проверка на мультиколлинеарность
- Фактор инфляции дисперсии (Variance inflation factor, VIF)
Как рассчитывается VIF
Мы должны оценить какую долю изменчивости конкретного предиктора могут объяснить другие предикторы (т.е. насколько предикторы независимы)
Для каждого предиктора:
- Строим регрессионную модель данного предиктора от всех остальных \[x_1 = c_0 + c_1x_2 +c_2x_3 + .... + c_{p - 1}x_p\]
- Находим \(R^2\) модели
- Вычисляем фактор инфляции дисперсии \[VIF = \frac{1}{1-R^2}\]
Что делать, если мультиколлинеарность выявлена?
- Можно последовательно удалить из модели избыточные предикторы с VIF > 3 (иногда VIF > 2)
- подбираем модель
- считаем VIF
- удаляем предиктор с самым большим VIF
- повторяем 1-3, пока VIF не станет меньше 2
- Можно заменить исходные предикторы новыми независимыми друг от друга переменными, полученными с помощью метода главных компонент
Проверяем отсутствие мультиколлинеарности
Функция vif() из пакета car
vif(mod1)
## lcavol lweight age lbph svi lcp gleason pgg45 ## 2.05 1.36 1.32 1.38 1.96 3.10 2.47 2.97
В нашей модели сильной мультиколлинеарности нет.
Однако, возможно, что lcp - избыточный предиктор
Удалим из модели избыточный предиктор
mod2 <- update(mod1, ~ . - lcp) vif(mod2)
## lcavol lweight age lbph svi gleason pgg45 ## 1.64 1.36 1.29 1.37 1.62 2.47 2.65
Мультиколлинеарность осталась. Стоит удалить pgg45
mod3 <- update(mod2, ~ . - pgg45) vif(mod3)
## lcavol lweight age lbph svi gleason ## 1.64 1.36 1.28 1.37 1.46 1.34
Теперь мультиколлинеарности нет.
В этой модели осталось много незначимых предикторов
summary(mod3)
## ## Call: ## lm(formula = lpsa ~ lcavol + lweight + age + lbph + svi + gleason, ## data = prost) ## ## Residuals: ## Min 1Q Median 3Q Max ## -1.8084 -0.4056 0.0057 0.4410 1.5954 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 0.2709 1.0957 0.25 0.8053 ## lcavol 0.5421 0.0786 6.90 7e-10 *** ## lweight 0.4533 0.1698 2.67 0.0090 ** ## age -0.0170 0.0110 -1.55 0.1247 ## lbph 0.1069 0.0583 1.83 0.0702 . ## svi 0.6943 0.2110 3.29 0.0014 ** ## gleason 0.1106 0.1159 0.95 0.3425 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.708 on 90 degrees of freedom ## Multiple R-squared: 0.648, Adjusted R-squared: 0.624 ## F-statistic: 27.6 on 6 and 90 DF, p-value: <2e-16
Что дальше?
Два варианта действий:
- Оставить все как есть. Если значение коэффициента при предикторе не отличается значимо от нуля, значит, этот предиктор не влияет на объем опухоли. ПРОВЕРЬТЕ ВЫПОЛНЕНИЕ УСЛОВИЙ ПРИМЕНИМОСТИ ДЛЯ ФИНАЛЬНОЙ МОДЕЛИ!
- Провести пошаговый подбор оптимальной модели (об этом позднее)
Пока оставим все как есть и попытаемся выяснить, какие предикторы влияют сильнее всего
Сравнение вклада предикторов
Какой из предиктов оказывает наиболее сильное влияние?
Для ответа на этот вопрос надо “уравнять” шкалы, всех предикторов, то есть стандартизовать их.
Коэффициенты при стандартизованных предикторах покажут, насколько сильно меняется отклик при изменении предиктора на одно стандартное отклонение.
Для стандартизации используем функцию scale()
mod3_scaled <- lm(lpsa ~ scale(lcavol) + scale(lweight) + scale(age) +
scale(lbph) + scale(svi) + scale(lcp) + scale(gleason) + scale(pgg45),
data = prost)
Какой из предиктов оказывает наиболее сильное влияние?
summary(mod3_scaled)
## ## Call: ## lm(formula = lpsa ~ scale(lcavol) + scale(lweight) + scale(age) + ## scale(lbph) + scale(svi) + scale(lcp) + scale(gleason) + ## scale(pgg45), data = prost) ## ## Residuals: ## Min 1Q Median 3Q Max ## -1.733 -0.371 -0.017 0.414 1.638 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 2.4784 0.0719 34.46 < 2e-16 *** ## scale(lcavol) 0.6919 0.1036 6.68 0.0000000021 *** ## scale(lweight) 0.2257 0.0844 2.67 0.0090 ** ## scale(age) -0.1462 0.0832 -1.76 0.0823 . ## scale(lbph) 0.1553 0.0848 1.83 0.0704 . ## scale(svi) 0.3172 0.1011 3.14 0.0023 ** ## scale(lcp) -0.1475 0.1273 -1.16 0.2496 ## scale(gleason) 0.0326 0.1137 0.29 0.7751 ## scale(pgg45) 0.1276 0.1247 1.02 0.3089 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.708 on 88 degrees of freedom ## Multiple R-squared: 0.655, Adjusted R-squared: 0.623 ## F-statistic: 20.9 on 8 and 88 DF, p-value: <2e-16
Визуализация модели
# Создаем искусственный датафрейм, где будет визуализироваться самый важный предиктор, а остальные будут рассматриваться, как средние
MyData <- data.frame(lcavol = seq(min(prost$lcavol), max(prost$lcavol), 0.5),
lweight = mean(prost$lweight),
age = mean(prost$age),
lbph = mean(prost$lbph),
svi = mean(prost$svi),
lcp = mean(prost$lcp),
gleason = mean(prost$gleason),
pgg45 = mean(prost$pgg45) )
MyData$Predict <- predict(mod3, newdata = MyData)
ggplot(MyData, aes(x = lcavol, y = Predict)) + geom_line() + geom_point(data = prost, aes(x = lcavol, y = lpsa))
Визуализация модели с использовнием двух пердикторов
MyData <- expand.grid(lcavol = seq(min(prost$lcavol), max(prost$lcavol), 0.5),
lweight = mean(prost$lweight),
age = mean(prost$age),
lbph = mean(prost$lbph),
svi = unique(prost$svi),
lcp = mean(prost$lcp),
gleason = mean(prost$gleason),
pgg45 = mean(prost$pgg45) )
MyData$Predict <- predict(mod3, newdata = MyData)
ggplot(MyData, aes(x = lcavol, y = Predict, color = factor(svi))) + geom_line() + geom_point(data = prost, aes(x = lcavol, y = lpsa)) + facet_wrap(~svi)
Визуализация модели с использовнием трех пердикторов
MyData <- expand.grid(lcavol = seq(min(prost$lcavol), max(prost$lcavol), 0.5),
lweight = seq(min(prost$lweight), max(prost$lweight), 0.1),
age = mean(prost$age),
lbph = mean(prost$lbph),
svi = unique(prost$svi),
lcp = mean(prost$lcp),
gleason = mean(prost$gleason),
pgg45 = mean(prost$pgg45) )
MyData$Predict <- predict(mod3, newdata = MyData)
ggplot(MyData, aes(x = lcavol, y = Predict, color = lweight, group = lweight)) + geom_line() + geom_point(data = prost, aes(x = lcavol, y = lpsa)) + facet_wrap(~svi) + scale_color_continuous(high = "red", low = "yellow") + ylab("lpsa")
Вопрос
Какая доля суммарной дисперсии зависимой переменной описывается данной моделью?
Adjusted R-squared - скорректированный коэффициет детерминации
Применяется если необходимо сравнить две модели с разным количеством параметров
\[ R^2_{adj} = 1- (1-R^2)\frac{n-1}{n-p}\]
\(p\) - количество параметров в модели
Вводится штраф за каждый новый параметр
Какой из предиктов оказывает наиболее сильное влияние?
- Сильнее всего c концентрацией простат-специфичного антигена связан логарифм объема опухоли (
lcavol) и оценка поражения семенных пузырьков (svi) - При изменении логарифма объема опухоли на 1SD логарифм концентрации антигена изменяется на 0.69
- При изменении логарифма оценки поражения семенных пузырьков на 1SD, логарифм концентрации антигена изменяется на 0.32
Include or Don’t include? That is the question…
В рассмотренной модели мы не пытались выяснить есть ли взаимодействия предикторов. То есть мы изначально заложили в модель идею, что влияние на зависимую переменную каждого из предикторов не зависит от других предикторов.
Include or Don’t include? That is the question…
Вопрос о включении в модель взаимодействия предикторов совсем непростой
Существует несколько подходов:
- Не включать взаимодействия в модель. Но если при валидации модели в остатках появляется явный паттерн, то это может быть следствием наличия взаимодействия предикторов
- Основываясь на априорных знаниях свойств объектов включить только те взаимоотношения, которые имеют биологический смысл, либо взаимодействия с наиболее важными переменными (теми, ради которых была затеяна работа)
- Включать в модель все взаимодействия, потом пошагово выбросить недостоверные (Model selection - об этом позднее)
Include or Don’t include? That is the question…
Включать в анализ и обсуждать все взаимодействия “дорого” (неудобно):
- Взаимодействия высоких порядков сложно интерпретировать
- Каждое взаимодействие — это коэффициент в модели, или несколько, если это взаимодействие с дискретной переменной. Чтобы подобрать модель нужно много данных — по 20-40 наблюдений в расчете на каждый коэффициент.
Summary
- При построении множественной регрессии важно, помимо других условий, проверить модель на наличие мультиколлинеарности
- Если модель построена на основе стандартизированных значений предикторов, то можно сравнивать влияние этих предикторов
- В модель можно (а иногда и нужно) включать взаимодействия предикторов
Что почитать
- Кабаков Р.И. R в действии. Анализ и визуализация данных на языке R. М.: ДМК Пресс, 2014.
- Quinn G.P., Keough M.J. (2002) Experimental design and data analysis for biologists, pp. 92-98, 111-130
- Diez D. M., Barr C. D., Cetinkaya-Rundel M. (2014) Open Intro to Statistics., pp. 354-367.
- Logan M. (2010) Biostatistical Design and Analysis Using R. A Practical Guide, pp. 170-173, 208-211
- Zuur, A.F. et al. 2009. Mixed effects models and extensions in ecology with R. - Statistics for biology and health. Springer, New York, NY. pp. 538-552.