• Линейные модели c непрерывными и дискретными предикторами

Мы уже умеем описывать реальность с помощью моделей.
Например, модель движения автомобиля выглядит так:

\[S=Vt\] Можем ли мы построить модель движения машины, у которой может быть два разных водителя, предпочитающих ездить с разными скоростями?

Путь, проделанный машиной за некоторое время \(t\), будет зависеть от того, кто сидел за рулем.

Это будет множественная модель, включающая помимо непрерывного предиктора \(t\) еще и дискретный предиктор “Водитель”, Driver.

Дискретные предикторы по традиции называют “факторами”.

У фактора Driver есть две градации (levels): Driver_1 и Driver_2

График модели будет разным для двух градаций дискретного фактора.

От значения фактора Driver зависит угловой коэффициент модели.

Модель можно представить в виде формулы, но старая формула \(S=Vt\) не годится!

Новая формула:

\[S = (V_1 + b_2 \cdot I_{Driver_2})t\]

• \(V_1\) – Скорость водителя №1
  
• \(t\) – Время
• \(I_{Driver_2}\) – Переменная-индикатор
• \(b_2\) – Поправочный коэффициент, показывающий на сколько скорость второго водителя отличается от скорости первого.

Ту же самую модель можно записать в привычных обозначениях, использованных ранее:

\[y = b_1 x + b_2 \cdot I_{Driver_2} \cdot x\]

• \(y\) — зависимая величина (путь)
  
• \(x\) — предиктор (время)
• \(b_1\) — коэффициент, равный скорости Водителя_1
• \(b_2\) — коэффициент, показывающий насколько скорость второго водителя отличается от скорости первого
• \(I_{Driver_2}\) — переменная-индикатор

\[y = b_1 x + b_2 \cdot I_{Driver_2}\cdot x\]

\(I_{Driver_2}\) — переменная-индикатор:

• \(I_{Driver_2} = 0\) если рассматриваем водителя №1
• \(I_{Driver_2} = 1\) если рассматриваем водителя №2

Переменные-индикаторы, или переменные-болванки (dummy variable), “переключают” модель между разными уровнями дискретного фактора.

У дискретных факторов всегда есть базовый уровень.

Базовый уровень — это градация дискретного предиктора, относительно которой вычисляются коэффициенты модели.

В нашей модели в качестве базового уровня фактора Driver выбрана градация Driver_1.

\[y = b_1 x + b_2 \cdot I_{Driver_2}\cdot x\]

\(b_2\) — поправочный коэффициент, который показывает, на сколько единиц (км/ч) отличается скорость Driver_2 от скорости на базовом уровне (Driver_1).

Пусть скорость, которую развивает второй водитель, будет на 20 км/ч выше (\(b_2 = 20\)), чем скорость первого водителя (\(b_1\)).

Модель для Driver_1 (базовый уровень):

\[ I_{Driver_2}=0 \qquad\qquad y = b_1 x + 20 \cdot 0 \cdot x = b_1 x \qquad\quad \]

Модель для Driver_2

\[ I_{Driver_2} = 1 \qquad\qquad y = b_1 x + 20 \cdot 1\cdot x = (b_1 + 20)x \]

Предположим, вы бронируете номер в гостинице, где могут быть номера с двумя кроватями (Twin room), с одной большой кроватью (Double room) и семейные номера (Family room).

Модель затрат на оплату проживания в зависимости от времени должна включать дискретный предиктор Room_type с тремя градациями: Twin, Double, Family.

Все три графика параллельны, но одни выше, другие – ниже.

\(y = b_1 x + b_2 \cdot I_{Double} + b_3 \cdot I_{Family}\)

• \(y\) — зависимая переменная (потраченная сумма)
  
• \(x\) — непрерывный предиктор (время)
  
• \(b_1\) — цена одной ночи в номере Twin (базовый уровень)
  
