Статистические модели

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

Построение моделей является одной из центральных тем статистики как науки. Существует большое число видов статистических моделей, различающихся как по лежащим в их основе математическим принципам, так и по преимущественным областям применения. Однако несмотря на все свое разнообразие, статистические модели сходны в том, что они описывают взаимосвязь между случайными переменными. Как именно это происходит? Постараемся разобраться...

Пример простейшей статистической модели

Существует мало вещей и явлений, в отношении которых мы можем быть совершенно уверены - практически всегда присутствует некоторый уровень неопределенности, вызванный естественной изменчивостью, влиянием внешних факторов и погрешностью измерений. Предположим, что мы имеем дело с некоторой случайной количественной переменной \(Y\). Отдельные выборочные значения этой переменной будем обозначать как \(y_i\) . Индекс \(i\) изменяется от 1 до \(n\), где \(n\) - это объем выборки. В качестве примера рассмотрим систолическое кровяное давление у людей (выражается в мм ртутного столба). Очевидно, что уровень кровяного давления не является одинаковым у всех людей - при обследовании случайно сформированной выборки мы почти всегда будем наблюдать определенный разброс значений этой переменной, хотя некоторые значения будут встречаться чаше других. На рис. 1 приведен пример возможного распределения значений давления крови, измеренных у 100 случайно отобранных человек (различавшихся по возрасту, полу, массе тела и, возможно, каким-то другим характеристикам):

y <- c(
109.14, 117.55, 106.76, 115.26, 117.13, 125.39, 121.03,
114.03, 124.83, 113.92, 122.04, 109.41, 131.61, 103.93,
116.64, 117.06, 111.73, 120.41, 112.98, 101.20, 120.19,
128.53, 120.14, 108.70, 130.77, 110.16, 129.07, 123.46,
130.02, 130.31, 135.06, 129.17, 137.08, 107.62, 139.77,
121.47, 130.95, 138.15, 114.31, 134.58, 135.86, 138.49,
110.01, 127.80, 122.57, 136.99, 139.53, 127.34, 132.26,
120.85, 124.99, 133.36, 142.46, 123.58, 145.05, 127.83,
140.42, 149.64, 151.01, 135.69, 138.25, 127.24, 135.55,
142.76, 146.67, 146.33, 137.00, 145.00, 143.98, 143.81,
159.92, 160.97, 157.45, 145.68, 129.98, 137.45, 151.22,
136.10, 150.60, 148.79, 167.93, 160.85, 146.28, 145.97,
135.59, 156.62, 153.12, 165.96, 160.94, 168.87, 167.64,
154.64, 152.46, 149.03, 159.56, 149.31, 153.56, 170.87,
163.52, 150.97)
 
library(ggplot2)
ggplot(data = Data, aes(x = y)) + geom_histogram() + 
       ylab("Частота") + xlab("Давление, мм рт. ст.")

Рис. 1. Распределение 100 симулированных значений систолического кровяного давления (см. объяснения в тексте).

Откуда взялись данные, отображенные на рис. 1, мы узнаем чуть позже. Пока же примем, что это реальные значения систолического давления крови, измеренные у реальных людей. Распределение этих 100 значений по форме напоминает колокол, указывая на то, что, вероятно, мы имеем дело с нормально распределенной переменной (для пущей уверенности мы, конечно, могли бы применить и более формальные методы, чем простая визуальная оценка гистограммы, но пока допустим, что это распределение действительно является нормальным). Как известно, нормальное распределение полностью задается двумя параметрами: \(\mu\) - математическим ожиданием (=арифметическим средним), и \(\sigma\) - стандартным отклонением. То обстоятельство, что переменная \(Y\) распределена согласно нормальному закону, принято обозначать следующим образом: \(y_i \sim N(\mu, \sigma)\). Выражение \(y_i \sim N(\mu, \sigma)\) представляет собой пример простейшей статистической модели. Подобно другим параметрическим моделям, эта модель включает две составные части: систематическую и случайную. Систематическая часть представлена средним значением \(\mu\), которое отражает центральную тенденцию в данных. В свою очередь, стандартное отклонение характеризует степень вариации значений \(Y\), т.е. их разброс относительно среднего. Допуская, что наша выборка из 100 значений действительно происходит из нормально распределенной генеральной совокупности, мы можем использовать эту модель для предсказания давления крови: для всех людей предсказанное значение окажется одинаковым и будет равно \(\mu\). Однако чему именно равно это значение?

