Автор: Владимир Шитиков
Введение
Современные исследования приобретают все более и более обобщающий и стратегический характер, а глубокая стратегия никогда не ограничивается рассмотрением какой-то одной идеи, гипотезы или модели. Принцип "множественности моделей", сформулированный еще в 1890 г. Т. Чемберленом, предполагает формирование набора альтернативных научных гипотез \(H_1, H_2, \dots, H_r\), для каждой из которых подбирается адекватная математическая модель. В итоге вместо того, чтобы находить по Фишеру соотношение вероятностей основной и нулевой гипотез \(H_0 | H_1\), оценивается относительная сила обоснованности (strength of evidence) каждой из рассматриваемых гипотетических моделей \(g_1, g_2, \dots, g_r\). Эта новая парадигма оформилась в современную методологию "Model selection and Multimodel inference" (Burnham, Anderson, 2002), которая базируется на основных принципах теории информации Кульбака-Лейблера (Kullback-Leibler, 1951) и включает ранжирование моделей с последующим формированием статистического вывода на основе этих нескольких моделей.
Пакет MuMIn (от "Multi-Model Inference"), разработанный К. Бартоном, содержит набор функций, которые, используя информационные критерии, реализуют ранжирование и отбор статистических моделей различного типа и их последующее обобщение (model averaging) с целью получения коллективного решения. Ансамбль \(G_r\) моделей, включаемых в анализ, может формироваться либо автоматически (например, как все возможные комбинации подмножеств предикторов заданной "глобальной" модели), либо задается исследователем в виде набора конкретных моделей (т.е. математически выраженных гипотез, нуждающихся в проверке). Для всех анализируемых моделей выполняется подгонка их коэффициентов по эмпирическим наборам данных и рассчитываются основные статистики и информационные критерии \(IC\) качества аппроксимации. Полный список из \(r\) моделей сортируют по уменьшению адекватности на основе заданного \(IC\) и устанавливают порог, согласно которому некоторое количество "оптимальных" моделей далее будет использовано для формирования окончательного коллективного решения. Такие коллективные решения могут быть получены с использованием как традиционных взвешенных средних, так и других специальных алгоритмов (метод Бейтса-Гренджера, бутстреп, метод "складного ножа", адаптивная регрессия и др.).
В этой статье мы рассмотрим первую часть описанной процедуры, т.е. построение ранжированного списка моделей на основе информационных критериев с использованием пакета MuMIn.
Информационные критерии
Энтропийная мера Кульбака-Лейблера (K-L) трактуется как доля потерянной информации при использовании модели \(g_i\) по сравнению с отображаемой этой моделью "полной реальностью" (full reality). Таким образом, для выбора оптимальной модели из набора \(G_r\) необходимо подобрать такой критерий, который соответствует минимуму информационных потерь. В своей ключевой работе Хиротугу Акаике (Akaike, 1973) объединил понятия информационной и статистической теории и пришел к выводу фундаментальной важности: поиск оптимальной модели можно выполнить путем минимизации величины
\[\log(L(\hat{\theta}|data)) - Const,\] где \(L(\hat{\theta}|data)\) - максимум функции правдоподобия модели \(g_i\) с коэффициентами \(\hat{\theta}\), подогнанными по выборочным данным \(data\), а \(Const\) - некоторая константа, корректирующая смещение и зависящая от числа степеней свободы. Акаике предложил конкретное определение соответствующего информационного критерия: \[ AIC = -2 \log(L(\hat{\theta}|data)) + 2K, \] где K - количество параметров модели.
Далее последовали достаточно многочисленные попытки "подправить" выражение для \(Const\). В частности, для малых выборок при \(n/K < 40\) рекомендуется использовать скорректированный \(AIC\), который имеет вид (Sugiura, 1978): \[ AIC_{c} = AIC + 2K(K + 1) / (n - K - 1). \] Если объем выборки \(n \gg K\), то \(AIC_{c}\) сводится к \(AIC\). В случае обычной регрессии, решаемой методом наименьших квадратов, \(\log(L(\hat{\theta}|data)) = -(n/2) \log(RSS/n)\) и \[AIC_{c} = n \log(RSS/n) + 2K + 2K(K + 1) / (n - K - 1), \] где \(RSS\) - сумма квадратов остатков подогнанной модели.
Весьма популярен байесовский информационный критерий, который позволяет выбирать модели на основе принципа максимального апостериорного правдоподобия (a posteriori most likely). При аппроксимации формулы Байеса, выполненной Шварцем (Schwarz, 1978), этот критерий имеет вид: \[BIC = -2 \log(L(\hat{\theta}|data)) + K \log(n).\] В литературе (Claeskens, Hjort, 2008) подробно описаны и другие версии информационных критериев: Takeuchi's IC (\(TIC\)), Deviance IC (\(DIC\)), Hannan and Quinn criterion (\(HQ\)), Focused IC (\(FIC\)), информационные критерии для счетных данных с высокой дисперсией (\(QAIC, QAIC_{c}\)), \(CAICF\) (где \(C\) соответствует слову "Consistent", а \(F\) означает, что используется информационная матрица Фишера), а также близкие по сути статистики \(C_p\) Мэллоуза (Mallows) и \(ICOMP\) (Informational Complexity).
Для вычисления перечисленных критериев в пакетах stats, MuMIn, arm и GEE имеются соответствующие функции. Подлежащий использованию тип критерия выбирается из сложно формализуемых соображений и задается в представленных ниже примерах параметром rank = <имя IC-функции>, причем по умолчанию используется \(AIC_{c}\).
Рассмотрим процесс построения и селекции моделей на трех примерах с разной степенью их детализации.
\[\log(L(\hat{\theta}|data)) - Const,\] где \(L(\hat{\theta}|data)\) - максимум функции правдоподобия модели \(g_i\) с коэффициентами \(\hat{\theta}\), подогнанными по выборочным данным \(data\), а \(Const\) - некоторая константа, корректирующая смещение и зависящая от числа степеней свободы. Акаике предложил конкретное определение соответствующего информационного критерия: \[ AIC = -2 \log(L(\hat{\theta}|data)) + 2K, \] где K - количество параметров модели.
Далее последовали достаточно многочисленные попытки "подправить" выражение для \(Const\). В частности, для малых выборок при \(n/K < 40\) рекомендуется использовать скорректированный \(AIC\), который имеет вид (Sugiura, 1978): \[ AIC_{c} = AIC + 2K(K + 1) / (n - K - 1). \] Если объем выборки \(n \gg K\), то \(AIC_{c}\) сводится к \(AIC\). В случае обычной регрессии, решаемой методом наименьших квадратов, \(\log(L(\hat{\theta}|data)) = -(n/2) \log(RSS/n)\) и \[AIC_{c} = n \log(RSS/n) + 2K + 2K(K + 1) / (n - K - 1), \] где \(RSS\) - сумма квадратов остатков подогнанной модели.
Весьма популярен байесовский информационный критерий, который позволяет выбирать модели на основе принципа максимального апостериорного правдоподобия (a posteriori most likely). При аппроксимации формулы Байеса, выполненной Шварцем (Schwarz, 1978), этот критерий имеет вид: \[BIC = -2 \log(L(\hat{\theta}|data)) + K \log(n).\] В литературе (Claeskens, Hjort, 2008) подробно описаны и другие версии информационных критериев: Takeuchi's IC (\(TIC\)), Deviance IC (\(DIC\)), Hannan and Quinn criterion (\(HQ\)), Focused IC (\(FIC\)), информационные критерии для счетных данных с высокой дисперсией (\(QAIC, QAIC_{c}\)), \(CAICF\) (где \(C\) соответствует слову "Consistent", а \(F\) означает, что используется информационная матрица Фишера), а также близкие по сути статистики \(C_p\) Мэллоуза (Mallows) и \(ICOMP\) (Informational Complexity).
Для вычисления перечисленных критериев в пакетах stats, MuMIn, arm и GEE имеются соответствующие функции. Подлежащий использованию тип критерия выбирается из сложно формализуемых соображений и задается в представленных ниже примерах параметром rank = <имя IC-функции>, причем по умолчанию используется \(AIC_{c}\).
Рассмотрим процесс построения и селекции моделей на трех примерах с разной степенью их детализации.
Теплота гидратации цемента
Сначала обратимся к хрестоматийному примеру, который основан на наборе данных Cement и подробно комментируется в основной литературе (см. Burnham, Anderson, 2002, стр. 100). Задача заключается в нахождении зависимости теплоты гидратации y (кал/г), выделяющейся при затвердевании различных марок цемента, от химического состава этих марок. Четыре независимых переменных соответствуют содержанию различных оксидных составляющих: X1: \(3 CaO \times Al_2O_3\), X2: \(3CaO \times SiO_2\), X3: \(4CaO \times Al_2O_3 \times Fe_2O_3\), и X4: \(2CaO \times SiO_2\). Построим линейную регрессионную модель:library(MuMIn) # можно изменить значение по умолчанию "na.omit", чтобы # не учитывать пропущенные значения options(na.action = "na.fail") # осуществляем построение глобальной модели fm1 <- lm(y ~ ., data = Cement) summary(fm1) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 62.4054 70.0710 0.891 0.3991 X1 1.5511 0.7448 2.083 0.0708 . X2 0.5102 0.7238 0.705 0.5009 X3 0.1019 0.7547 0.135 0.8959 X4 -0.1441 0.7091 -0.203 0.8441 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.446 on 8 degrees of freedom Multiple R-squared: 0.9824, Adjusted R-squared: 0.9736 F-statistic: 111.5 on 4 and 8 DF, p-value: 4.756e-07
Функция dredge() осуществляет построение всех возможных моделей из различных комбинаций независимых переменных полной линейной модели. В данном примере общее количество таких моделей составляет \(r = 2^K = 16.\)
(ms1 <- dredge(fm1)) Model selection table (Intrc) X1 X2 X3 X4 df logLik AICc delta weight 4 52.58 1.468 0.6623 4 -28.156 69.3 0.00 0.566 12 71.65 1.452 0.4161 -0.2365 5 -26.933 72.4 3.13 0.119 8 48.19 1.696 0.6569 0.2500 5 -26.952 72.5 3.16 0.116 10 103.10 1.440 -0.6140 4 -29.817 72.6 3.32 0.107 14 111.70 1.052 -0.4100 -0.6428 5 -27.310 73.2 3.88 0.081 15 203.60 -0.9234 -1.4480 -1.5570 5 -29.734 78.0 8.73 0.007 16 62.41 1.551 0.5102 0.1019 -0.1441 6 -26.918 79.8 10.52 0.003 13 131.30 -1.2000 -0.7246 4 -35.372 83.7 14.43 0.000 7 72.07 0.7313 -1.0080 4 -40.965 94.9 25.62 0.000 9 117.60 -0.7382 3 -45.872 100.4 31.10 0.000 3 57.42 0.7891 3 -46.035 100.7 31.42 0.000 11 94.16 0.3109 -0.4569 4 -45.761 104.5 35.21 0.000 2 81.48 1.869 3 -48.206 105.1 35.77 0.000 6 72.35 2.312 0.4945 4 -48.005 109.0 39.70 0.000 5 110.20 -1.2560 3 -50.980 110.6 41.31 0.000 1 95.42 2 -53.168 111.5 42.22 0.000 Models ranked by AICc(x)
В полученной таблице представлены коэффициенты каждой из построенных моделей, значения логарифмов максимального правдоподобия и скорректированные значения информационного критерия Акаике, по убыванию которых выполнялась сортировка моделей. Важным показателем является величина \(\Delta_{i} = AIC_{c_i} - AIC_{c_{min}}\), которая соответствует потере информации Кульбака-Лейблера или, эквивалентно, увеличению расстояния от текущей до оптимальной модели. Шкала \(\Delta\) дает возможность исследователю выбрать допустимый порог информационных потерь и ограничить подмножество отбираемых моделей-претендентов по некоторому уровню адекватности.
Применительно к имеющимся данным, относительное правдоподобие \(l_i\), оцененное для каждой модели \(g_i\), соответствует формализованному уровню ее обоснованности (strength of evidence) по сравнению с наилучшей моделью и легко вычисляется как \(l_i = \exp(-0.5 \Delta_{i})\). Это выражение позволяет сформировать хорошо интерпретируемые утверждения, такие как "Модель \(g_1\), основанная на X1 и X2, в 192 раза лучше объясняет результаты эксперимента, чем полная модель \(g_{16}\)" (или "гипотеза \(H_4\) в 192 раза более вероятна, чем \(H_{16}\)"). Нетрудно вычислить, что соотношение обоснованности гипотез для \(\Delta = 2\) равно 2.7, а при \(\Delta = 4\) оно достигает 7.4. Отметим, что в рамках обсуждаемой методологии не рекомендуется использовать термин "статистическая значимость" в том смысле, в каком он используется в рамках классической дихотомической процедуры Фишера.
Относительная вероятностная мера каждой модели \(g_i\) в ансамбле из \(G_r\) моделей, или ее вес (weight), может быть вычислена как нормированное значение силы обоснованности \[w_i = \exp(-0.5 \Delta_{i}) / \sum \exp(-0.5 \Delta_{i}) .\] Найденные веса позволяют также оценить относительную важность каждой независимой переменной как сумму весов тех моделей, в которых соответствующий предиктор присутствует:
# Вычисляем степень важности отдельных предикторов importance(ms1) X1 X2 X4 X3 Importance: 0.99 0.81 0.32 0.21 N containing models: 8 8 8 8 # Выводим на экран график таблицы сгенерированных моделей: par(mar = c(3,5,6,4)) plot(ms1, labAsExpr = TRUE)
Влияние курения на заболеваемость раком
Пакет MuMIn, разумеется, не ограничивается построением моделей обычной линейной регрессии. Список поддерживаемых моделей включает более 30 типов функций из различных пакетов (glm, gam, lmer, polr, MCMCglmm, survreg, zeroinfl, betareg и проч.). Рассмотрим пример селекции моделей логистической регрессии с учетом парных взаимодействий факторов.Набор данных esoph включает результаты обследования различных групп населения одного из департаментов Франции на предмет заболеваемости раком пищевода. Три порядковых переменных определяли интенсивность приема алкоголя (alcgp) и табакокурения (tobgp) в баллах от 1 до 4, а также градацию возраста обследуемых (agegp) от 1 до 6. Переменные ncases, ncontrols определяли число случаев заболевания раком и объем контрольной группы соответственно.
Хорошим способом познакомиться со структурой данных является отображение ее в виде мозаичной диаграммы:
## Реорганизация данных для мозаичной диаграммы ttt <- table(esoph$agegp, esoph$alcgp, esoph$tobgp) ttt[ttt == 1] <- esoph$ncases tt1 <- table(esoph$agegp, esoph$alcgp, esoph$tobgp) tt1[tt1 == 1] <- esoph$ncontrols tt <- array(c(ttt, tt1), c(dim(ttt),2), c(dimnames(ttt), list(c("Болен", "Здоров")))) mosaicplot(tt, main="Заболеваемость раком", color = TRUE)
str(esoph) 'data.frame': 88 obs. of 5 variables: $ agegp : Ord.factor w/ 6 levels "25-34"<"35-44"<..: 1 1 1 1 1 1 1 1 1 1 ... $ alcgp : Ord.factor w/ 4 levels "0-39g/day"<"40-79"<..: 1 1 1 1 2 2 2 2 3 3 ... $ tobgp : Ord.factor w/ 4 levels "0-9g/day"<"10-19"<..: 1 2 3 4 1 2 3 4 1 2 ... $ ncases : num 0 0 0 0 0 0 0 0 0 0 ... $ ncontrols: num 40 10 6 5 27 7 4 7 2 1 ... fm2 <- glm(cbind(ncases, ncontrols) ~ agegp + tobgp + alcgp, data = esoph, family = binomial()) summary(fm2) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.77997 0.19796 -8.992 < 2e-16 *** agegp.L 3.00534 0.65215 4.608 4.06e-06 *** agegp.Q -1.33787 0.59111 -2.263 0.02362 * agegp.C 0.15307 0.44854 0.341 0.73291 agegp^4 0.06410 0.30881 0.208 0.83556 agegp^5 -0.19363 0.19537 -0.991 0.32164 tobgp.L 0.59448 0.19422 3.061 0.00221 ** tobgp.Q 0.06537 0.18811 0.347 0.72823 tobgp.C 0.15679 0.18658 0.840 0.40071 alcgp.L 1.49185 0.19935 7.484 7.23e-14 *** alcgp.Q -0.22663 0.17952 -1.262 0.20680 alcgp.C 0.25463 0.15906 1.601 0.10942 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Null deviance: 227.241 on 87 degrees of freedom Residual deviance: 53.973 on 76 degrees of freedom AIC: 225.45
Модель показывает, насколько возрастает риск заболеть раком при перемещении в следующую по порядку градацию каждого фактора. Построим серию моделей с разными комбинациями состава потенциальных предикторов:
ms2 <- dredge(fm2)) Model selection table (Intrc) agegp alcgp tobgp df logLik AICc delta weight 8 -1.780 + + + 12 -100.727 229.6 0.00 0.798 4 -1.911 + + 9 -106.026 232.4 2.75 0.202 6 -1.951 + + 9 -133.754 287.8 58.20 0.000 7 -1.297 + + 7 -139.482 294.4 64.75 0.000 3 -1.402 + 4 -143.135 294.8 65.14 0.000 2 -2.139 + 6 -143.297 299.6 70.02 0.000 5 -1.417 + 4 -178.506 365.5 135.88 0.000 1 -1.584 1 -187.361 376.8 147.15 0.000 Models ranked by AICc(x)
Согласно полученной сводной таблице, можно установить, что оптимальной является модель со всеми тремя факторами, которая идентична исходной полной модели (это, впрочем, вполне соответствует медицинским представлениям). Здесь знаки "+" лишь показывают, какие факторы следует использовать в моделях, а конкретные значения подогнанных коэффициентов для каждой градации фактора скрыты из соображений компактности представления результатов. Полную информацию о коэффициентах моделей можно получить, выполнив команду (attr(ms2,"coefTables").
Проведем небольшой статистический эксперимент, превратив порядковые категориальные переменные alcgp и tobgp в действительные числа и построив модели, включающие полученные целочисленные переменные, а также их парные взаимодействия:
fm3 <- glm(cbind(ncases, ncontrols) ~ agegp + unclass(tobgp) * unclass(alcgp),data = esoph, family = binomial()) summary(fm3) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -4.70284 0.49575 -9.486 < 2e-16 *** agegp.L 2.92868 0.65126 4.497 6.89e-06 *** agegp.Q -1.34088 0.58964 -2.274 0.02296 * agegp.C 0.12258 0.44850 0.273 0.78462 agegp^4 0.06880 0.30787 0.223 0.82317 agegp^5 -0.22286 0.19530 -1.141 0.25382 unclass(tobgp) 0.60346 0.20014 3.015 0.00257 ** unclass(alcgp) 0.94297 0.17916 5.263 1.41e-07 *** unclass(tobgp):unclass(alcgp) -0.14099 0.07572 -1.862 0.06262 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Null deviance: 227.241 on 87 degrees of freedom Residual deviance: 55.833 on 79 degrees of freedom AIC: 221.31 (ms3 <- dredge(fm3)) Model selection table (Int) agg unc(alc) unc(tbg) unc(alc):unc(tbg) df logLik AICc delta weight 16 -4.703 + 0.9430 0.6035 -0.1410 9 -101.657 223.6 0.00 0.612 8 -4.011 + 0.6531 0.2616 8 -103.379 224.6 0.96 0.379 4 -3.570 + 0.6824 7 -108.395 232.2 8.57 0.008 6 -2.766 + 0.3293 7 -134.536 284.5 60.85 0.000 15 -4.200 0.9903 0.5544 -0.1445 4 -139.416 287.3 63.69 0.000 7 -3.491 0.6920 0.2048 3 -141.396 289.1 65.46 0.000 3 -3.165 0.7221 2 -144.846 293.8 70.21 0.000 2 -2.139 + 6 -143.297 299.6 76.01 0.000 5 -2.152 0.2957 2 -179.350 362.8 139.22 0.000 1 -1.584 1 -187.361 376.8 153.15 0.000 Models ranked by AICc(x)
Приведенные расчеты имеют отношение к обсуждению такой нетривиальной проблемы, как статистическая природа повсеместно используемых балльных показателей, которые следует интерпретировать как порядковые категориальные переменные. В частности, по типологии Стивенса операции суммы или разности к ранговым переменным неприменимы, и поэтому некорректно вычислять средний балл, применяемый, например, в фигурном катании. Тем не менее, в нашем примере модели с "деклассированными" переменными оказались существенно более эффективными по AICc-критерию, чем модели, основанных на исходных данных. Кроме того, добавление парных взаимодействий показывает, что много пьющие и курящие люди имеют более низкую вероятность заболеть раком, нежели это предсказывается моделью, включающей только аддитивные эффекты.
Содержание озона в атмосферном воздухе
Набор данных airquality, входящий в базовую комплектацию R, содержит результаты ежедневных метеорологических наблюдений с 1 мая по 31 октября в окрестностях Нью-Йорка. Он включает такие показатели, как солнечная радиация Solar.R, средняя скорость ветра Wind, температура атмосферного воздуха Temp и содержание озона Ozone, которое интерпретируется как зависимая переменная. Чтобы оценить степень нелинейности связи между переменными, выполним визуализацию основных зависимостей с использованием функций удобного пакета visreg:
data(airquality) library(visreg) airquality$MonthF <- as.factor(airquality$Month) fit <- lm(Ozone ~ Solar.R * MonthF, data=airquality) visreg(fit, "MonthF", by="Solar.R", breaks=1, points=list(cex=1, pch=16)) fit <- lm(Ozone ~ Wind * MonthF, data=airquality) visreg(fit, "MonthF", by="Wind", breaks=1, points=list(cex=1, pch=16))
library(splines) fit <- lm(Ozone ~ Solar.R +ns(Wind, df=2)*ns(Temp, df=2), data=airquality) visreg2d(fit, "Wind", "Temp", plot.type="persp", col = 3)
Нетрудно заметить, что имеющиеся предикторы не являются независимыми: между ними существуют выраженные взаимодействия, имеющие нелинейный характер. Поэтому построим две исходные регрессионные модели для прогнозирования концентрации озона, первая из которых будет учитывать неаддитивность парных взаимодействий предикторов, а вторая - нелинейность их влияния на отклик. Конечно, можно было бы построить и единую модель, включающую сразу все необходимые трансформации базовых переменных, но мы хотим показать, как работает функция merge(), выполняющая объединение двух списков частных моделей.
Перед началом расчетов выполним очистку набора данных от строк, содержащих пропущенные значения:
airquality <- airquality[complete.cases(airquality),] airquality$Month.Day <- airquality$Month + airquality$Day/31 mod1 <- lm(Ozone ~ Month.Day + (Solar.R + Wind + Temp )^2, data=airquality) summary(mod1) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -1.421e+02 6.404e+01 -2.218 0.02872 * Month.Day -1.690e+00 1.368e+00 -1.235 0.21963 Solar.R -1.978e-01 2.114e-01 -0.936 0.35158 Wind 1.056e+01 4.280e+00 2.467 0.01527 * Temp 2.527e+00 8.472e-01 2.982 0.00357 ** Solar.R:Wind -7.228e-03 6.671e-03 -1.083 0.28113 Solar.R:Temp 4.594e-03 2.468e-03 1.862 0.06549 . Wind:Temp -1.612e-01 5.881e-02 -2.741 0.00723 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 19.12 on 103 degrees of freedom Multiple R-squared: 0.6909, Adjusted R-squared: 0.6699 F-statistic: 32.89 on 7 and 103 DF, p-value: < 2.2e-16
На основе этой "полной" модели выполним построение ансамбля моделей с использоваением всех возможных комбинаций сформированных переменных. Поскольку число исходных параметров \(m = 7\), то наш ансамбль будет содержать \(2^7 = 128\) моделей. Отсортируем его по убыванию BIC-критерия и с использованием функции subset() выберем подмножество наиболее адекватных моделей с delta < 4:
ms1 <- dredge(mod1, rank=BIC) subset(ms1,delta < 4) Global model call: lm(formula = Ozone ~ Month.Day + (Solar.R + Wind + Temp)^2, data = airquality) --- Model selection table (Int) Mnt.Day Slr.R Tmp Wnd Slr.R:Tmp Slr.R:Wnd Tmp:Wnd df logLik BIC delta weight 95 -136.8 -0.35310 2.451 11.15 0.005717 -0.1863 7 -482.309 997.6 0.00 0.435 79 -245.1 0.06599 3.914 14.38 -0.2280 6 -485.303 998.9 1.28 0.230 96 -138.1 -1.690 -0.32480 2.656 11.15 0.005250 -0.1862 8 -481.502 1000.7 3.10 0.093 111 -232.6 0.18250 3.476 12.96 -0.010670 -0.1839 7 -483.930 1000.8 3.24 0.086 80 -235.4 -2.165 0.05826 4.022 14.05 -0.2235 7 -484.010 1001.0 3.40 0.079 127 -140.8 -0.22600 2.322 10.55 0.005061 -0.007231 -0.1613 8 -481.689 1001.1 3.47 0.077 Models ranked by BIC(x)
Для построения второй модели сформируем дополнительные переменные как квадраты исходных предикторов, а для отражения сезонности в данных воспользуемся функцией синусоиды для верхнего полупериода.
airquality$Solar.R.2 <- airquality$Solar.R^2 airquality$Wind.2 <- airquality$Wind^2 airquality$Temp.2 <- airquality$Temp^2 norm.0_1 <- function(x) (x-min(x))/(max(x)-min(x)) airquality$sin.period <- sin(pi*norm.0_1(airquality$Month.Day)) mod2 <- lm(Ozone ~ Solar.R + Wind + Temp + Solar.R.2 + Wind.2 + Temp.2 + sin.period, data=airquality) summary(mod2) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.541e+02 1.046e+02 3.384 0.00101 ** Solar.R 1.810e-01 9.041e-02 2.001 0.04798 * Wind -1.330e+01 2.280e+00 -5.831 6.37e-08 *** Temp -8.007e+00 2.812e+00 -2.847 0.00532 ** Solar.R.2 -3.316e-04 2.608e-04 -1.272 0.20637 Wind.2 4.599e-01 9.986e-02 4.605 1.18e-05 *** Temp.2 5.936e-02 1.805e-02 3.289 0.00137 ** sin.period 1.326e+01 7.046e+00 1.882 0.06264 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 18.08 on 103 degrees of freedom Multiple R-squared: 0.7236, Adjusted R-squared: 0.7048 F-statistic: 38.52 on 7 and 103 DF, p-value: < 2.2e-16
На основе второй "полной" модели проведем построение еще одного ансамбля, также включающего 128 моделей. С использованием функции merge() объединим оба списка моделей и получим объединенный ансамбль из 164 моделей (в обоих списках модели, состоящие из базовых переменных Solar.R + Wind + Temp, совпадают и в объединенном списке их дубли исключаются). С помощью функции subset() выберем подмножество из 17 наиболее адекватных моделей с накопленной суммой весов cumsum(weight) < 0.95 (в следующем сообщении мы покажем, что эта процедура связана с оценкой доверительных интервалов).
allmod <- merge(ms1, dredge(mod2, rank=BIC)) dim(allmod) [1] 164 17 ansmod <- subset(allmod, cumsum(weight) < 0.95) dim(ansmod) [1] 17 17
Функции dredge(), merge(), subset() и используемая ниже функция model.sel() возвращают специальный объект класса model.selection, весьма похожий на обычную таблицу класса data.frame. Мы может использовать компоненты этого объекта для вывода нужных нам блоков информации. Предварительно, разумеется, стоит ознакомиться с "внутренностями" объекта, выполнив команду str():
# Извлекаем список формул моделей ModList <- as.character(attr(ansmod,"model.calls")) # Из списка формул выделяем только предикторную часть ModList <- substring(ModList, regexpr("~", ModList)+1, regexpr(",", ModList) - 5) # Для компактности удаляем пробелы ModList <-gsub(" ", "", ModList, fixed = TRUE) # Объединяем данные data.frame(logLik=ansmod$logLik, BIC=ansmod$BIC, delta=ansmod$delta, weight=ansmod$weight, Formula=ModList) logLik BIC delta weight Formula 1 -476.8 986.7 0.0000 0.2688 Solar.R+Temp+Temp.2+Wind+Wind.2 2 -479.6 987.6 0.8917 0.1721 Solar.R+Temp.2+Wind+Wind.2 3 -477.7 988.5 1.7847 0.1101 Solar.R.2+Temp+Temp.2+Wind+Wind.2 4 -475.5 988.7 2.0011 0.0988 sin.period+Solar.R+Temp+Temp.2+Wind+Wind.2 5 -480.2 988.8 2.0989 0.0941 Solar.R.2+Temp.2+Wind+Wind.2 6 -481.0 990.4 3.6889 0.0425 Solar.R+Temp+Wind+Wind.2 7 -476.5 990.7 4.0253 0.0359 Solar.R+Solar.R.2+Temp+Temp.2+Wind+Wind.2 8 -476.7 991.2 4.5070 0.0282 sin.period+Solar.R.2+Temp+Temp.2+Wind+Wind.2 9 -479.1 991.3 4.5954 0.0270 sin.period+Solar.R+Temp.2+Wind+Wind.2 10 -481.6 991.6 4.8897 0.0233 Solar.R.2+Temp+Wind+Wind.2 11 -474.6 991.7 4.9815 0.0222 sin.period+Solar.R+Solar.R.2+Temp+Temp.2+Wind+Wind.2 12 -484.1 991.9 5.1737 0.0202 Temp.2+Wind+Wind.2 13 -479.5 992.0 5.3085 0.0189 Solar.R+Solar.R.2+Temp.2+Wind+Wind.2 14 -479.9 992.8 6.1169 0.0126 sin.period+Solar.R.2+Temp.2+Wind+Wind.2 15 -482.3 992.8 6.1473 0.0124 Temp+Temp.2+Wind+Wind.2 16 -480.5 994.0 7.3394 0.0068 sin.period+Solar.R+Temp+Wind+Wind.2 17 -485.4 994.4 7.7390 0.0056 Temp+Wind+Wind.2
Отметим, что, рассматривая информационные потери по Кульбаку, полученная выше оптимальная модель в \(\exp((1012 - 986.7)/2) = 347000 \) раз лучше, чем линейная модель с тремя базовыми предикторами, и в \( exp((1101 - 986.7)/2) = 8 \times 10^{24}\) раз эффективнее по сравнению с нуль-моделью без параметров (К. Бёрнхэм и Д. Андерсон считают, что в том же соотношении находятся и вероятности соответствующих статистических гипотез). Далее заметим, что в окончательный список вошли только предискторы из второй модели, а относительная важность отдельных предикторов приобрела следующие значения:
importance(ansmod)
Wind Wind.2 Temp.2 Solar.R Temp Solar.R.2 sin.period
Importance: 1.00 1.00 0.92 0.69 0.65 0.35 0.20
N containing models: 17 17 13 9 11 8 6
Для полноты изложения покажем, как к сформированному ансамблю ansmod добавить обобщенную аддитивную модель, в которой варьирование трех основных предикторов аппроксимируется в виде сглаживающих p-сплайнов:
library(mgcv) mgam <- gam(Ozone ~ s(Solar.R) + s(Wind) + s(Temp), data = airquality) summary(mgam) Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 42.099 1.663 25.32 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Approximate significance of smooth terms: edf Ref.df F p-value s(Solar.R) 2.760 3.447 3.967 0.00769 ** s(Wind) 2.910 3.657 13.768 1.52e-08 *** s(Temp) 3.833 4.753 11.769 7.40e-09 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 R-sq.(adj) = 0.723 Deviance explained = 74.7% GCV score = 338.9 Scale est. = 306.83 n = 111 BIC(mgam) [1] 993.7625 gansmod <- merge(ansmod, model.sel(mgam, rank=BIC))
К сожалению, без тщательного подбора характеристик сглаживающих функций модель GAM продемонстрировала мало впечатляющую эффективность и в списке моделей заняла скромное 16 место. Однако наша задача заключалась лишь в иллюстрации работы функции model.sel(), которая создает объект model.selection, объединяя произвольное количество частных моделей.
Рекомендуемая литература
- Akaike H (1973) Information theory as an extension of the maximum likelihood principle. In: Petrov BN, Csaki F (eds) Second international symposium on information theory. Akademiai Kiado, Budapest, pp 267-281
- Burnham KP, Anderson DR (2002) Model selection and multimodel inference: a practical information-theoretic approach, 2nd edn. Springer, New York
- Claeskens G, Hjort NL (2008) Model Selection and Model Averaging. Cambridge University Press, New York
- Kullback S, Leibler RA (1951) On information and sufficiency. Ann Math Stat 22: 79-86
Отправить комментарий