• \(b_2\) — на сколько единиц цена номера Double отличается от базового уровня
  
• \(b_3\) — на сколько единиц цена номера Family отличается от базового уровня
  
• \(I_{Double}, I_{Family}\) — переменные-индикаторы

Общее уравнение:

\[\qquad y = b_1 x_i + b_2 \cdot I_{Double} + b_3 \cdot I_{Family}\]

• Уравнение для проживания в номере типа Twin (базовый уровень):

\[y = b_1 x + b_2 \cdot 0 + b_3 \cdot 0 = b_1 \cdot x\]

• Уравнение для проживания в номере типа Double:

\[y = b_1 x + b_2 \cdot 1 + b_3 \cdot 0 = b_1 x + b_2\]

• Уравнение для проживания в номере типа Family:

\[y = b_1 x + b_2 \cdot 0 + b_3 \cdot 1 = b_1 x + b_3\]

Если характер связи зависимой переменной с непрерывным предиктором одинаков для разных градаций дискретного фактора (прямые параллельны), то говорят, что между предикторами отсутствует взаимодействие.

Изменение размера счета в банке зависит от уровня доходов человека, от уровня его трат и т.п.

Сумма на банковском счете у людей, занимающих разные должности, разная.

Поэтому модель, описывающая изменение суммы на счете, должна включать дискретный фактор Position (Должность).

Графики не параллельны – у них различаются как значения свободного члена (intercept), так и угловые коэффициенты (slope).

В такой ситуации говорят, что между дискретным фактором (Position, должность) и непрерывным предиктором (X, время) существует взаимодействие.

\[ y = b_0 + b_1x + b_2 \cdot I_{Position_2} + b_3 \cdot I_{Position_2} \cdot x \]

• \(y\) — зависимая переменная (сумма на банковском счете)
• \(x\) — предиктор (время)
• \(b_0\) — состояние счета работника на должности №1 в начальный момент времени (\(X=0\))
• \(b_1\) — коэффициент, равный доходам (руб/мес) сотрудника на должности №1
• \(I_{Position_2}\) — переменная-индикатор
• \(b_2\) — коэффициент, показывающий, на сколько рублей банковский счет на должности №2 отличается от счета сотрудника на должности №1 в начальный момент времени (\(x = 0\))
• \(b_3\) — на сколько единиц (руб/мес) доходы сотрудника на должности №2 отличаются от доходов сотрудника на должности №1

Общее уравнение:

\[y = b_0 + b_1x + b_2 \cdot I_{Position_2} + b_3 \cdot I_{Position_2} \cdot x\]

• Уравнение для должности №1 (базовый уровень):

\[\begin{aligned}y &= b_0 + b_1x + b_2 \cdot 0 + b_3 \cdot 0 \cdot x = \\ &= b_0 + b_1x\end{aligned}\]

• Уравнение для должности №2:

\[\begin{aligned} y &= b_0 + b_1x + b_2 \cdot 1 + b_3 \cdot 1 \cdot x = \\ &= b_0 + b_1x + b_2 + b_3 \cdot x = \\ &= (b_0 + b_2) + (b_1 + b_3)x \end{aligned}\]

У нас еще будет детальный разговор о взаимодействии дискретных и непрерывных предикторов.

Мы рассмотрели простейшие детерминистские модели.

В регрессионном анализе все аналогично!

Но есть две особенности:

• В моделях добавляется случайная часть. Они становятся стохастическими.
  
• Мы заранее не можем сказать, есть ли в изучаемой системе взаимодействие предикторов или его нет. Выявление взаимодействия между предикторами — специальная задача.

Как связан прирост массы коз с начальным весом животного и интенсивностью профилактики паразитарных заболеваний?

library(readxl)
goat <- read_excel("data/goats.xlsx", sheet = 1)
head(goat)

## # A tibble: 6 × 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)

