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

• Как все эти изменения закодировать в уравнении линейной модели?

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

\(Y = b_0 + b_1X\)

\(Y = b_0 + b_1X + b_2 Group2\)

\(Group2 = 1\) если наблюдение из группы 2, и \(Group2 = 0\) в остальных случаях

Уравнение для группы 1: \(Y = b_0 + b_1X\)

Уравнение для группы 2: \(Y = (b_0 + b_2) + b_1X\)

\(Y = b_0 + b_1X + b_2X Group2\)

\(Group2 = 1\) если наблюдение из группы 2, и \(Group2 = 0\) в остальных случаях

Уравнение для группы 1: \(Y = b_0 + b_1X\)

Уравнение для группы 2: \(Y = b_0 + (b_1 + b_2)X\)

\(Y = b_0 + b_1X + b_2 Group2 + b_3X Group2\)

\(Group2 = 1\) если наблюдение из группы 2, и \(Group2 = 0\) в остальных случаях

Уравнение для группы 1: \(Y = b_0 + b_1X\)

Уравнение для группы 2: \(Y = (b_0 + b_2) + (b_1 + b_3)X\)

Не называйте взаимодействие факторов “корреляцией” или “связью”!

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

Т.е. в данном случае можно сказать: “Для группы 1 и группы 2 зависимость Y от Х выглядит по-разному”

Пуромицин - антибиотик пуринового ряда, ингибитор синтеза белков. Эти данные — о том, как меняется активность фермента галактозил трансферазы под воздействием пуромицина (Treloar 1974). Измеряли скорость реакции в зависимости от концентрации субстрата на мембранах аппарата Гольджи из богатых мембранами фракций из печени крыс.

• conc — концентрация пуромицина
• rate — скорость химической реакции
• state — индикатор того, обработаны ли клетки пуромицином

library(readxl)
Puromycin <- read_excel("data/Puromycin.xlsx")
head(Puromycin)

## # A tibble: 6 x 3
##    conc  rate state  
##   <dbl> <dbl> <chr>  
## 1  0.02    76 treated
## 2  0.02    47 treated
## 3  0.06    97 treated
## 4  0.06   107 treated
## 5  0.11   123 treated
## 6  0.11   139 treated

sapply(Puromycin, class)

##        conc        rate       state 
##   "numeric"   "numeric" "character"

colSums(is.na(Puromycin))

##  conc  rate state 
##     0     0     0

nrow(Puromycin)

## [1] 23

table(Puromycin$state, Puromycin$conc)

##            
##             0.02 0.06 0.11 0.22 0.56 1.1
##   treated      2    2    2    2    2   2
##   untreated    2    2    2    2    2   1

library(ggplot2)
ggplot(Puromycin, aes(x = state, y = conc)) + geom_boxplot()

Постройте график зависимости скорости химической реакции от концентрации

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

library(ggplot2)
ggplot(Puromycin, aes(x = conc, y = rate, colour = state)) + 
  geom_point()
ggplot(Puromycin, aes(x = log(conc), y = rate, colour = state)) + 
  geom_point()

Response ~ Continuous + Categorical + Continuous : Categorical (полная запись)

Response ~ Continuous * Categorical(сокращенная запись)

Continuous : Categorical — Кодовое обозначение взаимодействия. Добавление взаимодействия непрерывного и дискретного факторов в модель означает, что возможен разный угол наклона прямых для разных групп.

Puromycin$lc <- log(Puromycin$conc)
M1 <- lm(rate ~ lc + state + lc:state, data = Puromycin)
# M1 <- lm(rate ~ lc * state, data = Puromycin) # То же самое
summary(M1)

## 
## Call:
## lm(formula = rate ~ lc + state + lc:state, data = Puromycin)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -17.02  -6.38  -1.00   7.62  13.32 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         209.19       4.45   46.97  < 2e-16 ***
## lc                   37.11       1.97   18.86  9.2e-14 ***
## stateuntreated      -44.61       6.81   -6.55  2.9e-06 ***
## lc:stateuntreated   -10.13       2.94   -3.45   0.0027 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.15 on 19 degrees of freedom
## Multiple R-squared:  0.968,  Adjusted R-squared:  0.963 
## F-statistic:  191 on 3 and 19 DF,  p-value: 2.27e-14

• Взаимодействие достоверно!

# Данные для графиков остатков
M1_diag <- fortify(M1)

1) График расстояния Кука

ggplot(M1_diag, aes(x = 1:nrow(M1_diag), y = .cooksd)) +
  geom_bar(stat = "identity")

2) График остатков от предсказанных значений

gg_resid <- ggplot(data = M1_diag, aes(x = .fitted, y = .stdresid)) +
  geom_point() + geom_hline(yintercept = 0)
grid.arrange(gg_resid, gg_resid + geom_smooth(),
             gg_resid + geom_smooth(method = "lm"), nrow = 1)

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

3) Графики остатков от предикторов в модели и не в модели

library(gridExtra)
grid.arrange(gg_resid + aes(x = lc), 
             ggplot(M1_diag, aes(x = state, y = .stdresid)) + geom_boxplot(), 
             nrow = 1)

4) Квантильный график остатков

library(car)
qqPlot(M1)

## [1]  2 10

Что не так с этим графиком? Это наша модель?

ggplot(Puromycin, aes(x = lc, y = rate, colour = state)) + 
  geom_point() + 
  geom_smooth(method = "lm")

• Это не наша модель — здесь разный угол наклона, а у нашей модели график должен быть другой, с одинаковым углом наклона. (По умолчанию ggplot рисует отдельные регрессии для каждого сабсета)

