Евин И. А. Введение в теорию сложных сетей [Текст] // Компьютерные исследования и моделирование. — 2010. — Т. 2, No 2. — С.
Читать далее
Заметки о жизни и об исследованиях

loglm() работает со специфической формой исходных данных — массивом (на языке R — array), отличающимся от обычной таблицы. Задача преобразовать обычную таблицу (на языке R — data.frame) в array оказалась не вполне тривиальной.
> lxx=read.csv('путь к файлу',header=T,na.strings=' ')
> names <- list(age=c('15-19','20-29','30-39','40-49','50.and.older','NA'),
family.status.legal=c('single','married','divorsed','NA'),
family.status.actual=c('alone','parent','with.man','with.woman','NA'),
clubA=c(0,1,'NA'),
partiesA=c(0,1,'NA'),
cruisingA=c(0,1,'NA'),
internet=c(0,1,'NA'),
newspaper=c(0,1,'NA'),
teletext=c(0,1,'NA'),
on.street=c(0,1,'NA'))
> lx=array(c(age=lxx$age,
family.status.legal=lxx$family.status.legal,
family.status.actual=lxx$family.status.actual,
clubA=lxx$clubA,
partiesA=lxx$partiesA,
cruisingA=lxx$cruisingA,
internet=lxx$internet,
newspaper=lxx$newspaper,
teletext=lxx$teletext,
on.street=lxx$on.street),dim=c(6,4,5,3,3,3,3,3,3,3),dimnames=names)
> data(Titanic)
> dimnames(Titanic)
$Class
[1] "1st" "2nd" "3rd" "Crew"
$Sex
[1] "Male" "Female"
$Age
[1] "Child" "Adult"
$Survived
[1] "No" "Yes" > margin.table(Titanic)
[1] 2201 > 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 > (344/(344+126)) / (367/(367+1364))
[1] 3.452165 ### начало текста скрипта
> 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!
}
### конец скрипта 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 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 все четыре переменные связаны друг с другом > 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 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 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 > 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 > 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 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 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 > 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 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 glm() можно посмотреть на страничке помощи, набрав ?glm. 