Греческие буквы \(\mu\) и \(\sigma\) обозначают среднее значение и стандартное отклонение генеральной совокупности, которые, как правило, нам неизвестны. Тем не менее, мы можем оценить значения параметров генеральной совокупности по соответствующим выборочным статистикам. Так, в случае с данными из рис. 1 среднее значение и стандартное отклонение составляют 135.16 мм рт. ст. и 16.96 мм рт. ст. соответственно, и мы можем записать нашу модель в виде \(y_i \sim N(135.16, 16.96)\). Эквивалентным способом записи этой модели является следующий:

\[y_i = 135.16 + \epsilon_i,\]

где \(\epsilon_i\) - это остатки модели, имеющие нормальное распределение со средним значением 0 и стандартным отклонением 16.96:  \(\epsilon_i \sim N(0, 16.96)\). Остатки рассчитываются как разница между реально наблюдаемыми значениями переменной \(Y\) и значениями, предсказанными моделью (в рассматриваемом примере \(\epsilon_i = y_i - 135.16\)). Из этого второго способа записи модели видно, что она представляет собой ни что иное, как линейную регрессионную модель, у которой нет ни одного предиктора (подробнее см. здесь). Модели, у которых нет предикторов, часто называют "нулевыми"  (англ. null models).

Симуляции на основе статистических моделей

В сущности, статистическая модель - это упрощенное математическое представление процесса, который, как мы полагаем, привел к генерации наблюдаемых значений изучаемой переменной. Это значит, что мы можем использовать модель для симуляции - процедуры, имитирующей действие моделируемого процесса и позволяющей тем самым искусственно сгенерировать новые значения изучаемой переменной, которые, как мы надеемся, будут обладать свойствами реальных данных. Причин, по которым мы хотели бы сгенерировать такие "искусственные данные", можно назвать много, но основные из них обычно сводятся к решению следующих практических задач:

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

С подробными примерами использования симуляций для решения этих и других задач (включая R-код) можно ознакомиться в великолепной книге Эндрю Гелмана (Andrew Gelman) и Дженнифер Хилл (Jennifer Hill) о регрессионном анализе с применением многоуровневых моделей (также известных как модели со смешанными эффектами и иерархические модели): Gelman A., Hill J. (2007) Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press. Здесь мы обратимся к последней из перечисленных выше задач и посмотрим, насколько хорошо значения кровяного давления, генерируемые нашей нулевой моделью, согласуются с данными из рис. 1.

В R новые данные на основе этой простой модели можно легко сгенерировать при помощи функции rnorm() (подробнее об этой функции см. здесь):

set.seed(101) # для воспроизводимости результата
y.new.1 <- rnorm(n = 100, mean = 135.16, sd = 16.96)

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

set.seed(101)
y.new.2 <- 135.16 + rnorm(n = 100, mean = 0, sd = 16.96)
 
# идентичны ли оба вектора?
all(y.new.1 == y.new.2)
[1] TRUE

Теперь необходимо вспомнить о том, что параметры нашей нулевой модели являются лишь точечными оценками истинных параметров, и что всегда будет присутствовать неопределенность в отношении того, насколько точны эти точечные выборочные оценки. В приведенных выше командах эта неопределенность не была учтена: при создании векторов y.new.1 и y.new.2 выборочные оценки среднего значения и стандартного отклонения кровяного давления рассматривались как истинные. В зависимости от стоящей задачи, такой подход может оказаться достаточным. Однако мы сделаем еще один шаг и учтем неопределенность в отношении точечных оценок параметров модели.

Для целей симуляции мы воспользуемся функцией lm(), которая предназначена для подгонки линейных регрессионных моделей. Ничего удивительного здесь нет - ведь мы уже знаем, что нашу простую модель кровяного давления можно рассматривать как линейную регрессионную модель, у которой нет ни одного предиктора:

y.lm <- lm(y ~ 1)
 
summary(y.lm)
Call:
lm(formula = y ~ 1) # "1" значит, что оценивается только свободный член
 
Residuals:
   Min     1Q Median     3Q    Max 
-33.96 -13.26   0.41  12.04  35.71 
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  135.157      1.696   79.69   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 16.96 on 99 degrees of freedom

