• Линейные модели 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 x 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 Treatmentstandard 
##                 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

Привес коз при интенсивной обработке от глистов на 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