12 июня 2013

Дисперсионный анализ: структура модельных объектов



При подгонке линейных и других моделей в R мы, как правило, сохраняем результаты вычислений в виде самостоятельного объекта - модельного объекта (англ. model object). С точки зрения представления данных в нем, такой объект обычно является списком. Соответственно, при необходимости мы без труда можем извлечь отдельные элементы из такого списка и использовать их в последующих вычислениях. Такой подход к выводу и хранению результатов подгонки моделей является одним из существенных преимуществ R в сравнении с большинством других статистических программ, которые просто выводят результаты вычислений на экран, без возможности их непосредственного использования в дальнейшем анализе. Рассмотрим, как устроены модельные объекты в R, и как извлекать из их отдельные элементы.

В качестве примера используем данные по эффективности инсектицидов и применим однофакторный дисперсионный анализ для проверки нулевой гипотезы об отсутствии разницы между этими препаратами. Как уже отмечалось ранее, однофакторный дисперсионный анализ в R можно выполнить двумя способами - при помощи функции aov() и, в более общем виде, при помощи функции lm(). Начнем с aov():

M1 <- aov(count ~ spray, data = InsectSprays)

В ходе выполнения приведенной команды внешне ничего не происходит. На самом же деле создается достаточно сложно устроенный объект (M1), строение которого мы можем выяснить при помощи функции str() (от structure - структура):

str(M1)
List of 13
 $ coefficients : Named num [1:6] 14.5 0.833 -12.417 -9.583 -11 ...
  ..- attr(*, "names")= chr [1:6] "(Intercept)" "sprayB" "sprayC" "sprayD" ...
 $ residuals    : Named num [1:72] -4.5 -7.5 5.5 -0.5 -0.5 ...
  ..- attr(*, "names")= chr [1:72] "1" "2" "3" "4" ...
 $ effects      : Named num [1:72] -80.6 22.1 -24.2 -19.9 -34.2 ...
  ..- attr(*, "names")= chr [1:72] "(Intercept)" "sprayB" "sprayC" "sprayD" ...
 $ rank         : int 6
 $ fitted.values: Named num [1:72] 14.5 14.5 14.5 14.5 14.5 ...
  ..- attr(*, "names")= chr [1:72] "1" "2" "3" "4" ...
 $ assign       : int [1:6] 0 1 1 1 1 1
 $ qr           :List of 5
  ..$ qr   : num [1:72, 1:6] -8.485 0.118 0.118 0.118 0.118 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:72] "1" "2" "3" "4" ...
  .. .. ..$ : chr [1:6] "(Intercept)" "sprayB" "sprayC" "sprayD" ...
  .. ..- attr(*, "assign")= int [1:6] 0 1 1 1 1 1
  .. ..- attr(*, "contrasts")=List of 1
  .. .. ..$ spray: chr "contr.treatment"
  ..$ qraux: num [1:6] 1.12 1.05 1.06 1.07 1.09 ...
  ..$ pivot: int [1:6] 1 2 3 4 5 6
  ..$ tol  : num 1e-07
  ..$ rank : int 6
  ..- attr(*, "class")= chr "qr"
 $ df.residual  : int 66
 $ contrasts    :List of 1
  ..$ spray: chr "contr.treatment"
 $ xlevels      :List of 1
  ..$ spray: chr [1:6] "A" "B" "C" "D" ...
 $ call         : language aov(formula = count ~ spray, data = InsectSprays)
 $ terms        :Classes 'terms', 'formula' length 3 count ~ spray
  .. ..- attr(*, "variables")= language list(count, spray)
  .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:2] "count" "spray"
  .. .. .. ..$ : chr "spray"
  .. ..- attr(*, "term.labels")= chr "spray"
  .. ..- attr(*, "order")= int 1
  .. ..- attr(*, "intercept")= int 1
  .. ..- attr(*, "response")= int 1
  .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..- attr(*, "predvars")= language list(count, spray)
  .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "factor"
  .. .. ..- attr(*, "names")= chr [1:2] "count" "spray"
 $ model        :'data.frame': 72 obs. of  2 variables:
  ..$ count: num [1:72] 10 7 20 14 14 12 10 23 17 20 ...
  ..$ spray: Factor w/ 6 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 1 1 ...
  ..- attr(*, "terms")=Classes 'terms', 'formula' length 3 count ~ spray
  .. .. ..- attr(*, "variables")= language list(count, spray)
  .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:2] "count" "spray"
  .. .. .. .. ..$ : chr "spray"
  .. .. ..- attr(*, "term.labels")= chr "spray"
  .. .. ..- attr(*, "order")= int 1
  .. .. ..- attr(*, "intercept")= int 1
  .. .. ..- attr(*, "response")= int 1
  .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. .. ..- attr(*, "predvars")= language list(count, spray)
  .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "factor"
  .. .. .. ..- attr(*, "names")= chr [1:2] "count" "spray"
 - attr(*, "class")= chr [1:2] "aov" "lm"


От выведенных результатов может закружиться голова, но переживать не стоит - я не стану объяснять, что значит каждый из элементов объекта M1. Главное то, что этот объект, как  хорошо видно, представляет собой список из 13 основных элементов (List of 13). Некоторые из этих элементов сами являются списками - отсюда такая сложная структура. Тем не менее, используя стандартные для R приемы индексации списков, мы можем легко извлечь любой желаемый элемент объекта M1. На практике наиболее интересными и важными для анализа являются, например, такие элементы, как coefficients (коэффициенты подогнанной линейной модели), residuals (остатки), fitted.values (предсказанные моделью значения), и некоторые другие. Используя оператор $, мы можем, например, извлечь остатки модели и предсказанные ею значения и сохранить их в виде самостоятельных векторов:

M1.res <- M1$residuals
M1.fit <- M1$fitted.values

Эти два новых вектора мы можем использовать далее для других целей. Например, хранящиеся в них значения можно отобразить графически для проверки условия гетероскедастичности в отношении модели M1:

plot(M1.fit, M1.res, pch = 19, col = 4, 
     xlab = "Предсказанные значения",
     ylab = "Остатки")



Аналогично сложно устроенный объект будет создан при подгонке линейной модели по данным InsectSpays при помощи функции lm():

M2 <- lm(count ~ spray, data = InsectSprays)
 
str(M2)
 
List of 13
 $ coefficients : Named num [1:6] 14.5 0.833 -12.417 -9.583 -11 ...
  ..- attr(*, "names")= chr [1:6] "(Intercept)" "sprayB" "sprayC" "sprayD" ...
 $ residuals    : Named num [1:72] -4.5 -7.5 5.5 -0.5 -0.5 ...
  ..- attr(*, "names")= chr [1:72] "1" "2" "3" "4" ...
 $ effects      : Named num [1:72] -80.6 22.1 -24.2 -19.9 -34.2 ...
  ..- attr(*, "names")= chr [1:72] "(Intercept)" "sprayB" "sprayC" "sprayD" ...
 $ rank         : int 6
 $ fitted.values: Named num [1:72] 14.5 14.5 14.5 14.5 14.5 ...
  ..- attr(*, "names")= chr [1:72] "1" "2" "3" "4" ...
 $ assign       : int [1:6] 0 1 1 1 1 1
 $ qr           :List of 5
  ..$ qr   : num [1:72, 1:6] -8.485 0.118 0.118 0.118 0.118 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:72] "1" "2" "3" "4" ...
  .. .. ..$ : chr [1:6] "(Intercept)" "sprayB" "sprayC" "sprayD" ...
  .. ..- attr(*, "assign")= int [1:6] 0 1 1 1 1 1
  .. ..- attr(*, "contrasts")=List of 1
  .. .. ..$ spray: chr "contr.treatment"
  ..$ qraux: num [1:6] 1.12 1.05 1.06 1.07 1.09 ...
  ..$ pivot: int [1:6] 1 2 3 4 5 6
  ..$ tol  : num 1e-07
  ..$ rank : int 6
  ..- attr(*, "class")= chr "qr"
 $ df.residual  : int 66
 $ contrasts    :List of 1
  ..$ spray: chr "contr.treatment"
 $ xlevels      :List of 1
  ..$ spray: chr [1:6] "A" "B" "C" "D" ...
 $ call         : language lm(formula = count ~ spray, data = InsectSprays)
 $ terms        :Classes 'terms', 'formula' length 3 count ~ spray
  .. ..- attr(*, "variables")= language list(count, spray)
  .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:2] "count" "spray"
  .. .. .. ..$ : chr "spray"
  .. ..- attr(*, "term.labels")= chr "spray"
  .. ..- attr(*, "order")= int 1
  .. ..- attr(*, "intercept")= int 1
  .. ..- attr(*, "response")= int 1
  .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..- attr(*, "predvars")= language list(count, spray)
  .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "factor"
  .. .. ..- attr(*, "names")= chr [1:2] "count" "spray"
 $ model        :'data.frame': 72 obs. of  2 variables:
  ..$ count: num [1:72] 10 7 20 14 14 12 10 23 17 20 ...
  ..$ spray: Factor w/ 6 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 1 1 ...
  ..- attr(*, "terms")=Classes 'terms', 'formula' length 3 count ~ spray
  .. .. ..- attr(*, "variables")= language list(count, spray)
  .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:2] "count" "spray"
  .. .. .. .. ..$ : chr "spray"
  .. .. ..- attr(*, "term.labels")= chr "spray"
  .. .. ..- attr(*, "order")= int 1
  .. .. ..- attr(*, "intercept")= int 1
  .. .. ..- attr(*, "response")= int 1
  .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. .. ..- attr(*, "predvars")= language list(count, spray)
  .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "factor"
  .. .. .. ..- attr(*, "names")= chr [1:2] "count" "spray"
 - attr(*, "class")= chr "lm"


Очевидно, что отдельные элементы объекта M2 могут быть извлечены тем же путем, что и в случае с объектом M1.

Для извлечения определенных элементов из модельных объектов, созданных при помощи aov(), lm() и многих других функций в R имеются соответствующие функции-экстракторы (англ. extractor functions). Наиболее полезны, в частности:
  • coefficients() или, в сокращенном виде, coef(): извлекает коэффициенты линейной модели.
  • residuals() или, в сокращенном виде, resid():  извлекает остатки модели.
  • fitted(): извлекает предсказанные моделью значения.
Так, приведенная ниже команда позволяет рассчитать коэффициент корреляции Пирсона между исходными данными и предсказанными моделью значениями:

cor.test(fitted(M2), InsectSprays$count)
 
 Pearsons product-moment correlation
 
data:  fitted(M2) and InsectSprays$count
t = 13.5657, df = 70, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.7716214 0.9044641
sample estimates:
      cor 
0.8511398


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

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