Вступительные замечания
Мы будем упражняться на встроенном в R наборе данных «Titanic». Он описывает последствия катастрофы лайнера «Титаник» в 1912 г. и представляет из себя массив с четырьмя столбцами...
> data(Titanic)
> dimnames(Titanic)
$Class
[1] "1st" "2nd" "3rd" "Crew"
$Sex
[1] "Male" "Female"
$Age
[1] "Child" "Adult"
$Survived
[1] "No" "Yes" Как видим, здесь есть четыре категориальных переменных, а именно «Class» описывает класс каюты путешествовавших пассажиров (значения «1st», «2nd», «3rd», «Crew»), «Sex» описывает пол пассажиров (значения «Male», «Female»), «Age» описывает их возраст по стратам («Child», «Adult») и, наконец, выжил ли пассажир или погиб в этой катастрофе, указано в переменной «Survived» («No», «Yes»). В этой таблице...
> margin.table(Titanic)
[1] 2201 ...записей, каждая из которых соответствует человеку. Эти данные были собраны British Board of Trade при расследовании катастрофы.
Логлинейный анализ (Log linear analysis) позволяет взглянуть на взаимоотношения между переменными в многомерной таблице сопряжённости. Подобно своему старшему брату, критерию хи-квадрат для двумерной таблицы, логлинейный анализ основывается на допущении, что все представленные в таблице наблюдения являются независимыми (возможно, это не совсем так в случае наших данных), а ожидаемые частоты в таблицах 2 на 2 будут достаточно высоки, как правило 5 и больше. До тех пор пока ожидаемые частоты составляют 1 или больше (при этом не более 20% из них меньше пяти), точность процедуры будет приемлемой, однако её мощность будет очень мала! (собственно, это такие же допущения, как и у пирсоновского критерия хи-квадрат).
Посмотрим на данные, собранные в переменных «Class» и «Age»...
> margin.table(Titanic, c(2,4)) # эта команда показывает столбцы 2 и 4
Survived
Sex No Yes
Male 1364 367
Female 126 344 Видно, что женщина имела значительно большую вероятность выжить, чем мужчина. Отношение шансов (в советской литературе иногда переводилось как отношение преобладаний) быть выжившим в той катастрофе составляет...
> (344/126) / (367/1364)
[1] 10.14697 ... величину большую 10 (т. е. десять к одному, 10:1). Отношение правдоподобия (relative likelihood или likelihood ratio, LR) соответственно...
> (344/(344+126)) / (367/(367+1364))
[1] 3.452165 ... (отношение правдоподобия — это доля выживших женщин, делённая на долю выживших мужчин). Итак, женщины более чем в 3.5 раз чаще выживали. В этом и состоит язык логлинейного анализа. Пирсоновский критерий хи-квадрат для этой таблицы показал бы, без сомнения, высокую значимость связи (в этом контексте мы называем её «взаимодействием» (interaction) между рассмотренными двумя факторами).
Немного теории
Хи-квадрат является основным кирпичиком логлинейного анализа, однако в несколько иной форме, называемой хи-квадрат отношения правдоподобия (likelihood ratio chi square). Её преимуществом является то, что хи-квадраты отношения правдоподобия обладают свойством аддитивности. Это значит, что хи-квадраты, полученные из отдельных действий, складываются вместе и дают значения сложного эффекта (или «модели»). R (насколько мне известно) не содержит теста, позволяющего вычислить хи-квадрат отношения правдоподобия для таблицы 2 на 2, однако нам ничего не мешает пользоваться этим тестом для такой простейшей таблицы, поэтому ниже я привожу скрипт:
### начало текста скрипта
> likelihood.test = function(x) {
nrows = dim(x)[1] # no. of rows in contingency table
ncols = dim(x)[2] # no. of cols in contingency table
chi.out = chisq.test(x,correct=F) # do a Pearson chi square test
table = chi.out[[6]] # get the OFs
ratios = chi.out[[6]]/chi.out[[7]] # calculate OF/EF ratios
sum = 0 # storage for the test statistic
for (i in 1:nrows) {
for (j in 1:ncols) {
sum = sum + table[i,j]*log(ratios[i,j])
}
}
sum = 2 * sum # the likelihood ratio chi square
df = chi.out[[2]] # degrees of freedom
p = 1 - pchisq(sum,df) # p-value
out = c(sum, df, p, chi.out[[1]]) # the output vector
names(out) = c("LR-chisq","df","p-value","Pears-chisq")
round(out,4) # done!
}
### конец скрипта Выделите, скопируйте и вставьте текст в окно терминала R. При необходимости нажмите клавишу Enter в конце. Это создаст функцию в вашем рабочем пространстве с именем
likelihood.test(). Посмотрим, как она применяется: > sex.survived = margin.table(Titanic, c(2,4)) ### создаём таблицу 2 на 2
likelihood.test(sex.survived)
LR-chisq df p-value Pears-chisq
434.4688 1.0000 0.0000 456.8742 Результат выполнения команды состоит из значения отношения правдоподобия, числа степеней свободы, значения p-критерия, и, для сравнения, значения пирсоновского критерия хи-квадрат.
Попытаемся объяснить частоты в таблице сопряжённости. В традиционном критерии хи-квадрат наша модель состоит из влияния каждого фактора, при этом между факторами нет взаимодействий. Если результаты позволяют отвергнуть нулевую гипотезу, то мы считаем, что избранная модель плохо описывает частоты в клетках таблицы и что, следовательно, нужно допустить наличие взаимодействия между факторами (т. е. они не независимы или связаны).
В традиционном критерии хи-квадрат таблицы 2 на 2 более простая модель обычно не имеет смысла, но для того, чтобы показать вам саму идею, допустим, что предлагается следующая модель (нулевая гипотеза): все частоты в клетках таблицы равны. Иными словами, нет не только каких-либо взаимодействий между двумя переменными, но также нет и действия какого-либо фактора. Ничего не мешает нам протестировать нулевую гипотезу с использованием критерия хи-квадрат и мы сделаем его, приняв ожидаемые частоты равными 2201/4=550.25. На языке логлинейного анализа эта модель, в которой все частоты равны, называется нулевой (null model).
С другой стороны, мы могли бы создать модель, включающую не только главные эффекты каждого фактора, но также все возможные взаимодействия между ними. Такая модель, конечно, прекрасно бы объяснила (предсказала) все частоты в ячейках (это имело бы следствием равенство нулю значение хи-квадрат при нулевом значении степеней свободы), однако у этой модели нет никакой объясняющей ценности (мы и так знаем, что все факторы, действуя так, как они действовали, приведут к тому, что произошло). Такая модель в терминах логлинейного анализа называется «насыщенной».
Логлинейный анализ
Наша цель — объяснить полученные частоты в массиве «Titanic» возможно более простой моделью.
Эффекты, что они означают?
Class в некоторых классах пассажиров больше, чем в других Sex пассажиров одного пола больше, чем пассажиров другого Age пассажиров одной возрастной группы больше, чем пассажиров другой Survived численности выживших и погибших в катастрофе не одинаковы Class × Sex переменные «Class» и «Sex» не независимы (напр., мужчины зарабатывают больше женщин, поэтому могут позволить себе путешествовать более дорогим классом) Class × Age переменные «Class» и «Age» не независимы (напр., члены команды корабля [значение «Crew»] ехали без детей) Class × Survived переменные «Class» и «Survived» не независимы (напр., члены команды корабля способны дольше продержаться на воде после катастрофы в силу профессиональной подготовки, пассажиры верхних палуб, т. е. более дорогого класса, находились в момент катастрофы ближе к шлюпкам и спасательным кругам) Sex × Age переменные «Sex» и «Age» не независимы (напр., среди путешествовавших взрослых больше мужчин) Sex × Survived переменные «Sex» и «Survived» не независимы (напр., в силу культурных традиций в шлюпки сажали прежде всего женщин или, наоборот, мужчины, будучи сильнее, могли первыми захватить спасательные средства) Age × Survived переменные «Age» и «Survived» не независимы (напр., взрослые выживают чаще, так как они сильнее и спасательные средства рассчитаны на них) Class × Sex × Age эти три переменные связаны друг с другом Class × Sex × Survived эти три переменные связаны друг с другом Class × Age × Survived эти три переменные связаны друг с другом Sex × Age × Survived эти три переменные связаны друг с другом Class × Sex × Age × Survived все четыре переменные связаны друг с другом Заметьте, мы не рассматриваем факторы как зависимые или независимые переменные. Мы это могли бы сделать, но не в этот раз. Важно, чтобы вы об этом помнили. Вполне объяснима склонность трактовать «Survived» как зависимую переменную (отклик), на которую влияют те или иные факторы. Но мы сейчас НЕ будем этого делать. Мы попытаемся смоделировать (объяснить, предсказать, оценить) частоты в таблице сопряжённости.
Начнём, пожалуй, с такого:
> summary(Titanic)
Number of cases in table: 2201
Number of factors: 4
Test for independence of all factors:
Chisq = 1637.4, df = 25, p-value = 0
Chi-squared approximation may be incorrect Вы только что сделали логлинейный анализ модели, которая включает в себя лишь главные эффекты, т. е. эффекты каждого из факторов без каких-либо взаимодействий между ними (первые четыре строчки в подразделе «Эффекты, что они означают»). Таким образом, эту модель (нулевую гипотезу) мы отбрасываем. В некоторых случаях существуют такие взаимодействия, которые мы должны найти, чтобы адекватно объяснить частоты в ячейках.
В R можно найти несколько функций для построения логлинейной модели. Среди них
loglin() в библиотеке "stats" (она загружается по умолчанию), loglm() в библиотеке "MASS", а также glm() в библиотеке "stats". Давайте немного поэкспериментируем... > library("MASS")
> loglm( ~ Sex + Survived, data=sex.survived)
Call:
loglm(formula = ~ Sex + Survived, data = sex.survived)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 434.4688 1 0
Pearson 456.8742 1 0 Мы видим, что в таблице 2 на 2 «Sex» от «Survived», когда мы предполагаем только чистые эффекты индивидуальных факторов, мы получаем те же самые результаты, какие у нас получились при использовании моей команды
likelihood.test(). Иными словами, мы воспользовались логлинейным моделированием, чтобы применить критерий хи-квадрат на независимость переменных в этой таблице сопряжённости. Модель, содержащая лишь чистые эффекты, оказалась неприемлемой (нулевая гипотеза отброшена, p очень близко к 0). Следовательно, тут есть связь между двумя переменными (мы это предполагали). Проделаем то же самое с полной таблицей...
> loglm(~ Class + Sex + Age + Survived, data=Titanic)
Call:
loglm(formula = ~Class + Sex + Age + Survived, data = Titanic)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 1243.663 25 0
Pearson 1637.445 25 0 Результат такой же, как и при использовании команды
summary(Titanic). Существенно, что мы применили критерий хи-квадрат на независимость четырёх переменных и отбросили нулевую гипотезу. Какие-то из этих факторов взаимодействуют друг с другоми, давая наблюдаемые частоты в ячейках. Насыщенная модель, с другой стороны, не даёт ничего полезного, поскольку она полностью объясняет (предсказывает, моделирует) наблюдаемые частоты: > loglm(~ Class * Sex * Age * Survived, data=Titanic)
Call:
loglm(formula = ~Class * Sex * Age * Survived, data = Titanic)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 0 0 1
Pearson NaN 0 1 Нулевая гипотеза не отбрасывается (p очень близко к 1). Истина лежит где-то посредине между этими крайностями.
Попытаемся обойтись без четырёхпеременного взаимодействия:
> loglm(~ Class * Sex * Age * Survived - Class:Sex:Age:Survived, data=Titanic)
Call:
loglm(formula = ~Class * Sex * Age * Survived - Class:Sex:Age:Survived,
data = Titanic)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 0.0002728865 3 0.9999988
Pearson NaN 3 NaN Критерий хи-квадрат для отношения правдоподобия показывает, что эта модель предсказывает частоты в ячейках так, что они существенно не отличаются от действительно наблюдаемых. Таким образом, четырёхпеременное взаимодействие является избыточным в нашей модели, т. к. оно никак не помогает понять, каковы отношения реальных частот.
С другой стороны...
> loglm(~ Class + Sex + Age + Survived + Age:Survived, data=Titanic)
Call:
loglm(formula = ~Class + Sex + Age + Survived + Age:Survived,
data = Titanic)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 1224.103 24 0
Pearson 1596.846 24 0 ... модель со всеми факторами плюс взаимодействием Age × Survived предсказывает частоты, которые значимо отличаются от реальных. Нужно добавить больше вариантов взаимодействий, чтобы получить «истинную правду» (т. е. чтобы создать адекватную статистическую модель того, что произошло на Титанике).
Тестирование конкретных гипотез
Предположим, сначала нам следует проверить гипотезу о том, что пол связан с выживаемостью в катастрофе на Титанике. Из критерия хи-квадрат для двух переменных, который мы сделали выше, следует, что это похоже на правду. Однако, переменная «Sex» связана также с переменными «Class» (p близко к нулю) и «Age» (p близко к нулю), и обе они значимо связаны с «Survived». Следовательно, взаимодействие переменных «Class» и «Age» может испытывать влияние взаимодействия переменных «Sex» и «Survived». Подобно множественной регрессии (multiple regression) для числовых переменных логлинейная регрессия позволит нам рассмотреть отдельный вклад этих эффектов.
Если мы удалим все обозначения взаимодействий, включающие как «Sex», так и «Survived», то модель по-прежнему будет адекватно описывать наблюдаемые частоты, следовательно можно сделать вывод, что пол и выживаемость не являются связанными...
> sat.model = loglm(~ Class * Sex * Age * Survived, data=Titanic)
sat.model
Call:
loglm(formula = ~Class * Sex * Age * Survived, data = Titanic)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 0 0 1
Pearson NaN 0 1
> ### already knew this...
> model2 = update(sat.model, ~.-(Class:Sex:Age:Survived+Sex:Age:Survived+
+ Class:Sex:Survived+Sex:Survived))
> model2
Call:
loglm(formula = ~Class + Sex + Age + Survived + Class:Sex + Class:Age +
Sex:Age + Class:Survived + Age:Survived + Class:Sex:Age +
Class:Age:Survived, data = Titanic)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 436.2715 8 0
Pearson NaN 8 NaN Модель плохо описывает наблюдаемые частоты, поскольку в ней нет взаимодействия между «Sex» и «Survived». Заметьте, что использование команды
update() даёт новую модель. Конечно, мы могли бы просто перенабрать формулу модели, но команда update() является более удобным средством. Использованный нами в этой команде синтаксис означает «обновить sat.model с использованием всех тех же предикторов (точка = „те же самые“), за исключением (знак минуса) следующих за знаком минуса». Давайте сделаем её по-другому, убрав только два взаимодействия между тремя переменными... > model3 = update(sat.model, ~.-(Class:Age:Sex:Survived+Sex:Age:Survived+
+ Class:Sex:Survived))
> model3
Call:
loglm(formula = ~Class + Sex + Age + Survived + Class:Sex + Class:Age +
Sex:Age + Class:Survived + Sex:Survived + Age:Survived +
Class:Sex:Age + Class:Age:Survived, data = Titanic)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 76.90406 7 5.884182e-14
Pearson NaN 7 NaN Неа! И эту модель отбрасываем. Очевидно, взаимодействие переменных «Sex» и «Survived» обуславливается (зависит от, изменяется с) ещё одним фактором.
Как вы можете представить, такой путь поиска логлинейной модели является очень утомительным и чреватым ошибками. R предоставляет нам способ сделать это автоматически.
Функция step()
Хотя я, в целом, не одобряю автоматизированное построение статистических моделей, в данном случае я склолнен сделать исключение...
> step(sat.model, direction="backward")
Start: AIC=64
~Class * Sex * Age * Survived
Df AIC
- Class:Sex:Age:Survived 3 58
64
Step: AIC=58
~Class + Sex + Age + Survived + Class:Sex + Class:Age + Sex:Age +
Class:Survived + Sex:Survived + Age:Survived + Class:Sex:Age +
Class:Sex:Survived + Class:Age:Survived + Sex:Age:Survived
Df AIC
- Sex:Age:Survived 1 57.685
58.000
- Class:Sex:Age 3 61.783
- Class:Age:Survived 3 89.263
- Class:Sex:Survived 3 117.013
Step: AIC=57.69
~Class + Sex + Age + Survived + Class:Sex + Class:Age + Sex:Age +
Class:Survived + Sex:Survived + Age:Survived + Class:Sex:Age +
Class:Sex:Survived + Class:Age:Survived
Df AIC
57.685
- Class:Sex:Age 3 71.953
- Class:Age:Survived 3 95.899
- Class:Sex:Survived 3 126.904
Call:
loglm(formula = ~Class + Sex + Age + Survived + Class:Sex + Class:Age +
Sex:Age + Class:Survived + Sex:Survived + Age:Survived +
Class:Sex:Age + Class:Sex:Survived + Class:Age:Survived,
data = Titanic, evaluate = FALSE)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 1.685479 4 0.7933536
Pearson NaN 4 NaN Эта процедура применяется для насыщенной модели (
"sat.model"), подвергая её тому, что называется «регрессия методом исключения» (backward regression). Функция step() применяется с опцией "direction=backward" (другие опции суть "forward" и "both"). Она упрощает модель, используя критерий Akaike's Information Criterion (AIC). Чем меньше значение AIC, тем лучше. В нашем случае, R приходит к модели, в которой только два взаимодействия удалены, а именно — четырёхпеременное и Sex:Age:Survived. Это значит, что отношение между переменными «Sex» и «Survived» обуславливаются переменной «Class». Давайте это проверим... > margin.table(Titanic, c(2,4,1))
, , Class = 1st
Survived
Sex No Yes
Male 118 62
Female 4 141
, , Class = 2nd
Survived
Sex No Yes
Male 154 25
Female 13 93
, , Class = 3rd
Survived
Sex No Yes
Male 422 88
Female 106 90
, , Class = Crew
Survived
Sex No Yes
Male 670 192
Female 3 20 Отношения шансов очень легко рассчитать вручную, исходя из данных этих таблиц...
> ### odds ratio 1st class
> (141/4) / (62/118)
[1] 67.08871
> ### odds ratio 2nd class
> (93/13) / (25/154)
[1] 44.06769
> ### odds ratio 3rd class
> (90/106) / (88/422)
[1] 4.071612
> ### odds ratio crew
> (20/3) / (192/670)
[1] 23.26389 Таким образом, во всех классах шансы выжить были у женщин выше, чем у мужчин, однако это особенно заметно для совокупности пассажиров первого класса.
Получение большего объёма информации из модели
Допустим, что мы хотим работать с моделью, созданной командой
step() как с моделью случившегося на Титанике. Вначале, я помещу модель в data object (я попросту скопирую и вставлю команду, продуцирующую эту модель из вывода команды step())... > loglm(formula = ~Class + Sex + Age + Survived + Class:Sex + Class:Age +
+ Sex:Age + Class:Survived + Sex:Survived + Age:Survived +
+ Class:Sex:Age + Class:Sex:Survived + Class:Age:Survived,
+ data = Titanic) -> step.model Далее последуют несколько извлекающих функций (extractor functions) для того, чтобы получить информацию из модели (см. help, если нужен полный перечень)...
> print(step.model) # то же самое, как если бы мы просто напечатали step.model
Call:
loglm(formula = ~Class + Sex + Age + Survived + Class:Sex + Class:Age +
Sex:Age + Class:Survived + Sex:Survived + Age:Survived +
Class:Sex:Age + Class:Sex:Survived + Class:Age:Survived,
data = Titanic)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 1.685479 4 0.7933536
Pearson NaN 4 NaN
>
> fitted(step.model) # ожидаемые частоты, вычисленные по нашей модели
Re-fitting to get fitted values
, , Age = Child, Survived = No
Sex
Class Male Female
1st 0.00000 0.00000
2nd 0.00000 0.00000
3rd 37.43281 14.56719
Crew 0.00000 0.00000
, , Age = Adult, Survived = No
Sex
Class Male Female
1st 118.0000 4.0000
2nd 154.0000 13.0000
3rd 384.5672 91.4328
Crew 670.0000 3.0000
, , Age = Child, Survived = Yes
Sex
Class Male Female
1st 5.00000 1.00000
2nd 10.98493 13.01507
3rd 10.56718 16.43282
Crew 0.00000 0.00000
, , Age = Adult, Survived = Yes
Sex
Class Male Female
1st 57.00000 140.00000
2nd 14.02291 79.97709
3rd 77.43281 73.56719
Crew 192.00000 20.00000
> resid(step.model) # стандартные остатки (standardized residuals)
Re-fitting to get frequencies and fitted values
, , Age = Child, Survived = No
Sex
Class Male Female
1st 0.000000e+00 0.000000e+00
2nd 0.000000e+00 0.000000e+00
3rd -4.020602e-01 6.208006e-01
Crew 0.000000e+00 0.000000e+00
, , Age = Adult, Survived = No
Sex
Class Male Female
1st 0.000000e+00 0.000000e+00
2nd 0.000000e+00 0.000000e+00
3rd 1.239264e-01 -2.555637e-01
Crew 0.000000e+00 0.000000e+00
, , Age = Child, Survived = Yes
Sex
Class Male Female
1st 2.142148e-08 -4.552313e-08
2nd 4.546268e-03 -4.178434e-03
3rd 7.221252e-01 -6.159455e-01
Crew 0.000000e+00 0.000000e+00
, , Age = Adult, Survived = Yes
Sex
Class Male Female
1st -6.825888e-08 6.182584e-08
2nd -6.118921e-03 2.561368e-03
3rd -2.779358e-01 2.820973e-01
Crew 0.000000e+00 0.000000e+00 Излишне отмечать, что в наших ожидаемых частотах есть те, которые равны нулю. Именно поэтому, существует опасность неточности нашей модели.
Использование glm() для логлинейного моделирования
Мы также можем выполнить логлинейный анализ, воспользовавшись командой
glm(), которая используется для общих линейных моделей. Для этого нам понадобится таблица с данными (data frame), а не таблица сопряжённости, однако нет ничего сложного в том, чтобы её получить: > ti = as.data.frame(Titanic)
> ti
Class Sex Age Survived Freq
1 1st Male Child No 0
2 2nd Male Child No 0
3 3rd Male Child No 35
4 Crew Male Child No 0
5 1st Female Child No 0
6 2nd Female Child No 0
7 3rd Female Child No 17
8 Crew Female Child No 0
9 1st Male Adult No 118
10 2nd Male Adult No 154
11 3rd Male Adult No 387
12 Crew Male Adult No 670
13 1st Female Adult No 4
14 2nd Female Adult No 13
15 3rd Female Adult No 89
16 Crew Female Adult No 3
17 1st Male Child Yes 5
18 2nd Male Child Yes 11
19 3rd Male Child Yes 13
20 Crew Male Child Yes 0
21 1st Female Child Yes 1
22 2nd Female Child Yes 13
23 3rd Female Child Yes 14
24 Crew Female Child Yes 0
25 1st Male Adult Yes 57
26 2nd Male Adult Yes 14
27 3rd Male Adult Yes 75
28 Crew Male Adult Yes 192
29 1st Female Adult Yes 140
30 2nd Female Adult Yes 80
31 3rd Female Adult Yes 76
32 Crew Female Adult Yes 20 После чего логлиненый анализ выполняется так:
> glm.model = glm(Freq ~ Class * Age * Sex * Survived, data=ti, family=poisson) Это — насыщенная модель. Команда
summary(glm.model) вывалила бы нам огромный и запутанный кусок информации, но я предпочитаю пользоваться... > anova(glm.model, test="Chisq")
Analysis of Deviance Table
Model: poisson, link: log
Response: Freq
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev P(>|Chi|)
NULL 31 4953.1
Class 3 475.8 28 4477.3 8.326e-103
Age 1 2183.6 27 2293.8 0.0
Sex 1 768.3 26 1525.4 4.169e-169
Survived 1 281.8 25 1243.7 3.078e-63
Class:Age 3 148.3 22 1095.3 6.048e-32
Class:Sex 3 412.6 19 682.7 4.126e-89
Age:Sex 1 6.1 18 676.6 1.363e-02
Class:Survived 3 180.9 15 495.7 5.634e-39
Age:Survived 1 25.6 14 470.2 4.237e-07
Sex:Survived 1 353.6 13 116.6 7.053e-79
Class:Age:Sex 3 4.0 10 112.6 0.3
Class:Age:Survived 3 35.7 7 76.9 8.825e-08
Class:Sex:Survived 3 75.2 4 1.7 3.253e-16
Age:Sex:Survived 1 1.7 3 4.237e-10 0.2
Class:Age:Sex:Survived 3 0.0 0 4.463e-10 1.0 Девиация (deviance) является терминологическим вариантом критерия хи-квадрат для отношения правдоподобия. Нулевая модель, в которой все частоты одинаковы, находится вверху. Очевидно, что нулевая модель не подходит для описания этого набора данных, поскольку добавление следующих эффектов значительно уменьшает девиацию. При переходе от верхних к нижним строкам в таблице добавляется другие предикторы и это уменьшает значение девиации, указанное в столбце Deviance. Значимость p-value в последнем столбце основывается на изменениях во втором столбце (Deviance) и Df (первый столбец), происходящих при добавлении новых эффектов в модель. Здесь мы видим значимое уменьшение девиаций, свидетельствующее о том, что каждый новый эффект помогает предсказывать наблюдаемые частоты.
Эта таблица даёт возможность предположить, что целесообразно исследовать исключение трёх сочетаний переменных: четырёхпеременное взаимодействие и два трёхпеременных. Вспомним, что команда
step() сохранила взаимодействие Class:Age:Sex ... > anova(update(glm.model,.~.-(Class:Age:Sex:Survived+Age:Sex:Survived
+ +Class:Age:Sex)),test="Chisq")
Analysis of Deviance Table
Model: poisson, link: log
Response: Freq
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev P(>|Chi|)
NULL 31 4953.1
Class 3 475.8 28 4477.3 8.326e-103
Age 1 2183.6 27 2293.8 0.0
Sex 1 768.3 26 1525.4 4.169e-169
Survived 1 281.8 25 1243.7 3.078e-63
Class:Age 3 148.3 22 1095.3 6.048e-32
Class:Sex 3 412.6 19 682.7 4.126e-89
Age:Sex 1 6.1 18 676.6 1.363e-02
Class:Survived 3 180.9 15 495.7 5.634e-39
Age:Survived 1 25.6 14 470.2 4.237e-07
Sex:Survived 1 353.6 13 116.6 7.053e-79
Class:Age:Survived 3 29.2 10 87.4 2.024e-06
Class:Sex:Survived 3 65.4 7 22.0 4.066e-14
> 1-pchisq(22,df=7)
[1] 0.002540414 В этой модели последняя строка содержит остаточное отклонение (residual deviance), равное 22 при 7 степенях свободы. Судя по значению p = .0025, это неподходящая модель. Поэтому один или больше удалённых эффектов должны быть возвращены назад.
Другие извлекающие функции для модели с использованием
glm() можно посмотреть на страничке помощи, набрав ?glm. Последнее замечание
Логлинейные модели рассчитываются с использованием итеративного алгоритма. Соответственно, этот метод можно рассматривать как «машинно-нагруженный». Различные программные пакеты по разному реализуют этот алгоритм или даже используют совсем другие алгоритмы. Поэтому результаты, полученные в разных программах, могут не совпадать с абсолютой точностью. К счастью, они всё-таки будут похожи!
Немає коментарів:
Дописати коментар