- Линейные модели c непрерывными и дискретными предикторами
Мы рассмотрим
Дискретные предикторы, или факторы
Дискретные предикторы, или факторы
Это перменные, которые принимают несколько кчачественно различных значений, между которыми нет переходов.
Разные занчения одного и того же фактора называются уровнями фактора. Уровней может быть два и более.
Примеры:
Если мы изучаем связь какой-либо величины (например, роста) с полом, то переменная Sex может принимать два значения: "F" и "M".
F и M - это уровни фактора Sex.
Если мы изучаем связь экспрессии некоторого гена с введением какого-то вещества, то переменная Treatment может принимать три значения: "Control", "Drag_1", "Drag_2"
"Control", "Drag_1", "Drag_2" - Это уровни фактора Treatment
Геометрическая интерпретация моделей с непрерывным и дискретным предиктором
В моделях с одним непрерывным предиктоом
Это одна линия, которая описывается уравнением \(y = b_0 + b_1x\)
В моделях с одним непрерывным предиктоом и одним дискретным предиктором
В простейшем случае, когда есть только две градации дискретного фактора, - это две линии
Если для каждой градации дискретного фактора характер связи между зависимой переменной \(y\) и непрерывным предиктором \(x\) один и тот же (нет статистически значимого взаимодействия), то это две параллельные прямые.
Эти прямые отличаются только значениями Intercept
Глистогонные и рост коз
Как связан прирост массы коз с начальным весом животного и интенсивностью профилактики паразитарных заболеваний?
Goat by Jennifer C. on Flickr
Treatment- обработка от глистов (стандартная, интенсивная)Weightgain- привес, кгInitial.wt- начальный вес, кг
Пример из библиотеки данных http://www.statlab.uni-heidelberg.de/data/ancova/goats.story.html
Читаем данные и знакомимся с ними
library(readxl)
goat <- read_excel("data/goats.xlsx", sheet = 1)
head(goat)
## # A tibble: 6 <U+00D7> 3 ## Treatment Weightgain `Initial wt` ## <chr> <dbl> <dbl> ## 1 standard 5 21 ## 2 standard 3 24 ## 3 standard 8 21 ## 4 standard 7 22 ## 5 standard 6 23 ## 6 standard 4 26
str(goat)
## Classes 'tbl_df', 'tbl' and 'data.frame': 40 obs. of 3 variables: ## $ Treatment : chr "standard" "standard" "standard" "standard" ... ## $ Weightgain: num 5 3 8 7 6 4 8 6 7 5 ... ## $ Initial wt: num 21 24 21 22 23 26 22 23 24 20 ...
Читаем данные и знакомимся с ними
colSums(is.na(goat))
## Treatment Weightgain Initial wt ## 0 0 0
# переименуем переменные для краткости
colnames(goat) <- c("Treatment", "Wt", "Init")
# объемы выборок
table(goat$Treatment)
## ## intensive standard ## 20 20
goat$Treatment <- factor(goat$Treatment)
Есть ли выбросы?
library(ggplot2) gg_dot <- ggplot(goat, aes(y = 1:nrow(goat))) + geom_point() gg_dot + aes(x = Wt)
gg_dot + aes(x = Init)
Строим модель
MG <- lm(Wt ~ Init + Treatment, data = goat)
В этой модели мы молчаливо считаем, что характер связи прироста коз с начальным весом будет одинаковым (нет взаимодействия предикторов). Но! Это надо специально проверять (об этом далее)
Проверяем условия применимости
Пока не смотрим в summary()!
Нет ли колинеарности между начальным весом и тритментом
vif(MG)
## Init Treatment ## 1 1
ggplot(goat, aes(x = Treatment, y = Init)) + geom_boxplot()
Выраженной коллинеарности нет
Проверяем условия применимости
Все хорошо
Нормальнсть распределения остатков
График модели
Результаты
summary(MG)
## ## Call: ## lm(formula = Wt ~ Init + Treatment, data = goat) ## ## Residuals: ## Min 1Q Median 3Q Max ## -2.972 -1.242 -0.034 0.988 3.223 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 14.9666 1.7526 8.54 0.00000000028 *** ## Init -0.3514 0.0742 -4.73 0.00003207789 *** ## Treatmentstandard -1.2649 0.5117 -2.47 0.018 * ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 1.62 on 37 degrees of freedom ## Multiple R-squared: 0.438, Adjusted R-squared: 0.408 ## F-statistic: 14.4 on 2 and 37 DF, p-value: 0.0000233
В summary() появиилась строка Treatmentstandard.
Почему не Treatment?
Информация в summary() необходима для построения уравнения модели
У нас две линии
Первая: \[Wt = 14.97 - 0.3514 Init\]
Эта прямая описывает свяь между Wt и Init для той градации фактора Treatment, которая взята за базовый уровень. По умолчанию в качестве базового уровня берется первая по алфавиту градация дискретного фактора. В данном случае - это градация intensive
Вторая линия описывает связь между Wt и Init для второй градации фактора Treatment. уравнение этой прямой будет таким:
\[ Wt = 14.97 -0.3514 Init -1.2649 \\ Wt = 13.7 -0.3514 Init \]
Смысл коэффициентов при дискретных предикторах
Смысл 1.
Коэффициенты при дискретных предикторах показывают на сколько единиц смещается Intercept у данного (указанного в summsry()) уровня по отношению к базовому уровню
Смысл 2.
На сколько единиц отличается среднее значение зависимой переменной для данного (указанного в summsry()) уровня от среднего значения зависимой переменной для базового уровня.
Уровень значимости в summary() говорит о статистической значимости этих отличий (при попарном сравнении; если уровней больше двух, то помните про проблемы множественности сравнений)
Общее уравнение модели
Уравнение преобретает такой вид \[ Wt = 14.97 - 0.35Init - 1.26I_{standard} \] В этом уравнении появилась перменная \(I\) - это переменная-селектор (dummy variable, переменная-болванка)
\(I = 0\) если мы рассматриваем базовый уровень дискретного предиктора. В данном случае уровень Treatment = "initial"
\(I = 1\) если мы рассматриваем другой уровень дискретного предиктора. В данном случае уровень Treatment = "standard"
Меняем базовый уровень
Это чисто формальная процедура от которой ничего не измеяется по сути, но это иногда необходимо для более удобной визуализации
goat$Treatment <- relevel(goat$Treatment, ref = "standard") levels(goat$Treatment)
## [1] "standard" "intensive"
Строим новую модель
MG1 <- lm(Wt ~ Init + Treatment, data = goat)
Результаты
summary(MG1)
## ## Call: ## lm(formula = Wt ~ Init + Treatment, data = goat) ## ## Residuals: ## Min 1Q Median 3Q Max ## -2.972 -1.242 -0.034 0.988 3.223 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 13.7017 1.7599 7.79 0.0000000026 *** ## Init -0.3514 0.0742 -4.73 0.0000320779 *** ## Treatmentintensive 1.2649 0.5117 2.47 0.018 * ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 1.62 on 37 degrees of freedom ## Multiple R-squared: 0.438, Adjusted R-squared: 0.408 ## F-statistic: 14.4 on 2 and 37 DF, p-value: 0.0000233
Что изменилось?
Уравнения
Ощее уравнение модели
\[
Wt = 13.7017 -0.3514 Init + 1.2649 I_{intensive}
\] Постройте уравнения для каждого из уровней фактора Treatment….
Как охарактеризовать влияние двух предикторов в целом
library(car) Anova(MG, type = 3)
## Anova Table (Type III tests) ## ## Response: Wt ## Sum Sq Df F value Pr(>F) ## (Intercept) 190.9 1 72.93 0.00000000028 *** ## Init 58.6 1 22.40 0.00003207789 *** ## Treatment 16.0 1 6.11 0.018 * ## Residuals 96.9 37 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Вывод:
Привес коз при интенсивной обработке от глистов на 1.26 кг выше, чем при стандартном методе обработки (ANOVA, p =0.018). Привес коз демонстрирует отрицательную связь с начальным весом животного (ANOVA, p < 0.01).
Применение моделей с дискретными предикторами
Когда это нужно?
- Для демонстрации связи с непрерывными и дискретными прдикторами одновременно в одной модели.
- Если в фокусе нашего внимания влияние дискретного предиктора, но это влияние зависит от дополнительной перменной - ковариаты (такой анализ называется ANCOVA)
ANCOVA
Влияет ли стаж работы на предприятиях, вырабатывающих кадмий, на жизненнй объем легких?
Пример взят из книги: P. Armitage and G. Berry (1987), Statistical Methods in Medical Research, 2nd ed., Blackwell, p.286.
Данные представлены в пакете ISwR
Переменные:
group - Группа 1: Более 10 лет в отрасли; Группа 2 - менее 10 лет; Группа 3 - не подвергались воздействию.
age - возраст
vital.capacity - объем легких (л).
Загружаем данные
vit <- read.table("data/vitcap2.csv", header = TRUE, sep = ";")
Немного преобразуем исходный датасет
vit$Group [vit$group == 1] <- "Long exposed" vit$Group [vit$group == 2] <- "Short exposed" vit$Group [vit$group == 3] <- "Not exposed"
Меняем порядок уровней
vit$Group <- factor(vit$Group, levels = c("Not exposed", "Short exposed", "Long exposed"))
levels(vit$Group)
## [1] "Not exposed" "Short exposed" "Long exposed"
В фокусе исследования дискретная переменная
Нас интересует зависит ли объем легких (vital.capacity) от стажа работы (Group)
Можно ли построить модель такого вида?
\[ vital.capacity = \huge {f} \normalsize (Group) \] то есть модель должна быть основана только на дискретных предикторах.
Моделями с дискретными предикторами занимается классический дисперсионный анализ
Модель с дискретным предиктором
M1 <- lm(vital.capacity ~ Group, data = vit) library(car) Anova(M1, type = 3)
## Anova Table (Type III tests) ## ## Response: vital.capacity ## Sum Sq Df F value Pr(>F) ## (Intercept) 876 1 1580.60 <2e-16 *** ## Group 3 2 2.48 0.09 . ## Residuals 45 81 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Геометрическая интерпретация модели с дискретным предиктором
Это будет график, отражающий средние значения зависимой переменной, вычисленные для каждой градации дискретного фактора
Геометрическая интерпретация модели с дискретным предиктором
MyData <- data.frame(Group = levels(vit$Group))
MyData$Group <- factor(MyData$Group, levels = c("Not exposed", "Short exposed", "Long exposed"))
MyData$Predicted <- predict(M1, newdata = MyData, se.fit = TRUE)$fit
MyData$SE <- predict(M1, newdata = MyData, se.fit = TRUE)$se.fit
Геометрическая интерпретация модели с дискретным предиктором
library(ggplot2) ggplot(MyData, aes(x = Group, y = Predicted)) + geom_bar(stat = "identity", aes(fill = Group)) + geom_errorbar(aes(ymin = Predicted - SE, ymax = Predicted + SE), width = 0.2)
Смотрим на результаты
summary(M1)
## ## Call: ## lm(formula = vital.capacity ~ Group, data = vit) ## ## Residuals: ## Min 1Q Median 3Q Max ## -1.7718 -0.4520 0.0981 0.5130 1.5708 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 4.46205 0.11223 39.76 <2e-16 *** ## GroupShort exposed 0.00974 0.17997 0.05 0.957 ## GroupLong exposed -0.51288 0.24245 -2.12 0.037 * ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.744 on 81 degrees of freedom ## Multiple R-squared: 0.0577, Adjusted R-squared: 0.0344 ## F-statistic: 2.48 on 2 and 81 DF, p-value: 0.0902
Куда делась одна градация фактора?
- В качестве базового уровня предиктора
Groupвзято значение "Not exposed"
Запись модели с дискретными предикторами
Общая формула для модели \[ vital.capacity = 4.46205 + 0.00974 \cdot I_{Short exposed} -0.51288 \cdot I_{Long exposed} \] \(I\) - это перменные-селекторы (dummy variable)
Запись модели с дискретными предикторами
Если мы хотим вычислить предсказание модели (значение объема легких) для любого индивида, имеющего Group = "Short exposed"
\[ vital.capacity = 4.46205 + 0.00974 \cdot 1 -0.51288 \cdot 0 = \\ = 4.46205 + 0.00974 = 4.47 \]
Запись модели с дискретными предикторами
Если мы хотим вычислить предсказание модели (значение объема легких) для индавида, имеющего Group = "Long exposed"
\[ vital.capacity = 4.46205 + 0.00974 \cdot 0 -0.51288 \cdot 1 = \\ = 4.46205 -0.51288 = 3.95 \]
Запись модели с дискретными предикторами
Если мы хотим вычислить предсказание модели (значение объема легких) для индавида, имеющего Group = "Not exposed"
\[ vital.capacity = 4.46205 + 0.00974 \cdot 0 -0.51288 \cdot 0 = \\ = 4.46205 \]
Это базовый уровень. Для него предсказанное значение равно Intercept
Задание
- Измените базовый уровень переменной
Groupна "Long exposed" - Постройте модель, аналогичную
M1 - Вычислите предсказанные моделью значения для каждой градации фактора
Group
Вернемся к результатам
Можно ли доверять полученным результатам?
M1_diag <- fortify(M1) qplot(vit$age, M1_diag$.stdresid) + geom_smooth(method = "lm")
Очевидный паттерн в остатках!
Необходимо включать еще одну переменную - ковариату
Analysis of covariance (ANCOVA)
Меняем модель
M3 <- lm(vital.capacity ~ Group + age , data = vit)
Диагностика модели
M3_diag <- fortify(M3) qplot(vit$age, M3_diag$.stdresid) + geom_smooth(method = "lm")
Паттерн исчез!
Результаты
summary(M3)
## ## Call: ## lm(formula = vital.capacity ~ Group + age, data = vit) ## ## Residuals: ## Min 1Q Median 3Q Max ## -1.3805 -0.3805 0.0132 0.3791 1.3705 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 6.04492 0.26802 22.55 < 2e-16 *** ## GroupShort exposed -0.07020 0.14867 -0.47 0.64 ## GroupLong exposed -0.11693 0.20924 -0.56 0.58 ## age -0.03978 0.00632 -6.29 0.000000016 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.613 on 80 degrees of freedom ## Multiple R-squared: 0.37, Adjusted R-squared: 0.346 ## F-statistic: 15.6 on 3 and 80 DF, p-value: 0.0000000432
Результаты
anova(M3)
## Analysis of Variance Table ## ## Response: vital.capacity ## Df Sum Sq Mean Sq F value Pr(>F) ## Group 2 2.75 1.37 3.66 0.03 * ## age 1 14.86 14.86 39.58 0.000000016 *** ## Residuals 80 30.03 0.38 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Противоречие с результатами summary()!
Поменяем порядок предикторов
anova(lm(formula = vital.capacity ~ age + Group, data = vit))
## Analysis of Variance Table ## ## Response: vital.capacity ## Df Sum Sq Mean Sq F value Pr(>F) ## age 1 17.44 17.44 46.47 0.0000000016 *** ## Group 2 0.16 0.08 0.22 0.81 ## Residuals 80 30.03 0.38 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
` Результат тестирования зависит от порядка предикторов. Почему?
Два варианта тестирования
- Последовательное
- Иерархическое
Вариант 1. Последовательное тестирование (SS type I)
Факторы тестируются в порядке включения в модель. Результат тестирования зависит от порядка включения.
anova(lm(formula = vital.capacity ~ Group + age, data = vit))
## Analysis of Variance Table ## ## Response: vital.capacity ## Df Sum Sq Mean Sq F value Pr(>F) ## Group 2 2.75 1.37 3.66 0.03 * ## age 1 14.86 14.86 39.58 0.000000016 *** ## Residuals 80 30.03 0.38 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
- Тот предиктор, который интересует исследователия надо поставить на второе место. Влияние этого фактора тестируется после удаления влияния ковариат.
anova(lm(formula = vital.capacity ~ age + Group, data = vit))
## Analysis of Variance Table ## ## Response: vital.capacity ## Df Sum Sq Mean Sq F value Pr(>F) ## age 1 17.44 17.44 46.47 0.0000000016 *** ## Group 2 0.16 0.08 0.22 0.81 ## Residuals 80 30.03 0.38 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Вариант 2. Иерархическое тестирование (SS type III)
Каждый из факторов по отношению к модели только без него, но со всеми остальными. Нужно, если много факторов и выборки разного размера. Тогда результат не будет зависеть от порядка включения факторов в модель.
library(car) Anova(M3, type = 3)
## Anova Table (Type III tests) ## ## Response: vital.capacity ## Sum Sq Df F value Pr(>F) ## (Intercept) 191.0 1 508.66 < 2e-16 *** ## Group 0.2 2 0.22 0.81 ## age 14.9 1 39.58 0.000000016 *** ## Residuals 30.0 80 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Take home messages
- Нет разницы между моделями с дискретными и непрерывными предикторами
- При анализе влияния дискретного фактора часто необходимо учитывать влияние ковариаты.
- Пр тестировании гипотезы надо помнить о существовании двух типов тестирования: последовательное vs иерархическое.
Что почитать
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