Как следует из приведенных результатов, свободный член подогнанной модели (Intercept) в точности совпадает со средним значением данных из рис. 1 (135.16 мм рт. ст.), а стандартное отклонение остатков модели (Residual standard error) совпадает со стандартным отклонением этих данных (16.96 мм рт. ст.). Важно, однако, то, что теперь мы знаем также стандартную ошибку среднего среднего значения (1.696; см. столбец Std. Error в строке (Intercept)).

По определению, стандартная ошибка параметра - это стандартное отклонение [нормального] распределения значений этого параметра, рассчитанных по выборкам одинакового размера из той же генеральной совокупности. Мы можем использовать это обстоятельство для учета неопределенности в отношении точечных оценок параметров модели при симуляции новых данных. Так, зная выборочные оценки параметров и их стандартные ошибки, мы 1) можем сгенерировать несколько возможных значений этих параметров (т.е. составить несколько реализаций той же модели, варьируя значения параметров) и 2) симулировать новые данные на основе каждой из этих альтернативных реализаций модели. Функция sim() из пакета arm (является приложением к упомянутой выше книге Gelman and Hill (2007)) предназначена именно для этой задачи. Эта функция принимает модельные объекты классов lm, glm и др., извлекает оценки параметров и их соответствующие стандартные ошибки и возвращает заданное пользователем число альтернативных реализаций модельных параметров. Например, следующие команды позволят сгенерировать 5 альтернативных реализаций среднего кровяного давления и его стандартного отклонения для нашей простой модели:

install.packages("arm")
library(arm)
 
set.seed(102) # для воспроизводимости результата
y.sim <- sim(y.lm, 5)
 
# y.sim - объект класса S4, который содержит
# слоты coef (коэффициенты модели) и sigma 
# (станд. отклонения остатков модели):
str(y.sim)
Formal class 'sim' [package "arm"] with 2 slots
  ..@ coef : num [1:5, 1] 136 134 137 136 137
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : NULL
  ..@ sigma: num [1:5] 16.8 18.9 17.3 16.7 15
 
# Извлекаем альтернативные реализации коэффициентов модели из y.sim:
y.sim@coef
         [,1]
[1,] 136.4775
[2,] 134.3283
[3,] 136.7074
[4,] 136.0771
[5,] 137.3246
 
# Извлекаем альтернативные реализации станд. отклонения остатков:

y.sim@sigma
[1] 16.829, 18.870, 17.303, 16.743, 15.006

Конечно, 5 реализаций модели - это немного, чтобы сделать какие-либо убедительные выводы. Увеличим это число до 1000:

set.seed(102) # для воспроизводимости результата
y.sim <- sim(y.lm, 1000)
 
# Инициализация пустой матрицы, в которой мы будем сохранять
# данные, сгенерированные на основе 1000 альтернативных 
# реализаций модели:
y.rep <- array(NA, c(1000, 100))
 
# Заполняем матрицу y.rep симулированными данными:
for(s in 1:1000){
  y.rep[s, ] <- rnorm(100, y.sim@coef[s], y.sim@sigma[s])
}

Чтобы лучше понять, что мы только что сделали, изобразим гистограммы значений кровяного давления, симулированных на основе, например, первых 20 реализаций нулевой модели:

par(mfrow = c(5, 4), mar = c(2, 2, 1, 1))
for(s in 1:20){ hist(y.rep[s, ], xlab = "", ylab = "", breaks = 20, main = "")}

Рис. 2. Гистограммы значений кровяного давления, симулированных на основе 20 альтернативных реализаций  нулевой модели.

Рис. 2 дает некоторое представление о том, насколько значения кровяного давления, симулированные на основе альтернативных реализаций нашей простой модели сходны с реально наблюдаемыми значениями: некоторые гистограммы более сходны с гистограммой на рис. 1, некоторые - менее сходны. Однако больше всего нас интересует наличие систематических расхождений между свойствами этих симулированных значений и свойствами реально наблюдаемых данных. Если такие расхождения имеются, значит нулевая модель недостаточно хороша для описания вариации реально наблюдаемых значений.

Имеется несколько подходов для обнаружения этих расхождений. Один из них - расчет каких-либо описательных статистик для каждого из симулированных распределений и сравнение их с теми же статистиками у реально наблюдаемых данных. Так, например, мы можем рассчитать интерквартильный размах (ИКР) для каждого симулированного набора данных и сравнить полученное распределение из 1000 значений с ИКР реальных данных. Для расчета ИКР в R служит функция IQR():