## tibble [40 × 3] (S3: tbl_df/tbl/data.frame)
##  $ Treatment : chr [1:40] "standard" "standard" "standard" "standard" ...
##  $ Weightgain: num [1:40] 5 3 8 7 6 4 8 6 7 5 ...
##  $ Initial wt: num [1:40] 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", "Stw")
# объемы выборок
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 = Stw)

Нас интересует, как зависит прирост массы коз (Wt) от непрерывного предиктора (Stw) и дискретного фактора Treatment.

Но! Мы не можем сразу решить есть ли в системе взаимодействие или его нет. Поэтому мы должны рассмотреть модель, в которой взаимодействие будет включено.

Зависимость в генеральной совокупности:

\[ Wt_i = \beta_0 + \beta_1Stw_i + \beta_2 \cdot I_{standart\,i} + \beta_3 \cdot I_{standart\,i} \cdot Stw_i + \varepsilon_i, \quad \varepsilon \sim N(0, \sigma) \]

Модель, подобранная по выборке:

\[ Wt_i = b_0 + b_1Stw_i + b_2 \cdot I_{standart\,i} + b_3 \cdot I_{standart\,i} \cdot Stw_i + e_i \]

Предсказанные значения:

\[ \widehat{Wt}_i = b_0 + b_1Stw_i + b_2 \cdot I_{standart\,i} + b_3 \cdot I_{standart\,i} \cdot Stw_i \]

Мы не знаем, взаимодействуют ли дискретный фактор Treatment и непрерывный предиктор Stw.

Возможна ситуация, при которой худые и жирные козы будут по-разному реагировать на разные типы антигельминтной терапии.

В полной модели мы должны учесть как влияние самих предикторов, так и влияние их взаимодействия.

Mod_goat_full <- lm(Wt ~ Stw + Treatment + Stw:Treatment, data = goat)

# Аналогичная, но более короткая запись
Mod_goat_full <- lm(Wt ~ Stw * Treatment, data = goat)

Базовым уровнем фактора будет значение Treatment = intensive. По умолчанию в R используется первое по алфавиту значение.

drop1(Mod_goat_full, test = 'F')

## Single term deletions
## 
## Model:
## Wt ~ Stw * Treatment
##               Df Sum of Sq  RSS  AIC F value Pr(>F)
## <none>                     96.5 43.2               
## Stw:Treatment  1     0.342 96.9 41.4    0.13   0.72

Значимого взаимодействия не наблюдается!

Формула сокращенной модели

\[ Wt_i = b_0 + b_1Stw_i + b_2 \cdot I_{standart\,i} + e_i \]

\[ \widehat{Wt}_i = b_0 + b_1Stw_i + b_2 \cdot I_{standart\,i} \]

Можно просто переписать формулу.

# Можно просто переписать формулу
Mod_goat_reduced <- lm(Wt ~ Stw + Treatment, data = goat)

Но есть более удобный способ.

Mod_goat_reduced <- update(Mod_goat_full, . ~ . - Stw:Treatment)

Пока не смотрим в summary()!

library(car)
vif(Mod_goat_reduced)

##       Stw Treatment 
##         1         1

ggplot(goat, aes(x = Treatment, y = Stw)) + 
  geom_boxplot()

Выраженной коллинеарности нет.

MG_diag <- fortify(Mod_goat_reduced)
  Pl_diag1 <- ggplot(MG_diag, aes(x = 1:nrow(MG_diag), y = .cooksd)) +
    geom_bar(stat = 'identity')
  
  Pl_diag2 <- ggplot(MG_diag, aes(x = .fitted, y = .stdresid)) +
    geom_point() + geom_hline(yintercept = 0)
  
  Pl_diag3 <- ggplot(MG_diag, aes(x = Stw, y = .stdresid)) +
    geom_point() + geom_hline(yintercept = 0)
  
  Pl_diag4 <- ggplot(MG_diag, aes(x = Treatment, y = .stdresid)) +
    geom_boxplot()

