14 января 2018

Работа с пакетом MuMIn. Часть 1: Селекция моделей



Автор: Владимир Шитиков

Введение

Современные исследования приобретают все более и более обобщающий и стратегический характер, а глубокая стратегия никогда не ограничивается рассмотрением какой-то одной идеи, гипотезы или модели. Принцип "множественности моделей", сформулированный еще в 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\) - сумма квадратов остатков подогнанной модели.

Весьма популярен байесовский информационный критерий, который позволяет выбирать модели на основе принципа максимального апостериорного правдоподобия (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

Комментариев нет :

Отправить комментарий