В предыдущем сообщении был приведен пример оценки параметров простой линейной регрессии при помощи стандартной R-функции lm(). Как и в случае с любыми другими выборочными оценками, всегда существует неопределенность в отношении того, насколько выборочные оценки параметров регрессионной модели близки к соответствующим истинным значениям (т.е. в генеральной совокупности). В данном сообщении будут рассмотрены несколько способов, позволяющих охарактеризовать эту неопределенность.
Стандартные ошибки и доверительные интервалы параметров модели, рассчитанные в соответствии с Центральной Предельной Теоремой
Согласно Центральной Предельной Теореме, выборочные оценки того или иного параметра линейной модели имеют приближенно нормальное распределение (при условии, что объем выборки "достаточно велик" и выполняются определенные условия в отношении свойств остатков модели; рассмотрению этих свойств будет посвящено отдельное сообщение). Помимо оценки значения самого параметра модели, функция lm() рассчитывает также и соответствующую стандартную ошибку. По определению, стандартная ошибка параметра - это стандартное отклонение [нормального] распределения значений этого параметра, рассчитанных по выборкам одинакового размера из той же генеральной совокупности. Соответственно, чем меньше значение стандартной ошибки параметра, тем более он точен.
В предыдущем сообщении была рассмотрена простая регрессионная модель вида \(y_i = \beta x_i + \epsilon_i\) для данных о скорости движения галактик, где \(\beta\) - подлежавшая оценке постоянная Хаббла. Подгонка этой модели к соответствующим данным дала следующие результаты:
В предыдущем сообщении была рассмотрена простая регрессионная модель вида \(y_i = \beta x_i + \epsilon_i\) для данных о скорости движения галактик, где \(\beta\) - подлежавшая оценке постоянная Хаббла. Подгонка этой модели к соответствующим данным дала следующие результаты:
library(gamair) data(hubble) M <- lm(y ~ x - 1, data = hubble) summary(M) Call: lm(formula = y ~ x - 1, data = hubble) Residuals: Min 1Q Median 3Q Max -736.5 -132.5 -19.0 172.2 558.0 Coefficients: Estimate Std. Error t value Pr(>|t|) x 76.581 3.965 19.32 1.03e-15 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 258.9 on 23 degrees of freedom Multiple R-squared: 0.9419, Adjusted R-squared: 0.9394 F-statistic: 373.1 on 1 and 23 DF, p-value: 1.032e-15
Как видно из представленных результатов, постоянная Хаббла была оценена (Estimate) в 76.581 км/с на мегапарсек. Стандартная ошибка (Std. Error) этого значения составила 3.965 км/с на мегапарсек. Другими словами, мы можем сказать, что истинное значение постоянной Хаббла принадлежит нормально распределенной совокупности, математическое ожидание и стандартное отклонение которой составляют около 76.581 и 3.965 км/с на мегапарсек соответственно. Заметьте, однако, следующее обстоятельство: хотя наиболее вероятное истинное значение постоянной Хаббла составляет около 77 км/с на мегапарсек (после округления значения 76.581), возможны также и другие значения. Чтобы охарактеризовать эту неопределенность, мы можем рассчитать доверительный интервал - непротиворечащий имеющимся данным диапазон значений, в котором истинное значение постоянной Хаббла находится с определенной вероятностью (например, 95%).
Поскольку, как уже было отмечено выше, выборочная оценка постоянной Хаббла в лимите имеет нормальное распределение, мы могли бы рассчитать 95%-ный доверительный интервал, вспомнив, что примерно 95% всех значений нормального распределения лежат в диапазоне +/- 2 стандартных отклонения относительно среднего значения. Для нашего примера получаем: \(76.581 - 2 \times 3.965 = 68.651\) и \(76.581 + 2 \times 3.965 = 84.511\) км/с на мегапарсек. Однако неопределенность имеется в отношении не только оценки постоянной Хаббла, но и оценки стандартного отклонения соответствующего нормального распределения. Не углубляясь в детали, отметим, что в связи с эти обстоятельством более точные значения границ доверительного интервала дадут вычисления, основанные на свойствах не нормального, а t-распределения Стьюдента. Так, границы 95%-го доверительного интервала для параметра \(\beta\) составят:
\[\beta \pm t_{0.975} SE_{\beta},\]
где \(t_{0.975}\) - 0.975-квантиль t-распределения с \(n - p\) числом степеней свободы (где \(n\) - объем выборки, а \(p\) - число параметров модели), \(SE_{\beta}\) - стандартная ошибка параметра \(\beta\).
Поскольку, как уже было отмечено выше, выборочная оценка постоянной Хаббла в лимите имеет нормальное распределение, мы могли бы рассчитать 95%-ный доверительный интервал, вспомнив, что примерно 95% всех значений нормального распределения лежат в диапазоне +/- 2 стандартных отклонения относительно среднего значения. Для нашего примера получаем: \(76.581 - 2 \times 3.965 = 68.651\) и \(76.581 + 2 \times 3.965 = 84.511\) км/с на мегапарсек. Однако неопределенность имеется в отношении не только оценки постоянной Хаббла, но и оценки стандартного отклонения соответствующего нормального распределения. Не углубляясь в детали, отметим, что в связи с эти обстоятельством более точные значения границ доверительного интервала дадут вычисления, основанные на свойствах не нормального, а t-распределения Стьюдента. Так, границы 95%-го доверительного интервала для параметра \(\beta\) составят:
\[\beta \pm t_{0.975} SE_{\beta},\]
где \(t_{0.975}\) - 0.975-квантиль t-распределения с \(n - p\) числом степеней свободы (где \(n\) - объем выборки, а \(p\) - число параметров модели), \(SE_{\beta}\) - стандартная ошибка параметра \(\beta\).
Получаем:
Используя эти оценки нижней и верхней границ 95%-ного доверительного интервала, мы можем рассчитать диапазон значений, в котором с вероятностью 95% находится истинный возраст Вселенной:
Uni.upper <- 1/(ci.lower*60^2*24*365.25/3.09e19) Uni.lower <- 1/(ci.upper*60^2*24*365.25/3.09e19) Uni.lower [1] 11549039573 # т.е. примерно 11.5 млрд. лет Uni.upper [1] 14319551383 # т.е. примерно 14.3 млрд. лет
Заметьте, что полученный таким способом доверительный интервал включает значение возраста Вселенной, рассчитанное на основе самых последних космологических данных (13.798 млрд. лет).
Стандартные ошибки и доверительные интервалы параметров модели, рассчитанные бутстреп-методом
Насколько мне известно, устоявшегося перевода термина "bootstrap" с английского языка на русский не существует. Используются разные варианты: "бутстреп", "бутстрэп", "бутстрап", "размножение выборок", "метод псевдовыборок" и даже "ресамплинг" (от англ. "resampling"). Несмотря на сложности с русскоязычным названием, суть метода, тем менее, весьма проста (подробнее см. оригинальную работу Efron 1979). Предположим, что у нас есть выборка некоторого ограниченного объема и мы имеем основания полагать, что эта выборка является репрезентативной (т.е. хорошо отражает свойства генеральной совокупности, из которой она была взята). Идея бутстреп-метода заключается в том, что мы можем рассматривать саму эту выборку в качестве "генеральной совокупности" и, соответственно, можем извлечь большое число случайных выборок из этой исходной совокупности для расчета интересующего нас параметра (или параметров). Очевидно, что благодаря случайному формированию этих новых выборок, будет наблюдаться определенная вариация значений оцениваемого параметра. Другими словами, мы получим некоторое распределение значений этого параметра. Рассчитав стандартное отклонение этого распределения, мы получим оценку стандартной ошибки параметра, которая при большом числе наблюдений будет асимптотически приближаться к истинной стандартной ошибке. Аналогично, найдя, например, 0.025- и 0.975-квантили этого распределения, мы получим оценки нижней и верхней границ 95%-ного доверительного интервала.
Применим бутстреп-метод для оценки стандартной ошибки и 95%-ного доверительного интервала для постоянной Хаббла из рассмотренной выше задачи. На рисунке ниже приведены примеры 6 случайных выборок объемом 24 наблюдения каждая, извлеченных из исходной совокупности данных. Эти выборки формируются "с возвратом" - это значит, что однажды попав в новую выборку, то или иное наблюдение "возвращают" в исходную совокупность и оно может быть выбрано снова (следовательно, определенные наблюдения могут повторяться в новой выборке несколько раз).
Как отмечено выше, при выполнении бутстреп-анализа формируется "большое" число повторных выборок из исходной совокупности. То, насколько "большим" должно быть это число, зависит от объема исходной совокупности и ее свойств (подробнее см. Efron and Tibshirani 1994). Для нашего примера мы сформируем 1000 случайных повторных выборок. Хотя бутстреп-регрессию можно реализовать в R, "вручную" прописав необходимые команды, проще будет воспользоваться специально предназначенной для этого функцией boot() из одноименного стандартного пакета. Эта функция имеет следующие обязательные аргументы:
- data: таблица с исходными данными;
- statistic: функция, выполняющая вычисление интересующего нас параметра(-ов);
- R: число повторных выборок, по которым рассчитывается этот параметр.
Теперь подадим regr() на функцию boot():
При выводе содержимого объекта results на экран получим:
В представленных результатах original - это значение постоянной Хаббла, оцененное по исходным данным (см. выше); bias (смещение) - разница между средним значением 1000 бутстреп-оценок постоянной Хаббла и исходной оценкой; std. error - бутстреп-оценка стандартной ошибки постоянной Хаббла. Обратите внимание на то, что полученная бутстреп-оценка стандартной ошибки выше, чем рассчитанная по исходным данным.
Для объектов класса boot (к которому принадлежит results) имеется метод plot, который позволит изобразить графически полученное распределение бутстреп-оценок постоянной Хаббла и одновременно проверить нормальность их распределения при помощи графика квантилей.
plot(results)
Как отмечено выше, мы можем оценить нижнюю и верхнюю границы 95%-ного доверительного интервала постоянной Хаббла, найдя 0.025- и 0.975-квантили изображенного выше распределения:
Используя эти значения, рассчитаем соответствующие границы значений возраста Вселенной:
U.lower <- 1/(85.73249*60^2*24*365.25/3.09e19) U.upper <- 1/(67.07360*60^2*24*365.25/3.09e19) U.lower [1] 11421129999 # т.е. примерно 11.4 млрд. лет U.upper [1] 14598320553 # т.е. примерно 14.6 млрд. лет
Следует отметить, что рассчитанные таким образом границы доверительного интервала могут оказаться смещенными. В состав пакета boot входит функция boot.ci(), которая позволяет рассчитать несколько других типов доверительных интервалов, включая интервалы с поправкой на смещение (аргумент type = "bca"; подробнее см. ?boot.ci):
Подробнее о реализации бутстреп-регрессии средствами R можно прочитать в работе Fox (2012).
Стандартные ошибки и доверительные интервалы параметров модели, рассчитанные путем симуляций
Как было показано ранее в сообщении, посвященном понятию "статистические модели", мы можем также использовать принципы байесовской статистики и оценить неопределенность в отношении параметров той или иной модели путем симуляций. В общем виде эта процедура включает следующие шаги:
- Подгонка определенной модели к имеющимся данным;
- Использование полученных оценок параметров модели для симуляции большого количества альтернативных, но в то же время правдоподобных (т.е. непротиворечащих данным) реализаций этих параметров.
Соответствующие симуляции легко выполнить при помощи функции sim() из пакета arm. Алгоритм ее работы отражает описанные выше шаги (подробнее см. Gelman & Hill 2007):
- При помощи обычного регрессионного анализа оцениваются вектор регрессионных коэффициентов \(\hat{\beta}\), ковариационная матрица \(V_{\hat{\beta}}\) и стандартное отклонение остатков \(\hat{\sigma}\).
- Создается большое число реализаций вектора \(\beta\) и стандартного отклонения остатков \(\sigma\). Для каждой реализации:
- симулируется значение \(\sigma = \hat{\sigma}\sqrt{(n-k)/X}\), где \(X\) - это значение статистики хи-квадрат, случайным образом извлеченное из соответствующего распределения с \(n-k\) степенями свободы.
- с учетом полученного значения \(\sigma\), симулируется вектор \(\beta\) путем случайного извлечения значений из многомерного нормального распределения со средним значением \(\hat{\beta}\) и ковариационной матрицей \(V_{\hat{\beta}}\).
Для нашего примера получаем:
Как видим, эта стандартная ошибка (как и в случае с бутстреп-анализом) также оказалась выше исходной оценки (см. summary(M)).
95%-ный доверительный интервал легко нахойти уже известным нам способом:
Сравните эти значения со значениями, полученными выше другими методами.
95%-ный доверительный интервал легко нахойти уже известным нам способом:
Сравните эти значения со значениями, полученными выше другими методами.
Заключение: Конечно, может возникнуть вопрос - какой из описанных выше методов следует применять на практике? Ответ на этот вопрос определяется несколькими факторами. Так, если модель адекватно описывает данные (с точки зрения ее структуры и выполнения требований к ее остаткам), то рассмотренный в самом начале классический подход даст хороший результат. Бутстреп-анализ и симуляции, возможно, окажутся слишком "тяжелой артилерией" для характеристики неопределенности в отношении одного-единственного кэффициента простой регрессии. В то же время эти методы незаменимы при оценке неопределенности в отношении величин, которые являются функцией (особенно нелинейной функцией) от каких-либо параметров модели (например, при выполнении бутстреп-анализа мы могли бы рассчитывать для каждой "псевдовыборки" не просто постоянную Хаббла, а сразу возраст Вселенной, изменив соответствующим образом функцию regr()).
Давно ищу информацию, которая доступно излагала бы суть Бутсрепа и Байасовских методов. Очень рад, что Вы дошли до них (хоть и вскользь) в своей просветительской работе.
Хочу узнать, не собираетесь ли Вы в будущем посвятить отдельные сообщения по этим темам? Было бы очень неплохо:)
Другие работы на русском языке, доступнее этой, найти вряд ли удастся (но могу, конечно, и ошибаться).
Отправить комментарий