library(cowplot)
plot_grid(Pl_diag1, Pl_diag2, Pl_diag3, Pl_diag4, nrow = 2)

Влиятельных наблюдений нет.
Гетерогенности дисперсий нет.

qqPlot(Mod_goat_reduced, id = FALSE) # функция из пакета car

gg_goat <- ggplot(data = goat, aes(y = Wt, x = Stw, colour = Treatment)) +
  geom_point()  +
  labs(x = 'Начальный вес, кг', y = 'Привес, кг') +
  scale_colour_discrete('Способ обработки',
                  breaks = c('intensive', 'standard'),
                  labels = c('Интенсивный', 'Стандартный'))
gg_goat

library(dplyr)
new_data <- goat %>% group_by(Treatment)%>%
  do(data.frame(Stw = seq(min(.$Stw), max(.$Stw), length.out = 10))) 
new_data$fit = predict(Mod_goat_reduced, newdata = new_data)

gg_goat + geom_line(data = new_data, aes(x = Stw, y = fit, colour = Treatment)) 

# Находим вариационно-ковариационную матрицу для модели Mod_goat_reduced
VC <- vcov(Mod_goat_reduced)

# Модельная матрица
X <- model.matrix(~ Stw + Treatment, data = new_data)

# Стандартные ошибки
new_data$se <- sqrt(diag(X %*% VC %*% t(X)))

# t для расчета 95% доверительного интервала 
t_crit <- qt(0.975, df = nrow(goat) - length(coef(Mod_goat_reduced))) 
# t для расчета доверит. интервала 

# Границы доверительных интервалов
new_data$upr <- new_data$fit + t_crit * new_data$se
new_data$lwr <- new_data$fit - t_crit * new_data$se

gg_goat + geom_line(data = new_data, aes(x = Stw, y = fit, colour = Treatment)) + 
  geom_ribbon(data = new_data, 
              aes(x = Stw, y = fit, ymin = lwr, ymax = upr, 
                  fill = Treatment), alpha = 0.3, size = 0) + 
  scale_fill_discrete('Способ обработки',
                  breaks = c('intensive', 'standard'),
                  labels = c('Интенсивный', 'Стандартный'))

Уравнение приобретает такой вид