# Расчет ИКР для каждого из 1000 симулированных
# распределений значений кровяного давления (подробнее о функции
# apply() см. здесь): 
test.IQR <- apply(y.rep, MARGIN = 1, FUN = IQR)
 
# Гистограмма 1000 значений ИКР (рис. 3):
hist(test.IQR, xlim = range(IQR(y), test.IQR),
     main = "ИКР", xlab = "", ylab = "Частота", breaks = 20)
lines(rep(IQR(y), 2), c(0, 100), col = "blue", lwd = 4)

Рис. 3. Гистограмма значений ИКР, рассчитанных для каждого из 1000 симулированных распределений кровяного давления. Вертикальная синяя линия показывает ИКР реально наблюдаемых значений кровяного давления (из рис. 1).

Из рис. 3 хорошо видно, что ИКР симулированных данных систематически занижены по сравнению с реальными данными. Это свидетельствует о том, что нулевая модель в целом недооценивает уровень вариации реальных значений кровяного давления. Причиной этому может быть то, что мы не учитываем влияние на кровяное давление каких-то важных факторов (например, возраст, пол, диета, состояние здоровья, и т.п.). Рассмотрим, как можно усложнить нашу нулевую модель, добавив в нее один из таких факторов.

Пример модели с одним количественным предиктором

Предположим, что помимо кровяного давления мы также измеряли у каждого испытуемого его/ее возраст (в годах):

# Значения возраста:
x <- rep(seq(16, 65, 1), each = 2)
 
# Объединяем значения возраста и давления крови в одну таблицу (для удобства):
Data <- data.frame(Age = x, BP = y)
 
ggplot(data = Data, aes(x = Age, BP)) + geom_point() + 
       geom_smooth(method = "lm", se = FALSE) +
       geom_rug(color = "gray70", sides = "tr") + 
       ylab("Частота") + xlab("Возраст, лет")

Рис. 4. Связь между возрастом и систолическим кровяным давлением. Небольшие отрезки светло-серого цвета по краям графика соответствуют наблюдаемым значениям возраста (вверху) и кровяного давления (справа). Подробнее об этом типе диаграмм см. здесь. Для визуализации тренда в данных добавлена регрессионная линия (синего цвета).

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

\[y_i \sim N(\beta_0 + \beta_1\times  Age_i, \sigma),\]

или, что то же самое:

\[y_i = \beta_0 + \beta_1\times  Age_i + \epsilon_i,\]

где \(\epsilon_i \sim N(0, \sigma). \)

Параметры \(\beta_0\) и \(\beta_1\) можно без труда оценить методом наименьших квадратов при помощи функции lm():

summary(lm(BP ~ Age, data = Data))
 
Call:
lm(formula = BP ~ Age, data = Data)
 
Residuals:
     Min       1Q   Median       3Q      Max 
-21.6626  -6.2509   0.0075   6.3125  17.3525 
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 94.85346    2.66731   35.56   <2e-16 ***
Age          0.99515    0.06204   16.04   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 8.953 on 98 degrees of freedom
Multiple R-squared:  0.7242, Adjusted R-squared:  0.7214 
F-statistic: 257.3 on 1 and 98 DF,  p-value: < 2.2e-16

Согласно полученным результатам, модель кровяного давления можно записать как

\[y_i \sim N(94.853 +0.995\times  Age_i, 8.953)\]

или

\[y_i = 94.853 +0.995\times  Age_i + \epsilon_i,\]

где \(\epsilon_i \sim N(0, 8.953).\)

Графически эта модель изображена выше на рис. 4 в виде линии тренда. Обратите внимание: помимо высокой значимости параметров подогнанной модели (P << 0.001 в обоих случаях), стандартное отклонение остатков составляет 8.853, что почти в 2 раза меньше, чем у нулевой модели (16.96, см. выше). Это указывает на то, что модель, включающая возраст как предиктор, гораздо лучше описывает вариацию значений кровяного давления у 100 обследованных испытуемых, чем наша исходная нулевая модель. Это заключение подтверждается также при выполнении симуляций вроде тех, что были описаны выше (код не привожу) - значение интерквартильного размаха данных из рис. 1 приходится на центр распределения симулированных значений ИКР, указывая на отсутствие систематических различий между симулированными и наблюдаемыми данными.

Рис. 5. Гистограмма значений ИКР каждого из 1000 распределений кровяного давления, симулированных на основе модели с возрастом в качестве предиктора. Вертикальная синяя линия показывает ИКР реально наблюдаемых значений кровяного давления (из рис. 1).