Постройте график предсказаний модели по этим данным

Данные для графика

library(plyr)
NewData <- ddply(
  Puromycin, .variables = .(state), summarise, 
  lc = seq(min(lc), max(lc), length = 100))

# предсказанные значения
Predictions <- predict(M1, newdata = NewData, se.fit = TRUE)
NewData$fit <- Predictions$fit

# стандартные ошибки
NewData$SE <- Predictions$se.fit

# доверительный интервал
NewData$upr <- NewData$fit + 1.96 * NewData$SE
NewData$lwr <- NewData$fit - 1.96 * NewData$SE

# Обратная трансформация предиктора
NewData$conc <- exp(NewData$lc)

ggplot(NewData, aes(x = conc, y = fit)) + 
  geom_ribbon(alpha = 0.2, aes(ymin = lwr, ymax = upr, group = state)) +
  geom_line(aes(colour = state)) +
  geom_point(data = Puromycin, aes(x = conc, y = rate, colour = state))

Запишите уравнение этой модели и уравнения для каждой группы

coef(M1)

##       (Intercept)                lc    stateuntreated 
##             209.2              37.1             -44.6 
## lc:stateuntreated 
##             -10.1

\(rate = 209.2 + 37.1 * lc - 44.6 * stateuntreated - 10.1 * lc * stateuntreated\)

\(stateuntreated = 1\) если \(state = untreated\) и ноль в остальных случаях

Уравнение для treated: \(rate = 209.2 + 37.1 * lc\)

Уравнение для untreated: \(rate = (209.2- 44.6) + (37.1 - 10.1) * lc\),
то есть \(rate = 164.6 + 27 * lc\)

Это делается как обычно при выборе моделей: при помощи anova() или drop1()

M2 <- lm(rate ~ lc + state, data = Puromycin)

# anova(M1, M2, test="F")
drop1(M1, test = "F")

## Single term deletions
## 
## Model:
## rate ~ lc + state + lc:state
##          Df Sum of Sq  RSS AIC F value Pr(>F)   
## <none>                1591 105                  
## lc:state  1       996 2587 115    11.9 0.0027 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Если убрать взаимодействие, то модель становится значительно хуже. Поэтому взаимодействие сохраняем

Как вес новорожденных зависит от возраста матери и того, курит ли она

• age — возраст матери
• lwt — вес матери до беременности
• race — раса (1-белые, 2-черные, 3-другие)
• smoke — курение во время беременности (1-да,2-нет)
• ptl — число предыдущих преждевременных родов
• ht — гипертензия
• ui — гипертонус матки
• ftv — число визитов к врачу в последний триместр
• bwt — вес новорожденного, г

wt <- read.table("data/birthwt.csv", header = TRUE, sep = ";")

• Исследуйте данные о весе новорожденных
• Постройте модель зависимости веса новорожденных от возраста матери и взаимодействия
• Проверьте условия применимости этой модели
• Упростите модель, если это возможно
• Напишите общее уравнение и отдельные уравнения модели
• Постройте график предсказаний модели

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: int  0 0 1 1 1 0 0 0 1 1 ...
##  $ 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 33 34 35
##   0  2  3  5  7  3  8 10  7  8  9 10 12  3  2  5  4  5  2  4  2  0  1
##   1  1  0  2  5  7  8  8  5  5  4  3  3  5  1  4  3  2  3  2  1  1  1
##    
##     36 45
##   0  2  1
##   1  0  0

table(wt$smoke, wt$race)

##    
##      1  2  3
##   0 44 16 55
##   1 52 10 12

table(wt$smoke, wt$ui)

##    
##       0   1
##   0 100  15
##   1  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 = lwt)
gg_dot + aes(x = bwt)

Проверка на колинеарность

wt_mod_1 <- lm(bwt ~ age + smoke, data = wt)
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")])

1) График расстояния Кука

ggplot(wt_mod_2_diag, aes(x = 1:nrow(wt_mod_2_diag), y = .cooksd)) +
  geom_bar(stat = "identity")

2) График остатков от предсказанных значений

gg_resid <- ggplot(data = wt_mod_2_diag, aes(x = .fitted, y = .stdresid)) +
  geom_point() + geom_hline(yintercept = 0)
gg_resid

3) Графики остатков от предикторов в модели и не в модели

library(gridExtra)
grid.arrange(gg_resid + aes(x = age), 
             gg_resid + aes(x = lwt), 
             nrow = 1)

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)

4) Квантильный график остатков

qqPlot(wt_mod_2)

## [1] 131 132

coef(wt_mod_2)

## (Intercept)         age      smoke1  age:smoke1 
##      2406.1        27.7       798.2       -46.6

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

bwt = 2406.06 + 27.73age + 798.17smoke1 - 46.57age:smoke1

• Для некурящих: \(N = 2406.06 + 27.73 * age\)
• Для курящих: \(N = 3204.23 - -18.84 * age\)

Постройте график предсказаний модели

ggplot(NewData, aes(x = age, y = fit)) +
  geom_ribbon(alpha = 0.2, aes(ymin = lwr, ymax = upr, group = smoke)) +
  geom_line(aes(colour = smoke))

ggplot(NewData, aes(x = age, y = fit)) +
  geom_ribbon(alpha = 0.2, aes(ymin = lwr, ymax = upr, group = smoke)) +
  geom_line(aes(colour = smoke)) +
  geom_point(data = wt, aes(x = age, y = bwt, colour = smoke))

• У дискретных и непрерывных факторов могут быть взаимодействия, поэтому анализ нужно начинать с более полной модели (чтобы можно было протестировать и исключить недостоверные взаимодействия)

• 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