\[ Wt = 14.97 - 0.35Stw - 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"

Mod_goat_reduced_2 <- lm(Wt ~ Stw + Treatment, data = goat)

Общее уравнение модели

\[ Wt = 13.7017 -0.3514 Stw + 1.2649 I_{intensive} \] Постройте уравнения для каждого из уровней фактора Treatment….

library(car)
Anova(Mod_goat_reduced, type = 3)

## Anova Table (Type III tests)
## 
## Response: Wt
##             Sum Sq Df F value        Pr(>F)    
## (Intercept)  190.9  1   72.93 0.00000000028 ***
## Stw           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

SS Type I

• SS(A) для фактора A.

• SS(B | A) для фактора B, при условии, что A уже в модели.

• SS(AB | B, A) для взаимодействия AB.

В этом случае проверяется главный эффект фактора “A”, затем главный эффект фактора “B” после главного эффекта “A”, а затем эффект взаимодействия “AB” после эффектов факторов “A” и “B”.

Из-за последовательного характера и того, что два основных фактора проверяются в определенном порядке, этот тип сумм квадратов даст различные результаты в зависимости от того, какой фактор рассматривается первым.

SS Type II

• SS(A | B) для фактора A.

• SS(B | A) для фактора B.

Этот тип проверяет каждый главный эффект после другого главного эффекта.

Обратите внимание, что предполагается отсутствие значимого взаимодействия!

Сначала следует проверить значимость взаимодействия (SS(AB | A, B)) и только если AB не является значимым, продолжить анализ главных эффектов.

Если взаимодействие действительно отсутствует, то тип II статистически более эффективен, чем тип III

SS Type III

• SS(A | B, AB) для фактора A.

• SS(B | A, AB) для фактора B

Этот тип сумм квадратов проверяет наличие главного эффекта при условии, что в модель уже включен второй главный эффект и взаимодействия. Поэтому такой тип суммы квадратов надо применять при наличии значимых взаимодействий.

Однако зачастую интерпретация главного эффекта при наличии взаимодействий не представляет интереса (вообще говоря, при наличии значимого взаимодействия эффекты факторов не должны подвергаться дальнейшему анализу).

Если взаимодействия не значимы, то тип II дает более мощный тест.

Поскольку знчимого взаимодействия не выявлено, то можно применить SS Type II.

library(car)
Anova(Mod_goat_reduced, type = 2)

## Anova Table (Type II tests)
## 
## Response: Wt
##           Sum Sq Df F value   Pr(>F)    
## Stw         58.6  1   22.40 0.000032 ***
## 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).

Мы уже рассмотрели довольно много статистических методов:

• Простая регрессия;
  
• Множественная регрессия;
  
• Регрессионный анализ с дискретными и непрерывными предикторами.

Все они являются реализациями одного и того же метода — общие линейные модели.

NB! “Общие линейные модели” (General linear models) нельзя путать с “обобщенными линейными моделями” (Generalized linear models).

Все разновидности линейных моделей описываются общей формулой

\[ \textbf{y} = \textbf{X} \textbf{b} + \textbf{e} \]

• \(\textbf{y}\) — вектор наблюдаемых значений
• \(\textbf{X}\) — модельная матрица, характеризующая “поведение” предикторов
• \(\textbf{b}\) — вектор коэффициентов модели
• \(\textbf{e}\) — вектор, описывающая случайные отклонения

В случае общих линейных моделей подбор параметров модели (вектор \(\textbf{b}\)) производится методом наименьших квадратов.

\[\mathbf{y} = \mathbf{X} \mathbf{b} + \mathbf{e}\]

\[ \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix} = \begin{pmatrix} 1 & x_{11} & \cdots & x_{p-1\ 1} \\ 1 & x_{12} & \cdots & x_{p-1\ 2} \\ \vdots & & & \\ 1 & x_{1n} & \cdots & x_{p-1\ n} \end{pmatrix} \times \begin{pmatrix} b _0 \\ b _1 \\ b _2 \\ \vdots \\ b_{p-1} \end{pmatrix} + \begin{pmatrix} e _1 \\ e _2 \\ \vdots \\ e _n \end{pmatrix} \]

• \(\mathbf{y}\) — вектор значений зависимой переменной, один столбец в n строк (n - число наблюдений)
• \(\mathbf{X}\) — модельная матрица с n строк. Ее первый столбец содержит единицы, далее по столбцу для каждой из переменных в модели. Количество столбцов равно количеству параметров модели.
• \(\mathbf{b}\) — вектор из \(p\) коэффициентов линейной модели
• \(\mathbf{e}\) — вектор остатков, один столбец в n строк

Разнообразие методов во многом определяется устройством модельной матрицы \(\textbf{X}\), то есть тем, какова природа предикторов.

• Непрерывные, дискретные и их взаимодействия - Общие линейные модели
  
• Только непрерывные - Регрессионный анализ
  
• Только дискретные - Дисперсионный анализ (ANOVA)
  
• Дискретные и непрерывные, без взаимодействия - Анализ ковариации (ANCOVA)

Что мы уже о них знаем:

• Принципиальных отличий этой формы линейных моделей от обычной множественной регрессии нет.
  
• Дискретные предикторы входят в модель, как переменные-индикаторы (dummy-variable).
  
• Между дискретными и непрерывными предикторами может существовать взаимодействие.
  
• Наличие взаимодействия необходимо специально оценивать.
• Визуализация таких моделей включает несколько линий регрессии, соответствующих разным градациям дискретного фактора.

