Пуассоновская регрессия
Максим Кузнецов
18.07.2019
Данные
Предположим, мы хотим оценить связь между количеством дней, проведенных в стационаре, и средней величиной систолического артериального давления у пациентов с артериальной гипертензией. Чтобы сделать это, мы проводим наблюдательное исследование, в котором в течение года ежемесячно измеряем артериальное давление (с последующим расчетом среднего значения за год) и регистрируем суммарное количество дней, проведенное каждым пациентом в стационаре за этот год.
Всего за время исследования мы пронаблюдали 100 пациентов с АГ (n = 100).
NB! Данные были сгенерированы и не имеют отношения к реальным пациентам (и к реальности как таковой).
set.seed(1000)
y <- rpois(100, lambda = 6)
x <- log(y) * 15 + 120 + rnorm(100, mean = 0, sd = 10)
df <- data.frame(x = x, y = y)Посмотрим, как выглядят полученные нами данные (первые 5 наблюдений из 100).
library(kableExtra)
df1 <- head(df, 5)
names(df1) <- c('Среднее САД', 'N дней в стационаре')
df1 %>%
kable(row.names = T) %>% kable_styling(full_width = F)| Среднее САД | N дней в стационаре | |
|---|---|---|
| 1 | 146.6587 | 5 |
| 2 | 140.2447 | 8 |
| 3 | 140.4556 | 3 |
| 4 | 139.2256 | 7 |
| 5 | 147.8822 | 6 |
Теперь посмотрим на распределения, полученных нами переменных.
par(mfrow = c(1,2))
y = 'Плотность'
x1 = 'N дней в стационаре'
x2 = 'Среднее САД'
hist(df$x, col = "steelblue2", main = NULL, xlab = x2, ylab = y, freq = F)
lines(density(df$x), col = 'red')
hist(df$y, col = "wheat2", main = NULL, xlab = x1, ylab = y, freq = F)
lines(density(df$y), col = 'red')
Посторим диаграмму рассеяния, чтобы визуально оценить связь между переменными.
plot(df$x, df$y, col = rgb(0.158, 0.168, 0.184, alpha = 0.3), pch = 20, cex = 2,
ylab = 'N дней в стационаре', xlab = 'Среднее САД')
На гистограммах выше мы видим, что значения среднего САД у наших пациентов имеют “нормальное” распределение; это можно дополнительно проверять, но здесь я не буду этого делать, поскольку я знаю точно, что 1) данные сгенерированы из нормального распределения и 2) в данном случае нам не особо важно распределение независимой переменной. Однако, нам очень важно распределение зависимой перменной – количества дней в стационаре, оно явно асимметричное (кстати, данная переменная совершенно точно имеет распределение отличное от нормального, так как эта выборка сгенерирована таким образом, на всякий случай, в тесте Шапиро-Уилка p = 0.0215687).
Что это означает? Мы не можем использовать линейную регрессию для моделирования связи (= ассоциации, прямого отношения к причино-следственным связям эти термины не имеют) между длительностью нахождения в стационаре и значением САД (это не совсем так, но мы опустим данный момент).
Что делать? В таком случае у нас есть достаточно много вариантов, но мы сосредоточимся на одном – пуассоновской регрессии.
Вот несколько интересных статей журнала “Квант”, посвященных распределению Пуассона: 1, 2.
Немного теории: пуассоновская регрессия является представителем большой (и крайне важной) группы регрессионных моделей, которая называется обобщенные линейные модели (GLM), в которую также входит логистическая регрессия. Все представители этого семейства имеют следующую структуру:
Случайный компонент – модель для зависимой переменной (в нашем случае для моделирования количества дней в стационаре мы будем использовать пуассоновское распределение, в общем случае – распределение зависимой переменной может быть смоделировано другим представителем экспоненциального семейства).
Систематический (неслучайный) компонент – всегда линейная комбинация предикторов, независимо от их распределения (в нашем случае он только один – САД). Для того, чтобы представить это, попробуйте вспомнить, что находится с правой стороны в простом линейном регрессионном уравнении.
Функция связи – функция, обеспечивающая линейную связь между 1-ым и 2-ым компонентом (в нашем случае это будет натуральный логарифм).
Такми образом, в общем виде модель можно представить следующим образом:
\[g(Y_i) = \sum^p_{k=1}{X_{i,k}\beta_k}\]
где \(g(\cdot)\) – функция связи.
Для нашего случая, регрессионное уравнение будет иметь следующий вид:
\[log(Y_i) = \beta_1 \times X_i + \beta_0\]
где \(\beta_o\), \(\beta_1\) – коэффициенты модели.
Зададим такой вопрос: а что в этой модели нас интересует? На самом деле, исходя из главной задачи нашего исследования, нас интересует значение коэффициента \(\beta_1\), который и будет описывать ассоциацию между САД и количеством дней в стационаре.
Теперь внимание (самый главный момент в использовании так называемых параметрических методов): чтобы использовать данную модель для оценки связи интересующих нас переменных (т.е. для того чтобы оценить с помощью данной модели интересующий нас параметр \(\beta_1\)) мы должны постулировать следующие предположения:
Предположение о независмости наблюдений: все наши переменные для каждого наблюдения должны быть независимы друг от друга и иметь одинаковое распределение (т.н. iid sampling assumption). Это предположение является валидным для случайной выборки.
Наша зависимая переменная должна иметь распределение Пуассона. Что это означает? Распределение Пуассона имеет один единственный параметр \(\lambda\) – rate parameter, причем, 1) \(\lambda \gt 0\); 2) матожидание (среднее значение в генеральной совокупности) = \(\lambda\) И дисперсия = \(\lambda\).
Как минимум, для того чтобы наша модель была валидна (и следовательно коэффициент \(\beta_1\), мы могли интерпретировать) нужно, чтобы мы не сомневались во 2-ом пункте. Попробуем проверить этот постулат на основе полученных нами данных (дадим оценки матожиданию и дисперсии):
Выборочное среднее числа дней в стационаре равно 5.79.
Выборочная (несмещенная) дисперсия равна 5.2180808.
Мы не видим больших отличий между выборочной дисперсией и средней и, наверное, можем использовать нашу модель.
Давайте построим ее!
fit <- glm(y ~ x, data = df, family = 'poisson')
summary(fit)$coef %>%
as.table() -> coeffs
colnames(coeffs) <- c('Оценка', 'Стандартная\nошибка', 'z-статистика', 'p-значение')
rownames(coeffs) <- c('$\\beta_0$', '$\\beta_1$')
coeffs %>%
kable(caption = 'Коэффициенты регрессионного уравнения', escape = F) %>% kable_styling(full_width = F)| Оценка | Стандартная ошибка | z-статистика | p-значение | |
|---|---|---|---|---|
| \(\beta_0\) | -1.0164527 | 0.5274047 | -1.927273 | 0.0539456 |
| \(\beta_1\) | 0.0187858 | 0.0035313 | 5.319753 | 0.0000001 |
NB! При использовании обобщенных линейных моделей (на самом деле, как и в случае линейной регрессии) все самое важное и интересное начинается с момента, когда мы построили модель.
Итак, мы построили модель, которая (возможно) описывает связь между САД и количеством дней, проведенных пациентом в стационаре. Теперь зададимся вопросом, а насколько валидна данная модель (чтобы иметь право интерпретации оценки коэффициента \(\beta_1\) и его статистической значимости).
- Нам стоит проверить, как модель в целом объясняет имеющиеся у нас данные. Для этого в случае пуассоновской регрессии можно использовать специальный график (rootogram), который позволяет сравнить эмпирические частоты значений зависимой переменной для наших наблюдений с частотами на основе предсказанных моделью значений зависимой переменной.
library(countreg)
library(latex2exp)
rootogram(fit, main = NA, ylab = TeX('\\sqrt{Частота}'),
max = 10, xlab = 'N дней в стационаре')
Нам стоит вспомнить о сделанных нами ранее предположениях и еще раз подумать, могут ли они нарушаться. Дело в том, что полученная нами оценка коэффициента \(\beta_1\) была получен потому что мы постулировали: независимость и одинаковое распределение наших наблюдений и параметрические предположения (зависимая перменная имеет пуассоновское распределение и линейность ассоциации с учетом функции связи).
Чтобы сделать статистический вывод о значимости коэффициента \(\beta_1\) нам нужно: во-первых снова обратиться к нашим предположениям (если коэффициент невалиден, то о каком выводе может идти речь) и обратить внимание на 3-ий столбец таблицы с регрессионными коэффициентами. Мы видим, что там находится z-статистика для коэффициентов, а это означает, что при осуществлении статистического вывода мы полагаемся на еще одно предположение о том, что оценка коэффициента \(\beta_1\) (она является случайной величиной, подумайте почему) имеет нормальное распределение, но почему мы можем сделать такое предположение? Дело в том, что в данном случае мы при статистическом выводе опираемся на центральную предельную теорему (т.е. оценка коэффициента \(\beta_1\) имеет асимптотически нормальное распределение), что в свою очередь означает, что мы должны иметь достаточно большое количество наблюдений (самое скверное, что мы не можем сказать, насколько оно должно быть большое).
Чтобы в целом оценить качество модели, которое, в том числе, связано с наличием/отсутствием нарушений наших предположений мы можем применить множество стратегий, но данная тема требует более глубокого погружения в статистику.
Снова посторим диаграмму рассеяния и добаваим линию регрессии (обратите внимание, что она не является прямой, как в случае линейной регрессии)
y_hat <- predict(fit, type = 'response', newdata = data.frame(x = seq(120, 200, by = 0.1)))
plot(df$x, df$y, col = rgb(0.158, 0.168, 0.184, alpha = 0.3), pch = 20, cex = 2,
ylab = 'N дней в стационаре', xlab = 'Среднее САД')
lines(seq(120, 200, by = 0.1), y_hat, col = 'red')