Теперь настало время раскрыть небольшой секрет: тот факт, что модель, включающая возраст, гораздо лучше описывает данные из рис. 1, неудивителен, поскольку эти данные были... сгенерированы на основе модели

\[y_i = 97.078 + 0.949\times  Age_i + \epsilon_i,\]

где \(\epsilon_i \sim N(0, 9.563),\)

следующим образом:

set.seed(101)
y <- rnorm(100, mean = 97.078 + 0.949*x, 9.563)

Эта последняя модель [придуманная в целях демонстрации принципа!] была использована в качестве "истинной" (англ. true model) в том смысле, что она описывала генеральную совокупность (т.е. если бы у нас была возможность одномоментно измерить давление у всех существующих людей, то полученные данные описывались бы именно этой "истинной" моделью). Конечно, в реальной жизни ни структура (систематическая часть + остатки), ни значения параметров истинной модели исследователю неизвестны. Все, что видит исследователь - это данные (вроде тех, что представлены на рис. 1). Имея эти данные и хорошее понимание изучаемого явления (в смысле того, какие предикторы считать достаточно важными для рассмотрения), он может только надеяться приблизиться к структуре истинной модели и оценить ее параметры с определенной точностью. В приведенном выше примере параметры, оцененные по методу наименьших квадратов, оказались очень близки к истинным, что вполне понятно, поскольку нам "посчастливилось" воспроизвести структуру истинной модели и выполнить оценку параметров по выборке достаточно большого размера. К сожалению, на практике, такой успех гарантирован далеко не всегда.


Зачем нужны статистические модели?

В общем виде истинную функцию зависимости переменной \(Y\) от набора предикторов \(X\) можно обозначить как \(Y = f(X) + \epsilon\). В сущности, статистическое моделирование (или "статистическое обучение", от англ. statistical learning, см. James et al. (2013) An Introduction to Statistical Learning. Springer) представляет собой набор подходов, позволяющих оценить \(f\) на основе ограниченной совокупности наблюдений. Существуют две основные причины для нахождения \(f\):

• Предсказание: часто у нас в распоряжении имеется набор значений предикторов (например, уровень доходов, уровень образования, опыт работы, пол, возраст, кредитная история), на основе которых мы хотели бы предсказать определенный отклик \(Y\) (например, риск невозврата кредита). Поскольку остатки модели в среднем равны 0 (см. выше), мы можем их проигнорировать и предсказать \(Y\) как \(\widehat{Y} = \hat{f}(X)\). Здесь \(\hat{f}\) - наша приближенная оценка  \(f\), а  \(\widehat{Y}\) - предсказанное значение \(Y\). Точность предсказания будет зависеть от двух составляющих - устранимая и неустранимая погрешности (англ. reducible error и irreducible error соответственно). Первая связана с тем, что наша модель \(\hat{f}\) никогда не является идеальной, но при использовании дополнительной информации об изучаемом явлении и применении более точных статистических методов мы можем ее снизить. В свою очередь, неустранимая погрешность обусловлена неточностью измерений и влиянием скрытых факторов, которые мы не можем учесть, потому что просто о них не знаем.
• Статистические выводы: в научных исследования постоянно возникает вопрос вроде "зависит ли переменная \(Y\) от факторов \(X_1\), \(X_2\), и т.д. ?" В таких случаях мы также оцениваем \(f\), но не для осуществления предсказаний, а для понимания характера взаимоотношений между переменной-откликом и предикторами (например, насколько силен отклик при изменении того или иного предиктора, какой из предикторов наиболее важен, и т.п.). Так, в примере выше мы выяснили, что возраст тесно связан с уровнем систолического давления - информация, которая может оказаться очень важной для практикующего врача.

Заключение

В этом сообщении я постарался раскрыть суть понятия "статистическая модель", но обсудил лишь некоторые аспекты построения и применения статистических моделей. Тема эта бесконечна. Например, ничего не было сказано о параметрическом и непараметрическом типах моделей, методе максимального правдоподобия как универсальном методе оценивания параметров моделей, линейных и нелинейных моделях, проверке адекватности модели (хотя эта тема частично была затронута на примере симуляций), выборе "оптимальной" модели, переобучении, перекрестной проверке, и др. Об всем этом я обязательно напишу в будущих сообщениях - подписывайтесь на рассылку, чтобы не пропустить ничего интересного :)