Иногда влияние дискретного фактора зависит от дополнительной непрерывной переменной — ковариаты.

Например, испытывая влияние двух типов удобрений на суммарный урожай на делянках разной площади, мы должны учесть площадь делянки, которую включим в анализ как ковариату.

Площадь делянки — это ковариата.

Важное ограничение: в ANCOVA ковариата не должна взаимодействовать с дискретным предиктором.

Различается ли объем легких разных возрастных групп пациентов, которых готовят к операции?

Загрузите с сайта данные tlc.csv и поместите файл в папку data в рабочей директории.

{Данные: Altman, 1991 \
Источник: пакет }

tlc <- read.table('data/tlc.csv', sep = ';', header = TRUE)
str(tlc)

## 'data.frame':    32 obs. of  4 variables:
##  $ age   : int  35 11 12 16 32 16 14 16 35 33 ...
##  $ sex   : int  1 1 2 1 1 1 1 2 1 1 ...
##  $ height: int  149 138 148 156 152 157 165 152 177 158 ...
##  $ tlc   : num  3.4 3.41 3.8 3.9 4 4.1 4.46 4.55 4.83 5.1 ...

Переменные:

• age – возраст.
  
• sex – пол (1 - женский; 2 - мужской)
• height – рост (см)
• tlc– объем легких (л)

# Создаем переменную, кодирующую возрастную группу
tlc$age_group[tlc$age < 20] <- 'teenagers'
tlc$age_group[tlc$age < 30 & tlc$age >= 20] <- 'young'
tlc$age_group[tlc$age >= 30] <- 'adult'

# Создаем фактор с удобным порядком градаций
tlc$age_group <- factor(tlc$age_group, levels = c('teenagers', 'young', 'adult'))

table(tlc$age_group, tlc$sex)

##            
##             1 2
##   teenagers 4 2
##   young     5 6
##   adult     7 8

Нас интересует, зависит ли объем легких (tlc) от возрастной группы (age_group).

То есть модель должна быть основана только на дискретном предикторе с тремя градациями.

Можно ли построить модель такого вида?

\[ tlc_i = b_0 + b_1 \cdot I_{young\,i} + b_3 \cdot I_{adult\,i} + e_i \]

то есть модель должна быть основана только на дискретном предикторе с тремя градациями.

Моделями с дискретными предикторами занимается дисперсионный анализ. С ним мы детально познакомимся позднее.

НО! Мы уже знаем, что принципиальной разницы между обычным регрессионным анализом и регрессионным анализом с переменными-индикаторами нет.

Mod <- lm(tlc ~ age_group, data = tlc)
summary(Mod)

## 
## Call:
## lm(formula = tlc ~ age_group, data = tlc)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6920 -0.6505 -0.0367  0.8625  2.2500 
## 
## Coefficients:
##                Estimate Std. Error t value     Pr(>|t|)    
## (Intercept)       4.037      0.496    8.13 0.0000000057 ***
## age_groupyoung    3.163      0.617    5.13 0.0000178258 ***
## age_groupadult    2.055      0.587    3.50       0.0015 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.22 on 29 degrees of freedom
## Multiple R-squared:  0.475,  Adjusted R-squared:  0.439 
## F-statistic: 13.1 on 2 and 29 DF,  p-value: 0.0000865

Мы помним, что модель \(y_i = b_0 + b_1x_i + e_i\) можно записать в матричном виде

\[ \mathbf{y} = \mathbf{X}{\mathbf b} + \mathbf{e} \] или

\[ \widehat{\mathbf{y}} = \mathbf{X}{\mathbf b} \]

Как устроена модельная матрица \(\mathbf{X}\), если в модели только дискретные предикторы?

Вектор коэффициентов

coef(Mod)

##    (Intercept) age_groupyoung age_groupadult 
##           4.04           3.16           2.06

Для объекта, относящегося к базовому уровню

X[2,]

##    (Intercept) age_groupyoung age_groupadult 
##              1              0              0

\[ \widehat{tlc} = 1 \cdot 4.04 + 0 \cdot 3.16 + 0 \cdot 2.06 = 4.04 \]

\[ \widehat{\mathbf{y}} = \mathbf{X}{\mathbf b} \]

Вектор коэффициентов

coef(Mod)

##    (Intercept) age_groupyoung age_groupadult 
##           4.04           3.16           2.06

Для объекта, относящегося к young

X[12,]

##    (Intercept) age_groupyoung age_groupadult 
##              1              1              0

\[ \widehat{tlc} = 1 \cdot 4.04 + 1 \cdot 3.16 + 0 \cdot 2.06 = 4.04 + 3.16=7.20 \]

\[ \widehat{\mathbf{y}} = \mathbf{X}{\mathbf b} \]

Вектор коэффициентов

coef(Mod)

##    (Intercept) age_groupyoung age_groupadult 
##           4.04           3.16           2.06

Для объекта, относящегося к adult

X[1, ] 

##    (Intercept) age_groupyoung age_groupadult 
##              1              0              1

\[ \widehat{tlc} = 1 \cdot 4.04 + 0 \cdot 3.16 + 1 \cdot 2.06 = 4.04 + 2.06=6.10 \]

\[ \widehat{\mathbf{y}} = \mathbf{X}{\mathbf b} \]

Это будет график средних значений зависимой переменной для каждой градации дискретного фактора.

# Создаем искусственный датасет со всеми возможными значениями предиктора
new_data <- data.frame(age_group = factor(levels(tlc$age_group), 
                                          levels = levels(tlc$age_group)))
# Предсказанные моделью средние значения зависимой переменной
new_data$fit <- predict(Mod, newdata = new_data, se.fit = TRUE)$fit

# Стандартные ошибки
new_data$se <- predict(Mod, newdata = new_data, se.fit = TRUE)$se.fit

# Критические значения t для расчета доверительных интервалов 
t_crit <- qt(0.975, df = nrow(tlc) - length(coef(Mod)))

# Границы доверительных интервалов
new_data$upr <- new_data$fit + t_crit * new_data$se
new_data$lwr <- new_data$fit - t_crit * new_data$se

library(ggplot2)
ggplot(new_data, aes(x = age_group, y = fit)) +
  geom_bar(stat = 'identity', aes(fill = age_group)) +
  geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.2) +
  scale_x_discrete(labels = c('Подростки', 'Молодые', 'Взрослые')) +
  guides(fill = 'none') +
  labs(x = 'Возрастная группа', y = 'Объем легких')

Проведите диагностику этой модели.

Mod_diag <- fortify(Mod) # Создаем датафрейм с диагностическими данными
Mod_diag$height <- tlc$height # Переменная, не вошедшая в модель

ggplot(Mod_diag, aes(x = height, y = .stdresid)) +
  geom_point() + geom_hline(yintercept = 0) +
  geom_smooth(method = 'lm')

## `geom_smooth()` using formula = 'y ~ x'

Виден очевидный паттерн в остатках!

Изученные пациенты отличались не только возрастом, но и размером тела (height).
Этот параметр может тоже влиять на объем легких.

Необходимо провести ANCOVA.

Mod_cov <- lm(tlc ~ age_group * height, data = tlc)

В этой модели мы заложили наличие взаимодействия между дискретным и непрерывным предиктором.

ANCOVA подразумевает отсутствие взаимодействия между дискретным предиктором, который находится в фокусе нашего исследования, и ковариатой.

Нужен тест на наличие взаимодействия.

drop1(Mod_cov, test = 'F')

## Single term deletions
## 
## Model:
## tlc ~ age_group * height
##                  Df Sum of Sq  RSS  AIC F value Pr(>F)
## <none>                        25.8 5.10               
## age_group:height  2      2.59 28.4 4.16     1.3   0.29

Значимого взаимодействия нет! Можно делать ANCOVA.

Модель в ANCOVA ничем не отличается от любой другой регрессионной модели, включающей непрерывные и дискретные предикторы без их взаимодействия.

Mod_cov_2 <- lm(tlc ~ age_group + height, data = tlc)

X <- model.matrix(Mod_cov_2)
head(X)

##   (Intercept) age_groupyoung age_groupadult height
## 1           1              0              1    149
## 2           1              0              0    138
## 3           1              0              0    148
## 4           1              0              0    156
## 5           1              0              1    152
## 6           1              0              0    157

Здесь помимо колонок с переменными-индикаторами (age_groupyoung и age_groupadult) еще присутствует колонка со значениями непрерывного предиктора (height).

Аналогичная модельная матрица будет в любых других регрессионных моделях, включающих непрерывные и дискретные предикторы без их взаимодействия.

# Датафрейм с диагностическими данными
Mod_cov_2_diag <- fortify(Mod_cov_2)

ggplot(Mod_cov_2_diag, aes(x = .fitted, y = .stdresid)) +
  geom_point() + geom_hline(yintercept = 0) +
  geom_smooth(method = 'lm')

## `geom_smooth()` using formula = 'y ~ x'

Есть намек на гетероскедастичность!

ggplot(Mod_cov_2_diag, aes(x = height, y = .stdresid)) +
  geom_point() +
  geom_hline(yintercept = 0) +
  geom_smooth(method = 'lm')

## `geom_smooth()` using formula = 'y ~ x'

ggplot(Mod_cov_2_diag, aes(x = age_group, y = .stdresid)) +
  geom_boxplot() +
  geom_hline(yintercept = 0)

summary(Mod_cov_2)

## 
## Call:
## lm(formula = tlc ~ age_group + height, data = tlc)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9247 -0.5238 -0.0729  0.5184  2.1978 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -6.9258     2.9293   -2.36  0.02523 *  
## age_groupyoung   1.8991     0.6107    3.11  0.00427 ** 
## age_groupadult   0.7197     0.6011    1.20  0.24124    
## height           0.0718     0.0190    3.78  0.00076 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.01 on 28 degrees of freedom
## Multiple R-squared:  0.653,  Adjusted R-squared:  0.615 
## F-statistic: 17.5 on 3 and 28 DF,  p-value: 0.00000132

Фрагмент вывода функции summary()

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -6.9258     2.9293   -2.36  0.02523 *  
age_groupyoung   1.8991     0.6107    3.11  0.00427 ** 
age_groupadult   0.7197     0.6011    1.20  0.24124    
height           0.0718     0.0190    3.78  0.00076 ***

За базовый уровень взято значение age_group = teenagers.

Коэффициенты модели показывают, что объем легких при попарном сравнении с базовым уровнем (teenagers) статистически значимо выше только у young.

Введение в модель ковариаты помогло избежать ложных выводов!

Можем воспользоваться частными F-тестами, чтобы протестировать значимость предикторов в целом.

Anova(Mod_cov_2)

## Anova Table (Type II tests)
## 
## Response: tlc
##           Sum Sq Df F value  Pr(>F)    
## age_group   13.8  2     6.8 0.00392 ** 
## height      14.5  1    14.3 0.00076 ***
## Residuals   28.4 28                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Модель включает помимо дискретного предиктора еще и непрерывную ковариату.

Предсказания такой модели невозможно представить в виде простой столбчатой диаграммы, отражающей обычные средние значения зависимой переменной для каждой градации фактора.

Необходимо вместо обычных средних использовать скорректированные средние (Adjusted means).

• Нет разницы между моделями с дискретными и непрерывными предикторами
• При анализе влияния дискретного фактора часто необходимо учитывать влияние ковариаты.
• Пр тестировании гипотезы надо помнить о существовании двух типов тестирования: последовательное vs иерархическое.

• Must read paper! , 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