library(tidyverse)
library(magrittr)
library(gnm)
library(broom)3 第3章
第3章は3元表に対する部分連関モデルを扱う．パッケージはtidyverse（データセットの処理のため），magrittr（パイプ演算子），DescTools（記述統計を求めるため），vcdパッケージ（カテゴリカルデータの分析のため），broom（回帰係数の整理），gnm（連関分析の処理のため），broom（tidyな回帰分析の結果の表示）を使用する．
またモデル適合度を表示するための関数をここでも準備しておく．すべてgnmによって推定を行う．
# 引数となるobjはgnmの結果
model.summary <- function(obj, Model = NULL){
if (sum(class(obj) == "gnm") != 1)
stop("estimate with gnm")
aic <- obj$deviance - obj$df * 2 # AIC(L2)
bic <- obj$deviance - obj$df * log(sum(obj$y)) #BIC(L2)
delta <- 100 * sum(abs(obj$y - obj$fitted.values)) / (2 * sum(obj$y))
p <- pchisq(obj$deviance, obj$df, lower.tail = FALSE)
result <- matrix(0, 1, 7)
if (is.null(Model)){
Model <- deparse(substitute(obj))
}
result <- tibble(
"Model Description" = Model,
"df" = obj$df,
"L2" = obj$deviance,
#"AIC(L2)" = aic,
"BIC" = bic,
"Delta" = delta,
"p" = p
)
return(result)
}3.1 表3.1A
表3.1Aをまずベクトルとして入力する．
# 表3.1
Freq <- c( 9, 5, 5, 1, 1, 6, 5, 1, 2, 2, 2, 1,
17,13, 7, 4, 13,22, 9, 1, 7,13, 6, 2,
8,14, 6, 0, 10,29,10, 0, 5,14, 6, 2,
20,38,24, 8, 23,72,34,10, 17,67,36,12,
4,21,12, 4, 7,30, 9, 1, 9,19,14, 2,
2, 9, 8, 3, 1,16,19, 2, 11,28,28,11,
0, 1, 5, 0, 2, 3, 3, 2, 2, 7, 6, 6)次に行変数，列変数，層変数を入力する．gl(n = 水準数, k = 繰り返しの数, length = 長さ)を注意しながら設定する．
polviews <- gl(n = 7, k = 4*3, length = length(Freq))
fefam <- gl(n = 4, k = 1, length = length(Freq))
natfare <- gl(n = 3, k = 4, length = length(Freq))対応しているのかを確認したければ，matrix形式で示すとよい．
matrix(polviews, nrow = 7, ncol = 12, byrow = TRUE) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
[2,] "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2" "2"
[3,] "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3" "3"
[4,] "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4" "4"
[5,] "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5" "5"
[6,] "6" "6" "6" "6" "6" "6" "6" "6" "6" "6" "6" "6"
[7,] "7" "7" "7" "7" "7" "7" "7" "7" "7" "7" "7" "7" matrix(fefam, nrow = 7, ncol = 12, byrow = TRUE) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,] "1" "2" "3" "4" "1" "2" "3" "4" "1" "2" "3" "4"
[2,] "1" "2" "3" "4" "1" "2" "3" "4" "1" "2" "3" "4"
[3,] "1" "2" "3" "4" "1" "2" "3" "4" "1" "2" "3" "4"
[4,] "1" "2" "3" "4" "1" "2" "3" "4" "1" "2" "3" "4"
[5,] "1" "2" "3" "4" "1" "2" "3" "4" "1" "2" "3" "4"
[6,] "1" "2" "3" "4" "1" "2" "3" "4" "1" "2" "3" "4"
[7,] "1" "2" "3" "4" "1" "2" "3" "4" "1" "2" "3" "4" matrix(natfare, nrow = 7, ncol = 12, byrow = TRUE) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,] "1" "1" "1" "1" "2" "2" "2" "2" "3" "3" "3" "3"
[2,] "1" "1" "1" "1" "2" "2" "2" "2" "3" "3" "3" "3"
[3,] "1" "1" "1" "1" "2" "2" "2" "2" "3" "3" "3" "3"
[4,] "1" "1" "1" "1" "2" "2" "2" "2" "3" "3" "3" "3"
[5,] "1" "1" "1" "1" "2" "2" "2" "2" "3" "3" "3" "3"
[6,] "1" "1" "1" "1" "2" "2" "2" "2" "3" "3" "3" "3"
[7,] "1" "1" "1" "1" "2" "2" "2" "2" "3" "3" "3" "3" tableのデータ（tab_3.1）とセルを単位としたデータ（freq_tab_3.1）を作成する．
tab_3.1 <- xtabs(Freq ~ polviews + fefam + natfare)
freq_tab_3.1 <- tibble(Freq, polviews, fefam, natfare)表3.2のモデルで推定を行う．モデル6から8まではスコアパラメータの正則化が必要であるが，適合度については問題なく求めることができる．
# Model 1 - 独立モデル
M1 <- freq_tab_3.1 |>
gnm(
Freq ~ polviews + fefam + natfare,
family = poisson,
data = _,
trace = F,
tolerance = 1e-12
)
# Model 2 - 完全2元交互作用モデル
M2 <- freq_tab_3.1 |>
gnm(
Freq ~ polviews * fefam + fefam * natfare + polviews * natfare,
family = poisson,
data = _,
trace = F,
tolerance = 1e-12
)
# Model 3 - 条件付き独立 (polviews)
M3 <- freq_tab_3.1 |>
gnm(
Freq ~ polviews * fefam + polviews * natfare,
family = poisson,
data = _,
trace = F,
tolerance = 1e-12
)
# Model 4 - 条件付き独立 (fefam)
M4 <- freq_tab_3.1 |>
gnm(
Freq ~ polviews * fefam + fefam * natfare,
family = poisson,
data = _,
trace = F,
tolerance = 1e-12
)
# Model 5 - 条件付き独立 (natfare)
M5 <- freq_tab_3.1 |>
gnm(
Freq ~ fefam * natfare + polviews * natfare,
family = poisson,
data = _,
trace = F,
tolerance = 1e-12
)
# Model 6 - RC(1)+RL(1) 部分連関
# 値の正則化が必要
M6 <- freq_tab_3.1 |>
gnm(
Freq ~ polviews + fefam + natfare +
Mult(1, polviews, fefam) +
Mult(1, polviews, natfare),
data = _,
family = poisson,
tolerance = 1e-12
)Initialising
Running start-up iterations..
Running main iterations.........................................................
.
Done# Model 7 - RC(1)+RL(1) 部分連関 (一貫した行スコア - polviews)
# 値の正則化が必要
M7 <- freq_tab_3.1 |>
gnm(
Freq ~ polviews + fefam + natfare +
Mult(1, polviews, fefam + natfare),
data = _,
family = poisson,
tolerance = 1e-12
)Initialising
Running start-up iterations..
Running main iterations.....................................................
Done# Model 8 - RC(1)+RL(1) 部分連関 (一貫した行スコア - polviews)
# 行スコアへの等値制約
# 値の正則化が必要
# 等値制約をかけたいセルの水準を同じにする
freq_tab_3.1 <- freq_tab_3.1 |>
mutate(polviews.c = case_when(polviews == 1 ~ 1, # mu1 = mu2
polviews == 2 ~ 1, # mu1 = mu2
polviews == 3 ~ 2,
polviews == 4 ~ 3, # mu4 = mu5
polviews == 5 ~ 3, # mu4 = mu5
polviews == 6 ~ 4,
polviews == 7 ~ 5) |>
factor())
M8 <- freq_tab_3.1 |>
gnm(
Freq ~ polviews + fefam + natfare +
Mult(1, polviews.c, fefam + natfare),
data = _,
family = poisson,
tolerance = 1e-12
)Initialising
Running start-up iterations..
Running main iterations........................................................
Doneモデル間の比較は次のように行う．
# Model 3 と Model 2の比較
pchisq(q = M3$deviance - M2$deviance,
df = M3$df.residual - M2$df.residual,
lower.tail = FALSE)[1] 0.06425806# Model 6 と Model 3の比較
pchisq(q = M6$deviance - M3$deviance,
df = M6$df.residual - M3$df.residual,
lower.tail = FALSE)[1] 0.1266967# Model 13 と Model 8の比較
## Model 13が推定できないため表の値を用いる
pchisq(q = 9.22,
df = 6,
lower.tail = FALSE)[1] 0.16157823.2 表3.2
表3.2を再現する．
## 表3.2
# リストを作成
M <- list()
# モデルの適合度をリストに順に格納
M[[1]] <- model.summary(M1, "1. Complete independence")
M[[2]] <- model.summary(M2, "2. Full two-way interaction")
M[[3]] <- model.summary(M3, "3. Conditional independence on POLVIEWS")
M[[4]] <- model.summary(M4, "4. Conditional independence on FEFAM")
M[[5]] <- model.summary(M5, "5. Conditional independence on NATFARE")
M[[6]] <- model.summary(M6, "6. RC(1) + RL(1) partial association")
M[[7]] <- model.summary(M7, "7. Model 6 plus consisitent row (POLVIEWS) score restrictions")
M[[8]] <- model.summary(M8, "8. Model 6 plus consisitent and equality restrictions on row (POLVIEWS) scores (mu_1 = mu_2, mu_4 = mu_5)")
# リストを行方向に合併する
M |> bind_rows()# A tibble: 8 × 6
`Model Description` df L2 BIC Delta p
<chr> <int> <dbl> <dbl> <dbl> <dbl>
1 1. Complete independence 72 168. -324. 14.1 1.37e-9
2 2. Full two-way interaction 36 35.4 -211. 5.72 4.99e-1
3 3. Conditional independence on POLVIEWS 42 47.3 -240. 7.30 2.67e-1
4 4. Conditional independence on FEFAM 48 87.3 -241. 10.3 4.50e-4
5 5. Conditional independence on NATFARE 54 91.0 -278. 10.1 1.21e-3
6 6. RC(1) + RL(1) partial association 57 68.6 -321. 8.83 1.40e-1
7 7. Model 6 plus consisitent row (POLVIEWS) sc… 62 72.8 -351. 9.07 1.65e-1
8 8. Model 6 plus consisitent and equality rest… 64 73.6 -364. 9.21 1.93e-13.3 表3.3
表3.3を再現するために，モデル6から8までのスコアパラメータの正則化を行う．
まずはモデル6である．
# Estimating standard error
# Model 6
# 変数と係数と係数の順番を表示
pickCoef(M6, "[.]polviews") Mult(1, ., fefam).polviews1 Mult(1, ., fefam).polviews2
14 15
Mult(1, ., fefam).polviews3 Mult(1, ., fefam).polviews4
16 17
Mult(1, ., fefam).polviews5 Mult(1, ., fefam).polviews6
18 19
Mult(1, ., fefam).polviews7 Mult(1, ., natfare).polviews1
20 26
Mult(1, ., natfare).polviews2 Mult(1, ., natfare).polviews3
27 28
Mult(1, ., natfare).polviews4 Mult(1, ., natfare).polviews5
29 30
Mult(1, ., natfare).polviews6 Mult(1, ., natfare).polviews7
31 32 pickCoef(M6, "[.]fefam")Mult(1, polviews, .).fefam1 Mult(1, polviews, .).fefam2
21 22
Mult(1, polviews, .).fefam3 Mult(1, polviews, .).fefam4
23 24 pickCoef(M6, "[.]natfare")Mult(1, polviews, .).natfare1 Mult(1, polviews, .).natfare2
33 34
Mult(1, polviews, .).natfare3
35 data.frame(var = names(M6$coefficients),
estimate = M6$coefficients) |>
mutate(estimate = estimate,
number = row_number()) var estimate number
(Intercept) (Intercept) 0.43666257 1
polviews2 polviews2 1.12712181 2
polviews3 polviews3 1.08818509 3
polviews4 polviews4 2.40474437 4
polviews5 polviews5 1.38673165 5
polviews6 polviews6 1.33444426 6
polviews7 polviews7 -0.05728770 7
fefam2 fefam2 0.98926014 8
fefam3 fefam3 0.45096067 9
fefam4 fefam4 -0.84928436 10
natfare2 natfare2 0.34302758 11
natfare3 natfare3 0.26650624 12
Mult(., polviews, fefam). Mult(., polviews, fefam). 0.63662574 13
Mult(1, ., fefam).polviews1 Mult(1, ., fefam).polviews1 -0.52430127 14
Mult(1, ., fefam).polviews2 Mult(1, ., fefam).polviews2 -1.37269517 15
Mult(1, ., fefam).polviews3 Mult(1, ., fefam).polviews3 -1.08358988 16
Mult(1, ., fefam).polviews4 Mult(1, ., fefam).polviews4 0.03445392 17
Mult(1, ., fefam).polviews5 Mult(1, ., fefam).polviews5 -0.12334864 18
Mult(1, ., fefam).polviews6 Mult(1, ., fefam).polviews6 1.41253860 19
Mult(1, ., fefam).polviews7 Mult(1, ., fefam).polviews7 2.18532146 20
Mult(1, polviews, .).fefam1 Mult(1, polviews, .).fefam1 -0.60991027 21
Mult(1, polviews, .).fefam2 Mult(1, polviews, .).fefam2 -0.12478151 22
Mult(1, polviews, .).fefam3 Mult(1, polviews, .).fefam3 0.37414329 23
Mult(1, polviews, .).fefam4 Mult(1, polviews, .).fefam4 0.67829418 24
Mult(., polviews, natfare). Mult(., polviews, natfare). 0.08906570 25
Mult(1, ., natfare).polviews1 Mult(1, ., natfare).polviews1 -16.37705391 26
Mult(1, ., natfare).polviews2 Mult(1, ., natfare).polviews2 -7.73309186 27
Mult(1, ., natfare).polviews3 Mult(1, ., natfare).polviews3 -4.67504815 28
Mult(1, ., natfare).polviews4 Mult(1, ., natfare).polviews4 1.14976253 29
Mult(1, ., natfare).polviews5 Mult(1, ., natfare).polviews5 -1.85945949 30
Mult(1, ., natfare).polviews6 Mult(1, ., natfare).polviews6 12.83349452 31
Mult(1, ., natfare).polviews7 Mult(1, ., natfare).polviews7 12.88783913 32
Mult(1, polviews, .).natfare1 Mult(1, polviews, .).natfare1 -0.47314801 33
Mult(1, polviews, .).natfare2 Mult(1, polviews, .).natfare2 -0.17809970 34
Mult(1, polviews, .).natfare3 Mult(1, polviews, .).natfare3 0.45798464 35各スコアの両端のカテゴリ，つまり上記推定値のリストの 14,20,21,24,26,32,33,35番目の8つの値についてスコアを求め，固定する．
# R POLVIEWS mu1[i], i = 1 to 7
mu1 <- getContrasts(model = M6,
set = pickCoef(M6, "fefam)[.]polviews"),
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit")
# C FEFAM nu[j], j = 1 to 4
nu <- getContrasts(model = M6,
set = pickCoef(M6, "[.]fefam"),
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit")
# R POLVIEWS mu2[i], i = 1 to 7
mu2 <- getContrasts(model = M6,
set = pickCoef(M6, "natfare)[.]polviews"),
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit")
# L NATFARE eta[k], k = 1 to 3
eta <- getContrasts(model = M6,
set = pickCoef(M6, "[.]natfare"),
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit")
# 確認
list(mu1 = mu1, nu = nu, mu2 = mu2, eta = eta)$mu1
Estimate Std. Error
Mult(1, ., fefam).polviews1 -0.18890704 0.17539812
Mult(1, ., fefam).polviews2 -0.45611586 0.11238050
Mult(1, ., fefam).polviews3 -0.36505971 0.11498236
Mult(1, ., fefam).polviews4 -0.01292237 0.08199631
Mult(1, ., fefam).polviews5 -0.06262362 0.11058185
Mult(1, ., fefam).polviews6 0.42111705 0.12367540
Mult(1, ., fefam).polviews7 0.66451154 0.10976226
$nu
Estimate Std. Error
Mult(1, polviews, .).fefam1 -0.7026886 0.06444123
Mult(1, polviews, .).fefam2 -0.2081705 0.09371704
Mult(1, polviews, .).fefam3 0.3004108 0.10956493
Mult(1, polviews, .).fefam4 0.6104483 0.09158453
$mu2
Estimate Std. Error
Mult(1, ., natfare).polviews1 -0.60584367 0.12394769
Mult(1, ., natfare).polviews2 -0.27518966 0.11469790
Mult(1, ., natfare).polviews3 -0.15821154 0.11316687
Mult(1, ., natfare).polviews4 0.06460260 0.07187236
Mult(1, ., natfare).polviews5 -0.05050796 0.10030828
Mult(1, ., natfare).polviews6 0.51153570 0.10405220
Mult(1, ., natfare).polviews7 0.51361453 0.12512860
$eta
Estimate Std. Error
Mult(1, polviews, .).natfare1 -0.6073482 0.07616839
Mult(1, polviews, .).natfare2 -0.1689209 0.11149916
Mult(1, polviews, .).natfare3 0.7762692 0.03533078正則化したスコアを求めることができたので， 各スコアの両端のカテゴリを取り出し，conとする．
con <- c(mu1$qvframe[,1][c(1,7)],
nu$qvframe[,1][c(1,4)],
mu2$qvframe[,1][c(1,7)],
eta$qvframe[,1][c(1,3)])conを用いて制約を課す．制約を課す対象となるパラメータはconstrainで指定する．
M6_SE <- freq_tab_3.1 |>
gnm(Freq ~ polviews + fefam + natfare +
Mult(1, polviews, fefam) + Mult(1, polviews, natfare),
constrain = c(14,20,21,24,26,32,33,35),
constrainTo = con,
data = _,
family = poisson,
tolerance = 1e-12)Initialising
Running start-up iterations..
Running main iterations.........................................................
.............
DoneM6_SE_coef <- tidy(M6_SE)Warning: The `tidy()` method for objects of class `gnm` is not maintained by the broom team, and is only supported through the `glm` tidier method. Please be cautious in interpreting and reporting broom output.
This warning is displayed once per session.M6_SE_coef |> print(n = Inf)# A tibble: 35 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 0.491 0.271 1.81 7.07e- 2
2 polviews2 1.03 0.229 4.52 6.06e- 6
3 polviews3 0.993 0.226 4.39 1.12e- 5
4 polviews4 2.33 0.212 11.0 3.35e-28
5 polviews5 1.32 0.218 6.08 1.24e- 9
6 polviews6 1.26 0.282 4.49 7.13e- 6
7 polviews7 -0.0882 0.355 -0.248 8.04e- 1
8 fefam2 1.01 0.188 5.40 6.70e- 8
9 fefam3 0.498 0.340 1.46 1.43e- 1
10 fefam4 -0.787 0.449 -1.75 7.96e- 2
11 natfare2 0.329 0.127 2.60 9.35e- 3
12 natfare3 0.222 0.314 0.706 4.80e- 1
13 Mult(., polviews, fefam). 1.98 0.744 2.66 7.73e- 3
14 Mult(1, ., fefam).polviews1 -0.189 NA NA NA
15 Mult(1, ., fefam).polviews2 -0.456 0.308 -1.48 1.39e- 1
16 Mult(1, ., fefam).polviews3 -0.365 0.285 -1.28 2.00e- 1
17 Mult(1, ., fefam).polviews4 -0.0129 0.184 -0.0703 9.44e- 1
18 Mult(1, ., fefam).polviews5 -0.0626 0.212 -0.296 7.67e- 1
19 Mult(1, ., fefam).polviews6 0.421 0.205 2.06 3.98e- 2
20 Mult(1, ., fefam).polviews7 0.665 NA NA NA
21 Mult(1, polviews, .).fefam1 -0.703 NA NA NA
22 Mult(1, polviews, .).fefam2 -0.208 0.152 -1.37 1.71e- 1
23 Mult(1, polviews, .).fefam3 0.300 0.185 1.62 1.05e- 1
24 Mult(1, polviews, .).fefam4 0.610 NA NA NA
25 Mult(., polviews, natfare). 1.57 0.405 3.87 1.10e- 4
26 Mult(1, ., natfare).polviews1 -0.606 NA NA NA
27 Mult(1, ., natfare).polviews2 -0.275 0.200 -1.38 1.68e- 1
28 Mult(1, ., natfare).polviews3 -0.158 0.191 -0.828 4.08e- 1
29 Mult(1, ., natfare).polviews4 0.0646 0.156 0.414 6.79e- 1
30 Mult(1, ., natfare).polviews5 -0.0505 0.177 -0.286 7.75e- 1
31 Mult(1, ., natfare).polviews6 0.512 0.211 2.43 1.53e- 2
32 Mult(1, ., natfare).polviews7 0.514 NA NA NA
33 Mult(1, polviews, .).natfare1 -0.607 NA NA NA
34 Mult(1, polviews, .).natfare2 -0.169 0.175 -0.967 3.34e- 1
35 Mult(1, polviews, .).natfare3 0.776 NA NA NA モデル7は少し難しい．モデル7は一貫した制約を課している．このモデルからまず一貫したpolviewsについて正則化した場合の値を求める．これが求まったら，次は一貫した制約を課していないモデルに，先ほどの正則化したスコアをRCとRLのそれぞれの部分に適用する．
# Model 7
# 変数と係数と係数の順番を表示
pickCoef(M7, "[.]polviews")Mult(1, ., fefam + natfare).polviews1 Mult(1, ., fefam + natfare).polviews2
14 15
Mult(1, ., fefam + natfare).polviews3 Mult(1, ., fefam + natfare).polviews4
16 17
Mult(1, ., fefam + natfare).polviews5 Mult(1, ., fefam + natfare).polviews6
18 19
Mult(1, ., fefam + natfare).polviews7
20 pickCoef(M7, "[.]fefam")Mult(1, polviews, . + natfare).fefam1 Mult(1, polviews, . + natfare).fefam2
21 22
Mult(1, polviews, . + natfare).fefam3 Mult(1, polviews, . + natfare).fefam4
23 24 pickCoef(M7, "[.]natfare")Mult(1, polviews, fefam + .).natfare2 Mult(1, polviews, fefam + .).natfare3
25 26 data.frame(var = names(M7$coefficients),
estimate = M7$coefficients) |>
mutate(estimate = estimate,
number = row_number()) var
(Intercept) (Intercept)
polviews2 polviews2
polviews3 polviews3
polviews4 polviews4
polviews5 polviews5
polviews6 polviews6
polviews7 polviews7
fefam2 fefam2
fefam3 fefam3
fefam4 fefam4
natfare2 natfare2
natfare3 natfare3
Mult(., polviews, fefam + natfare). Mult(., polviews, fefam + natfare).
Mult(1, ., fefam + natfare).polviews1 Mult(1, ., fefam + natfare).polviews1
Mult(1, ., fefam + natfare).polviews2 Mult(1, ., fefam + natfare).polviews2
Mult(1, ., fefam + natfare).polviews3 Mult(1, ., fefam + natfare).polviews3
Mult(1, ., fefam + natfare).polviews4 Mult(1, ., fefam + natfare).polviews4
Mult(1, ., fefam + natfare).polviews5 Mult(1, ., fefam + natfare).polviews5
Mult(1, ., fefam + natfare).polviews6 Mult(1, ., fefam + natfare).polviews6
Mult(1, ., fefam + natfare).polviews7 Mult(1, ., fefam + natfare).polviews7
Mult(1, polviews, . + natfare).fefam1 Mult(1, polviews, . + natfare).fefam1
Mult(1, polviews, . + natfare).fefam2 Mult(1, polviews, . + natfare).fefam2
Mult(1, polviews, . + natfare).fefam3 Mult(1, polviews, . + natfare).fefam3
Mult(1, polviews, . + natfare).fefam4 Mult(1, polviews, . + natfare).fefam4
Mult(1, polviews, fefam + .).natfare2 Mult(1, polviews, fefam + .).natfare2
Mult(1, polviews, fefam + .).natfare3 Mult(1, polviews, fefam + .).natfare3
estimate number
(Intercept) 0.670513524 1
polviews2 1.048946988 2
polviews3 0.954549677 3
polviews4 2.119241242 4
polviews5 1.146390461 5
polviews6 0.810062282 6
polviews7 -0.622997763 7
fefam2 1.005922049 8
fefam3 0.476830092 9
fefam4 -0.810480292 10
natfare2 0.355772837 11
natfare3 0.330454996 12
Mult(., polviews, fefam + natfare). -2.568730061 13
Mult(1, ., fefam + natfare).polviews1 -0.528532311 14
Mult(1, ., fefam + natfare).polviews2 -0.495517192 15
Mult(1, ., fefam + natfare).polviews3 -0.353118498 16
Mult(1, ., fefam + natfare).polviews4 0.008241641 17
Mult(1, ., fefam + natfare).polviews5 -0.092973326 18
Mult(1, ., fefam + natfare).polviews6 0.577520277 19
Mult(1, ., fefam + natfare).polviews7 0.707468904 20
Mult(1, polviews, . + natfare).fefam1 0.555791776 21
Mult(1, polviews, . + natfare).fefam2 0.209768260 22
Mult(1, polviews, . + natfare).fefam3 -0.085152903 23
Mult(1, polviews, . + natfare).fefam4 -0.265028472 24
Mult(1, polviews, fefam + .).natfare2 -0.171067190 25
Mult(1, polviews, fefam + .).natfare3 -0.627642391 26# mu[i], i = 1 to 7
mu <- getContrasts(model = M7,
set = pickCoef(M7, "[.]polviews")[1:7],
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit")
con <- c(mu$qvframe[,1][c(1:7)],
mu$qvframe[,1][c(1:7)])ここでは一貫した制約を課さないモデルを用い，制約で一貫したスコアを割り当てる．そして，制約のない列スコアと層スコアを正則化する．
M7_SE <- freq_tab_3.1 |>
gnm(Freq ~ polviews + fefam + natfare +
Mult(1, polviews, fefam) +
Mult(1, polviews, natfare),
constrain = c(14:20, 26:32),
constrainTo = con,
data = _,
family = poisson,
tolerance = 1e-12)Initialising
Running start-up iterations..
Running main iterations...
Done# nu[j], j = 1 to 4
nu <- getContrasts(model = M7_SE,
set = pickCoef(M7_SE, "[.]fefam"),
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit")
# eta[k], k = 1 to 3
eta <- getContrasts(model = M7_SE,
set = pickCoef(M7_SE, "[.]natfare"),
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit")
con <- c(mu$qvframe[,1][c(1:7)],
nu$qvframe[,1][c(1,4)],
mu$qvframe[,1][c(1:7)],
eta$qvframe[,1][c(1,3)])最後に一貫したスコア（行）と正則化したスコア（列と層）を用いて再度推定を行う．
M7_SE <- freq_tab_3.1 |>
gnm(Freq ~ polviews + fefam + natfare +
Mult(1, polviews, fefam) +
Mult(1, polviews, natfare),
constrain = c(14:20,21,24,26:32,33,35),
constrainTo = con,
data = _,
family = poisson,
trace = T,
tolerance = 1e-12)Initialising
Initial Deviance = 177.313021
Running start-up iterations
Start-up iteration 1. Deviance = 76.224232
Start-up iteration 2. Deviance = 73.047675
Running main iterations
Iteration 1. Deviance = 72.774263
Iteration 2. Deviance = 72.773834
Iteration 3. Deviance = 72.773834
Iteration 4. Deviance = 72.773834
Iteration 5. Deviance = 72.773834
Iteration 6. Deviance = 72.773834
Iteration 7. Deviance = 72.773834
Iteration 8. Deviance = 72.773834
Iteration 9. Deviance = 72.773834
Iteration 10. Deviance = 72.773834
Iteration 11. Deviance = 72.773834
Iteration 12. Deviance = 72.773834
Iteration 13. Deviance = 72.773834
Iteration 14. Deviance = 72.773834
Iteration 15. Deviance = 72.773834
Iteration 16. Deviance = 72.773834
DoneM7_SE_coef <- tidy(M7_SE)
M7_SE_coef |> print(n = Inf)# A tibble: 35 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 0.497 0.200 2.49 1.28e- 2
2 polviews2 1.06 0.184 5.78 7.52e- 9
3 polviews3 1.03 0.190 5.42 5.91e- 8
4 polviews4 2.34 0.200 11.7 1.01e-31
5 polviews5 1.33 0.202 6.59 4.46e-11
6 polviews6 1.27 0.287 4.42 9.73e- 6
7 polviews7 -0.107 0.338 -0.317 7.51e- 1
8 fefam2 0.983 0.0947 10.4 3.01e-25
9 fefam3 0.435 0.104 4.19 2.73e- 5
10 fefam4 -0.864 0.150 -5.77 7.72e- 9
11 natfare2 0.345 0.0853 4.04 5.28e- 5
12 natfare3 0.290 0.0869 3.33 8.54e- 4
13 Mult(., polviews, fefam). 1.95 0.374 5.21 1.92e- 7
14 Mult(1, ., fefam).polviews1 -0.413 NA NA NA
15 Mult(1, ., fefam).polviews2 -0.386 NA NA NA
16 Mult(1, ., fefam).polviews3 -0.269 NA NA NA
17 Mult(1, ., fefam).polviews4 0.0275 NA NA NA
18 Mult(1, ., fefam).polviews5 -0.0555 NA NA NA
19 Mult(1, ., fefam).polviews6 0.494 NA NA NA
20 Mult(1, ., fefam).polviews7 0.601 NA NA NA
21 Mult(1, polviews, .).fefam1 -0.726 NA NA NA
22 Mult(1, polviews, .).fefam2 -0.170 0.150 -1.14 2.55e- 1
23 Mult(1, polviews, .).fefam3 0.304 0.189 1.61 1.08e- 1
24 Mult(1, polviews, .).fefam4 0.593 NA NA NA
25 Mult(., polviews, natfare). 1.44 0.226 6.36 2.03e-10
26 Mult(1, ., natfare).polviews1 -0.413 NA NA NA
27 Mult(1, ., natfare).polviews2 -0.386 NA NA NA
28 Mult(1, ., natfare).polviews3 -0.269 NA NA NA
29 Mult(1, ., natfare).polviews4 0.0275 NA NA NA
30 Mult(1, ., natfare).polviews5 -0.0555 NA NA NA
31 Mult(1, ., natfare).polviews6 0.494 NA NA NA
32 Mult(1, ., natfare).polviews7 0.601 NA NA NA
33 Mult(1, polviews, .).natfare1 -0.580 NA NA NA
34 Mult(1, polviews, .).natfare2 -0.207 0.186 -1.12 2.64e- 1
35 Mult(1, polviews, .).natfare3 0.788 NA NA NA モデル8は1,1,2,3,3,4,5という制約があり，結局は1,2,3,4,5の5つのパラメータを推定することになる．しかし，あくまで1,1,2,3,3,4,5という7つのパラメータにウェイトをかける必要がある．そこでc(2,1,2,1,1)/7とすることで，等値制約をかけている部分は平均を2倍にし，c(2,1,2,1,1)とすることでウェイトも2倍にしている．
# Model 8
# mu[i], i = 1 to 7
mu <- getContrasts(model = M8,
set = pickCoef(M8, "[.]polviews.c"),
ref = c(2, 1, 2, 1, 1) / 7,
scaleRef = c(2, 1, 2, 1, 1) / 7,
scaleWeights = c(2, 1, 2, 1, 1))
mu <- getContrasts(model = M8,
set = pickCoef(M8, "[.]polviews.c"),
ref = c(2, 1, 2, 1, 1) / 7,
scaleRef = c(2, 1, 2, 1, 1) / 7,
scaleWeights = c(2, 1, 2, 1, 1))
mu;mu Estimate Std. Error
Mult(1, ., fefam + natfare).polviews.c1 0.402725149 0.04298356
Mult(1, ., fefam + natfare).polviews.c2 0.278981096 0.08172199
Mult(1, ., fefam + natfare).polviews.c3 0.001961895 0.04688021
Mult(1, ., fefam + natfare).polviews.c4 -0.491602614 0.08428340
Mult(1, ., fefam + natfare).polviews.c5 -0.596752569 0.08715205 Estimate Std. Error
Mult(1, ., fefam + natfare).polviews.c1 0.402725149 0.04298356
Mult(1, ., fefam + natfare).polviews.c2 0.278981096 0.08172199
Mult(1, ., fefam + natfare).polviews.c3 0.001961895 0.04688021
Mult(1, ., fefam + natfare).polviews.c4 -0.491602614 0.08428340
Mult(1, ., fefam + natfare).polviews.c5 -0.596752569 0.08715205con <- c(mu$qvframe[, 1][c(1, 1, 2, 3, 3, 4, 5)],
mu$qvframe[, 1][c(1, 1, 2, 3, 3, 4, 5)])一貫した行スコア（等値制約あり）を指定して再推定．
M8 <- freq_tab_3.1 |>
gnm(Freq ~ polviews + fefam + natfare +
Mult(1, polviews, fefam) +
Mult(1, polviews, natfare),
constrain = c(14:20,26:32),
constrainTo = con,
data = _,
family = poisson,
tolerance = 1e-12)Initialising
Running start-up iterations..
Running main iterations...
Donetidy(M8)# A tibble: 35 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 0.477 NA NA NA
2 polviews2 1.05 0.184 5.70 1.20e- 8
3 polviews3 1.02 NA NA NA
4 polviews4 2.34 NA NA NA
5 polviews5 1.34 NA NA NA
6 polviews6 1.29 NA NA NA
7 polviews7 -0.0810 NA NA NA
8 fefam2 0.993 0.0956 10.4 2.62e-25
9 fefam3 0.451 0.104 4.33 1.52e- 5
10 fefam4 -0.840 0.149 -5.64 1.72e- 8
# ℹ 25 more rowsglance(M8)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 NA NA -189. 412. 453. 73.6 67 84次は列スコアと層スコアの正則化をおこなう．
# nu[j], j = 1 to 4
nu <- getContrasts(model = M8,
set = pickCoef(M8, "[.]fefam"),
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit")
# eta[k], k = 1 to 3
eta <- getContrasts(model = M8,
set = pickCoef(M8, "[.]natfare"),
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit")
con <- c(mu$qvframe[,1][c(1,1,2,3,3,4,5)],
nu$qvframe[,1][c(1,4)],
mu$qvframe[,1][c(1,1,2,3,3,4,5)],
eta$qvframe[,1][c(1,3)])一貫した行スコア（等値制約あり）のすべて，列スコア（最初と最後），層スコア（最初と最後）を指定して再推定．
M8_SE <- freq_tab_3.1 |>
gnm(Freq ~ polviews + fefam + natfare +
Mult(1, polviews, fefam) +
Mult(1, polviews, natfare),
constrain = c(14:20,21,24,26:32,33,35),
constrainTo = con,
data = _,
family = poisson,
tolerance = 1e-12)Initialising
Running start-up iterations..
Running main iterations...
DoneM8_SE_coef <- tidy(M8_SE)
M8_SE_coef |> print(n = Inf)# A tibble: 35 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 0.494 0.199 2.49 1.29e- 2
2 polviews2 1.05 0.184 5.70 1.20e- 8
3 polviews3 1.02 0.189 5.38 7.41e- 8
4 polviews4 2.33 0.195 12.0 6.32e-33
5 polviews5 1.32 0.207 6.39 1.62e-10
6 polviews6 1.25 0.286 4.39 1.12e- 5
7 polviews7 -0.123 0.337 -0.366 7.14e- 1
8 fefam2 0.993 0.0956 10.4 2.62e-25
9 fefam3 0.451 0.104 4.33 1.52e- 5
10 fefam4 -0.840 0.149 -5.64 1.72e- 8
11 natfare2 0.346 0.0856 4.04 5.36e- 5
12 natfare3 0.302 0.0868 3.48 5.01e- 4
13 Mult(., polviews, fefam). 1.94 0.374 5.19 2.16e- 7
14 Mult(1, ., fefam).polviews1 0.403 NA NA NA
15 Mult(1, ., fefam).polviews2 0.403 NA NA NA
16 Mult(1, ., fefam).polviews3 0.279 NA NA NA
17 Mult(1, ., fefam).polviews4 0.00196 NA NA NA
18 Mult(1, ., fefam).polviews5 0.00196 NA NA NA
19 Mult(1, ., fefam).polviews6 -0.492 NA NA NA
20 Mult(1, ., fefam).polviews7 -0.597 NA NA NA
21 Mult(1, polviews, .).fefam1 0.733 NA NA NA
22 Mult(1, polviews, .).fefam2 0.162 0.150 1.08 2.78e- 1
23 Mult(1, polviews, .).fefam3 -0.313 0.191 -1.64 1.01e- 1
24 Mult(1, polviews, .).fefam4 -0.582 NA NA NA
25 Mult(., polviews, natfare). -1.41 0.226 -6.25 4.03e-10
26 Mult(1, ., natfare).polviews1 0.403 NA NA NA
27 Mult(1, ., natfare).polviews2 0.403 NA NA NA
28 Mult(1, ., natfare).polviews3 0.279 NA NA NA
29 Mult(1, ., natfare).polviews4 0.00196 NA NA NA
30 Mult(1, ., natfare).polviews5 0.00196 NA NA NA
31 Mult(1, ., natfare).polviews6 -0.492 NA NA NA
32 Mult(1, ., natfare).polviews7 -0.597 NA NA NA
33 Mult(1, polviews, .).natfare1 -0.569 NA NA NA
34 Mult(1, polviews, .).natfare2 -0.222 0.190 -1.17 2.42e- 1
35 Mult(1, polviews, .).natfare3 0.791 NA NA NA すべてのスコアを列方向に合併し，表3.3を再現．
bind_cols(
M6_SE_coef |> dplyr::select(term, Model6 = estimate),
M7_SE_coef |> dplyr::select(Model7 = estimate),
M8_SE_coef |> dplyr::select(Model8 = estimate)) |>
dplyr::filter(grepl("Mult", term)) |>
print(n = Inf)# A tibble: 23 × 4
term Model6 Model7 Model8
<chr> <dbl> <dbl> <dbl>
1 Mult(., polviews, fefam). 1.98 1.95 1.94
2 Mult(1, ., fefam).polviews1 -0.189 -0.413 0.403
3 Mult(1, ., fefam).polviews2 -0.456 -0.386 0.403
4 Mult(1, ., fefam).polviews3 -0.365 -0.269 0.279
5 Mult(1, ., fefam).polviews4 -0.0129 0.0275 0.00196
6 Mult(1, ., fefam).polviews5 -0.0626 -0.0555 0.00196
7 Mult(1, ., fefam).polviews6 0.421 0.494 -0.492
8 Mult(1, ., fefam).polviews7 0.665 0.601 -0.597
9 Mult(1, polviews, .).fefam1 -0.703 -0.726 0.733
10 Mult(1, polviews, .).fefam2 -0.208 -0.170 0.162
11 Mult(1, polviews, .).fefam3 0.300 0.304 -0.313
12 Mult(1, polviews, .).fefam4 0.610 0.593 -0.582
13 Mult(., polviews, natfare). 1.57 1.44 -1.41
14 Mult(1, ., natfare).polviews1 -0.606 -0.413 0.403
15 Mult(1, ., natfare).polviews2 -0.275 -0.386 0.403
16 Mult(1, ., natfare).polviews3 -0.158 -0.269 0.279
17 Mult(1, ., natfare).polviews4 0.0646 0.0275 0.00196
18 Mult(1, ., natfare).polviews5 -0.0505 -0.0555 0.00196
19 Mult(1, ., natfare).polviews6 0.512 0.494 -0.492
20 Mult(1, ., natfare).polviews7 0.514 0.601 -0.597
21 Mult(1, polviews, .).natfare1 -0.607 -0.580 -0.569
22 Mult(1, polviews, .).natfare2 -0.169 -0.207 -0.222
23 Mult(1, polviews, .).natfare3 0.776 0.788 0.791 3.4 表3.4
表3.4のデータを入力する．
# 表3.1
Freq <- c(76, 14, 15, 4,
32, 17, 7, 3,
64, 23, 28, 15,
41, 11, 27, 16,
15, 2, 7, 4,
27, 20, 9, 5,
57, 31, 24, 15,
27, 9, 22, 16,
13, 6, 13, 5,
12, 13, 10, 6,
46, 32, 75, 20,
54, 26, 58, 55,
7, 6, 7, 6,
7, 2, 3, 6,
12, 11, 31, 15,
52, 36, 80,101)行，列，層のカテゴリを入力し，matrixで確認する．
L <- gl(4, 16, 64)
R <- gl(4, 4, 64)
C <- gl(4, 1, 64)
# 確認する
L |> matrix(nrow = 16, ncol = 4, byrow = TRUE) [,1] [,2] [,3] [,4]
[1,] "1" "1" "1" "1"
[2,] "1" "1" "1" "1"
[3,] "1" "1" "1" "1"
[4,] "1" "1" "1" "1"
[5,] "2" "2" "2" "2"
[6,] "2" "2" "2" "2"
[7,] "2" "2" "2" "2"
[8,] "2" "2" "2" "2"
[9,] "3" "3" "3" "3"
[10,] "3" "3" "3" "3"
[11,] "3" "3" "3" "3"
[12,] "3" "3" "3" "3"
[13,] "4" "4" "4" "4"
[14,] "4" "4" "4" "4"
[15,] "4" "4" "4" "4"
[16,] "4" "4" "4" "4" R |> matrix(nrow = 16, ncol = 4, byrow = TRUE) [,1] [,2] [,3] [,4]
[1,] "1" "1" "1" "1"
[2,] "2" "2" "2" "2"
[3,] "3" "3" "3" "3"
[4,] "4" "4" "4" "4"
[5,] "1" "1" "1" "1"
[6,] "2" "2" "2" "2"
[7,] "3" "3" "3" "3"
[8,] "4" "4" "4" "4"
[9,] "1" "1" "1" "1"
[10,] "2" "2" "2" "2"
[11,] "3" "3" "3" "3"
[12,] "4" "4" "4" "4"
[13,] "1" "1" "1" "1"
[14,] "2" "2" "2" "2"
[15,] "3" "3" "3" "3"
[16,] "4" "4" "4" "4" C |> matrix(nrow = 16, ncol = 4, byrow = TRUE) [,1] [,2] [,3] [,4]
[1,] "1" "2" "3" "4"
[2,] "1" "2" "3" "4"
[3,] "1" "2" "3" "4"
[4,] "1" "2" "3" "4"
[5,] "1" "2" "3" "4"
[6,] "1" "2" "3" "4"
[7,] "1" "2" "3" "4"
[8,] "1" "2" "3" "4"
[9,] "1" "2" "3" "4"
[10,] "1" "2" "3" "4"
[11,] "1" "2" "3" "4"
[12,] "1" "2" "3" "4"
[13,] "1" "2" "3" "4"
[14,] "1" "2" "3" "4"
[15,] "1" "2" "3" "4"
[16,] "1" "2" "3" "4" 以上をtibbleでまとめ，データを作成する．
freq_tab_3.1 <- tibble(Freq, L, R, C) |> arrange(L, R, C)
freq_tab_3.1# A tibble: 64 × 4
Freq L R C
<dbl> <fct> <fct> <fct>
1 76 1 1 1
2 14 1 1 2
3 15 1 1 3
4 4 1 1 4
5 32 1 2 1
6 17 1 2 2
7 7 1 2 3
8 3 1 2 4
9 64 1 3 1
10 23 1 3 2
# ℹ 54 more rows行，列，層について整数スコアを作成する．
freq_tab_3.1 <- freq_tab_3.1 |>
mutate(Rscore = as.numeric(R),
Cscore = as.numeric(C),
Lscore = as.numeric(L))
freq_tab_3.1# A tibble: 64 × 7
Freq L R C Rscore Cscore Lscore
<dbl> <fct> <fct> <fct> <dbl> <dbl> <dbl>
1 76 1 1 1 1 1 1
2 14 1 1 2 1 2 1
3 15 1 1 3 1 3 1
4 4 1 1 4 1 4 1
5 32 1 2 1 2 1 1
6 17 1 2 2 2 2 1
7 7 1 2 3 2 3 1
8 3 1 2 4 2 4 1
9 64 1 3 1 3 1 1
10 23 1 3 2 3 2 1
# ℹ 54 more rowsそれではモデル1から15までを用いた推定を行う．ただしgnmではモデル3と5の推定はできない． また
# Model 1 - Complete Independence
model1 <- freq_tab_3.1 |>
gnm(
Freq ~ R + C + L,
data = _,
family = poisson,
trace = F,
tolerance = 1e-12
)
# Model 2 - Unrestricted RC(1)+RL(1)+CL(1)
model2 <- freq_tab_3.1 |>
gnm(
Freq ~ R + C + L + Mult(1, R, C) + Mult(1, R, L) + Mult(1, C, L),
data = _,
family = poisson,
trace = F,
tolerance = 1e-12
)Initialising
Running start-up iterations..
Running main iterations.........................................................
..............................
Done# Model 3 - Restricted RC(1)+RL(1)+CL(1) with consistent score
# It may not be fitted with gnm package!
# Model 4 - Model 2 with consistent cells fitted exactly
freq_tab_3.1 <- freq_tab_3.1 |>
mutate(consistent.cells = factor(ifelse((R == C) & (C == L), R, 0)))
model4 <- freq_tab_3.1 |>
gnm(
Freq ~ R + C + L + consistent.cells +
Mult(1, R, C) + Mult(1, R, L) + Mult(1, C, L),
data = _,
family = poisson,
trace = F,
tolerance = 1e-12
)Initialising
Running start-up iterations..
Running main iterations................................
Done# Model 5 - Model 3 with consistent cells fitted exactly
# It may not be fitted with gnm package!
# Model 6 - Full Two-way interaction
model6 <- freq_tab_3.1 |>
gnm(
Freq ~ R * C + R * L + C * L,
data = _,
family = poisson,
trace = F,
tolerance = 1e-12
)
summary(model6)
Call:
gnm(formula = Freq ~ R * C + R * L + C * L, family = poisson,
data = freq_tab_3.1, tolerance = 1e-12, trace = F)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.685582 -0.504109 -0.008581 0.460359 1.299438
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.26029 0.11066 38.499 < 2e-16 ***
R2 -0.79606 0.18519 -4.299 1.72e-05 ***
R3 -0.09579 0.14996 -0.639 0.52299
R4 -0.43476 0.16083 -2.703 0.00687 **
C2 -1.63953 0.22829 -7.182 6.89e-13 ***
C3 -1.39351 0.20261 -6.878 6.08e-12 ***
C4 -2.33707 0.28019 -8.341 < 2e-16 ***
R2:C2 0.90120 0.28232 3.192 0.00141 **
R3:C2 0.63274 0.25192 2.512 0.01202 *
R4:C2 0.33707 0.26138 1.290 0.19719
R2:C3 -0.06534 0.29098 -0.225 0.82232
R3:C3 0.66917 0.22004 3.041 0.00236 **
R4:C3 0.60286 0.22215 2.714 0.00665 **
R2:C4 0.34284 0.36220 0.947 0.34388
R3:C4 0.56412 0.29673 1.901 0.05729 .
R4:C4 1.27643 0.28261 4.517 6.29e-06 ***
L2 -1.48502 0.22186 -6.694 2.18e-11 ***
L3 -1.52970 0.20833 -7.343 2.10e-13 ***
L4 -2.08728 0.24553 -8.501 < 2e-16 ***
R2:L2 1.34634 0.28170 4.779 1.76e-06 ***
R3:L2 1.28703 0.24739 5.203 1.97e-07 ***
R4:L2 1.03853 0.26517 3.916 8.99e-05 ***
R2:L3 0.67005 0.28444 2.356 0.01849 *
R3:L3 1.21742 0.22757 5.350 8.81e-08 ***
R4:L3 1.59671 0.23321 6.847 7.56e-12 ***
R2:L4 0.17240 0.35528 0.485 0.62750
R3:L4 0.61348 0.27149 2.260 0.02384 *
R4:L4 2.16906 0.25643 8.459 < 2e-16 ***
C2:L2 0.35033 0.21463 1.632 0.10263
C3:L2 0.22827 0.20862 1.094 0.27389
C4:L2 0.49174 0.25839 1.903 0.05703 .
C2:L3 0.65539 0.20819 3.148 0.00164 **
C3:L3 1.09826 0.18400 5.969 2.39e-09 ***
C4:L3 1.11150 0.23123 4.807 1.53e-06 ***
C2:L4 0.87747 0.23892 3.673 0.00024 ***
C3:L4 1.29956 0.20753 6.262 3.80e-10 ***
C4:L4 1.79692 0.23862 7.530 5.06e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Residual deviance: 28.998 on 27 degrees of freedom
AIC: 397.32
Number of iterations: 5# Model 7 - Model 6 with consistent cells fitted exactly
model7 <- freq_tab_3.1 |>
gnm(
Freq ~ R * C + R * L + C * L + consistent.cells,
data = _,
family = poisson,
trace = F,
tolerance = 1e-12
)
summary(model7)
Call:
gnm(formula = Freq ~ R * C + R * L + C * L + consistent.cells,
family = poisson, data = freq_tab_3.1, tolerance = 1e-12, trace = F)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.505e+00 -4.123e-01 -1.490e-08 3.008e-01 1.562e+00
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.70806 0.31642 11.719 < 2e-16 ***
R2 -0.26471 0.33752 -0.784 0.432870
R3 0.43595 0.31473 1.385 0.166007
R4 0.04506 0.31538 0.143 0.886401
C2 -1.23787 0.31573 -3.921 8.83e-05 ***
C3 -0.95472 0.29544 -3.232 0.001231 **
C4 -2.00687 0.35086 -5.720 1.07e-08 ***
R2:C2 0.49806 0.36012 1.383 0.166653
R3:C2 0.33263 0.30439 1.093 0.274482
R4:C2 0.02956 0.30915 0.096 0.923826
R2:C3 -0.37252 0.33538 -1.111 0.266684
R3:C3 0.21202 0.29412 0.721 0.470997
R4:C3 0.33724 0.27290 1.236 0.216553
R2:C4 0.06257 0.39636 0.158 0.874574
R3:C4 0.29993 0.33550 0.894 0.371334
R4:C4 1.06332 0.33951 3.132 0.001737 **
L2 -1.05809 0.31923 -3.314 0.000918 ***
L3 -1.09080 0.29894 -3.649 0.000263 ***
L4 -1.74115 0.32163 -5.413 6.18e-08 ***
R2:L2 0.91597 0.36361 2.519 0.011766 *
R3:L2 0.93125 0.31450 2.961 0.003066 **
R4:L2 0.67771 0.32813 2.065 0.038888 *
R2:L3 0.33649 0.33467 1.005 0.314694
R3:L3 0.75119 0.30376 2.473 0.013399 *
R4:L3 1.28290 0.29072 4.413 1.02e-05 ***
R2:L4 -0.12414 0.39162 -0.317 0.751243
R3:L4 0.33826 0.31651 1.069 0.285199
R4:L4 1.91005 0.31369 6.089 1.14e-09 ***
C2:L2 0.17440 0.24791 0.704 0.481744
C3:L2 0.12343 0.21731 0.568 0.570049
C4:L2 0.38519 0.26377 1.460 0.144196
C2:L3 0.54622 0.21620 2.526 0.011522 *
C3:L3 0.78217 0.22762 3.436 0.000590 ***
C4:L3 0.98203 0.23833 4.120 3.78e-05 ***
C2:L4 0.75925 0.24670 3.078 0.002087 **
C3:L4 1.14125 0.21839 5.226 1.73e-07 ***
C4:L4 1.79914 0.31902 5.640 1.70e-08 ***
consistent.cells1 0.62268 0.33657 1.850 0.064303 .
consistent.cells2 0.25993 0.39723 0.654 0.512891
consistent.cells3 0.47364 0.26552 1.784 0.074458 .
consistent.cells4 -0.16248 0.31471 -0.516 0.605666
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Residual deviance: 21.934 on 23 degrees of freedom
AIC: 398.25
Number of iterations: 5# Model 8 - Unrestricted uniform and log-multiplicative association in all partial
# association
model8 <- freq_tab_3.1 |>
gnm(
Freq ~ R + C + L +
Rscore:Cscore +
Mult(1, R, C) +
Rscore:Lscore +
Mult(1, R, L) +
Cscore:Lscore +
Mult(1, C, L),
data = _,
family = poisson,
trace = F,
tolerance = 1e-12
)Initialising
Running start-up iterations..
Running main iterations.........................................................
.....................
Donesummary(model8)
Call:
gnm(formula = Freq ~ R + C + L + Rscore:Cscore + Mult(1, R, C) +
Rscore:Lscore + Mult(1, R, L) + Cscore:Lscore + Mult(1, C,
L), family = poisson, data = freq_tab_3.1, tolerance = 1e-12, trace = F)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.64348 -0.68031 -0.02712 0.48954 2.01814
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.598847 NA NA NA
R2 -0.980385 NA NA NA
R3 -0.767847 NA NA NA
R4 -1.592697 NA NA NA
C2 -1.611171 NA NA NA
C3 -2.094217 NA NA NA
C4 -3.433621 NA NA NA
L2 -1.452544 NA NA NA
L3 -2.267652 NA NA NA
L4 -3.811423 NA NA NA
Rscore:Cscore 0.129644 0.027279 4.752 2.01e-06 ***
Mult(., R, C). 0.238257 NA NA NA
Mult(1, ., C).R1 -0.276217 NA NA NA
Mult(1, ., C).R2 0.604742 NA NA NA
Mult(1, ., C).R3 0.267895 NA NA NA
Mult(1, ., C).R4 -0.349098 NA NA NA
Mult(1, R, .).C1 -0.775191 NA NA NA
Mult(1, R, .).C2 2.514695 NA NA NA
Mult(1, R, .).C3 -0.110731 NA NA NA
Mult(1, R, .).C4 -2.507478 NA NA NA
Rscore:Lscore 0.254496 0.028029 9.080 < 2e-16 ***
Mult(., R, L). -0.605608 NA NA NA
Mult(1, ., L).R1 2.789803 NA NA NA
Mult(1, ., L).R2 -2.835243 NA NA NA
Mult(1, ., L).R3 -2.228299 NA NA NA
Mult(1, ., L).R4 1.676905 NA NA NA
Mult(1, R, .).L1 -0.006301 NA NA NA
Mult(1, R, .).L2 0.268441 NA NA NA
Mult(1, R, .).L3 0.065107 NA NA NA
Mult(1, R, .).L4 -0.275010 NA NA NA
Cscore:Lscore 0.201725 0.025855 7.802 < 2e-16 ***
Mult(., C, L). -0.540930 NA NA NA
Mult(1, ., L).C1 0.147910 NA NA NA
Mult(1, ., L).C2 -0.122908 NA NA NA
Mult(1, ., L).C3 -0.862245 NA NA NA
Mult(1, ., L).C4 0.206794 NA NA NA
Mult(1, C, .).L1 -0.018366 NA NA NA
Mult(1, C, .).L2 -0.306632 NA NA NA
Mult(1, C, .).L3 0.549238 NA NA NA
Mult(1, C, .).L4 0.028881 NA NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Std. Error is NA where coefficient has been constrained or is unidentified
Residual deviance: 45.854 on 36 degrees of freedom
AIC: 396.17
Number of iterations: 78pickCoef(model8, "[.]R")Mult(1, ., C).R1 Mult(1, ., C).R2 Mult(1, ., C).R3 Mult(1, ., C).R4
13 14 15 16
Mult(1, ., L).R1 Mult(1, ., L).R2 Mult(1, ., L).R3 Mult(1, ., L).R4
23 24 25 26 pickCoef(model8, "[.]C")Mult(1, R, .).C1 Mult(1, R, .).C2 Mult(1, R, .).C3 Mult(1, R, .).C4
17 18 19 20
Mult(1, ., L).C1 Mult(1, ., L).C2 Mult(1, ., L).C3 Mult(1, ., L).C4
33 34 35 36 pickCoef(model8, "[.]L")Mult(1, R, .).L1 Mult(1, R, .).L2 Mult(1, R, .).L3 Mult(1, R, .).L4
27 28 29 30
Mult(1, C, .).L1 Mult(1, C, .).L2 Mult(1, C, .).L3 Mult(1, C, .).L4
37 38 39 40 data.frame(var = names(model8$coefficients),
estimate = model8$coefficients) |>
mutate(estimate = estimate,
number = row_number()) var estimate number
(Intercept) (Intercept) 3.59884710 1
R2 R2 -0.98038545 2
R3 R3 -0.76784741 3
R4 R4 -1.59269723 4
C2 C2 -1.61117117 5
C3 C3 -2.09421716 6
C4 C4 -3.43362111 7
L2 L2 -1.45254414 8
L3 L3 -2.26765245 9
L4 L4 -3.81142296 10
Rscore:Cscore Rscore:Cscore 0.12964400 11
Mult(., R, C). Mult(., R, C). 0.23825669 12
Mult(1, ., C).R1 Mult(1, ., C).R1 -0.27621726 13
Mult(1, ., C).R2 Mult(1, ., C).R2 0.60474220 14
Mult(1, ., C).R3 Mult(1, ., C).R3 0.26789527 15
Mult(1, ., C).R4 Mult(1, ., C).R4 -0.34909784 16
Mult(1, R, .).C1 Mult(1, R, .).C1 -0.77519131 17
Mult(1, R, .).C2 Mult(1, R, .).C2 2.51469493 18
Mult(1, R, .).C3 Mult(1, R, .).C3 -0.11073136 19
Mult(1, R, .).C4 Mult(1, R, .).C4 -2.50747823 20
Rscore:Lscore Rscore:Lscore 0.25449639 21
Mult(., R, L). Mult(., R, L). -0.60560791 22
Mult(1, ., L).R1 Mult(1, ., L).R1 2.78980334 23
Mult(1, ., L).R2 Mult(1, ., L).R2 -2.83524309 24
Mult(1, ., L).R3 Mult(1, ., L).R3 -2.22829891 25
Mult(1, ., L).R4 Mult(1, ., L).R4 1.67690495 26
Mult(1, R, .).L1 Mult(1, R, .).L1 -0.00630102 27
Mult(1, R, .).L2 Mult(1, R, .).L2 0.26844081 28
Mult(1, R, .).L3 Mult(1, R, .).L3 0.06510732 29
Mult(1, R, .).L4 Mult(1, R, .).L4 -0.27500999 30
Cscore:Lscore Cscore:Lscore 0.20172487 31
Mult(., C, L). Mult(., C, L). -0.54093005 32
Mult(1, ., L).C1 Mult(1, ., L).C1 0.14790966 33
Mult(1, ., L).C2 Mult(1, ., L).C2 -0.12290819 34
Mult(1, ., L).C3 Mult(1, ., L).C3 -0.86224538 35
Mult(1, ., L).C4 Mult(1, ., L).C4 0.20679354 36
Mult(1, C, .).L1 Mult(1, C, .).L1 -0.01836573 37
Mult(1, C, .).L2 Mult(1, C, .).L2 -0.30663170 38
Mult(1, C, .).L3 Mult(1, C, .).L3 0.54923804 39
Mult(1, C, .).L4 Mult(1, C, .).L4 0.02888076 40# mu1: 13-16
mu1 <- getContrasts(
model8,
pickCoef(model8, "[.]R")[1:4],
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit"
)
mu1 Estimate Std. Error
Mult(1, ., C).R1 -0.4291986 0.18120993
Mult(1, ., C).R2 0.6893016 0.09977232
Mult(1, ., C).R3 0.2616277 0.15931019
Mult(1, ., C).R4 -0.5217306 0.15976113# nu1: 17-20
nu1 <- getContrasts(
model8,
pickCoef(model8, "[.]C")[1:4],
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit"
)
nu1 Estimate Std. Error
Mult(1, R, .).C1 -0.15388700 0.15383900
Mult(1, R, .).C2 0.75746716 0.08227847
Mult(1, R, .).C3 0.03017965 0.15480459
Mult(1, R, .).C4 -0.63375980 0.11651692# mu2: 23-26
mu2 <- getContrasts(
model8,
pickCoef(model8, "[.]R")[5:8],
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit"
)
mu2 Estimate Std. Error
Mult(1, ., L).R1 0.6061444 0.08512883
Mult(1, ., L).R2 -0.5539702 0.08877502
Mult(1, ., L).R3 -0.4287935 0.09713147
Mult(1, ., L).R4 0.3766193 0.10554742# eta1: 27-30
eta1 <- getContrasts(
model8,
pickCoef(model8, "[.]L")[1:4],
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit"
)
eta1 Estimate Std. Error
Mult(1, R, .).L1 -0.04977493 0.08992116
Mult(1, R, .).L2 0.65658058 0.06234704
Mult(1, R, .).L3 0.13381443 0.08765352
Mult(1, R, .).L4 -0.74062008 0.04970264# nu2: 33-36
nu2 <- getContrasts(
model8,
pickCoef(model8, "[.]C")[5:8],
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit"
)
nu2 Estimate Std. Error
Mult(1, ., L).C1 0.35910324 0.2995712
Mult(1, ., L).C2 0.04079069 0.3606231
Mult(1, ., L).C3 -0.82820782 0.1075525
Mult(1, ., L).C4 0.42831389 0.3177968# eta2: 37-40
eta2 <- getContrasts(
model8,
pickCoef(model8, "[.]L")[5:8],
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit"
)
eta2 Estimate Std. Error
Mult(1, C, .).L1 -0.13230141 0.3592002
Mult(1, C, .).L2 -0.59941504 0.2494518
Mult(1, C, .).L3 0.78745843 0.1409196
Mult(1, C, .).L4 -0.05574198 0.3241035con <- c(mu1$qvframe[, 1][c(1,4)],
nu1$qvframe[, 1][c(1,4)],
mu2$qvframe[, 1][c(1,4)],
eta1$qvframe[, 1][c(1,4)],
nu2$qvframe[, 1][c(1,4)],
eta2$qvframe[, 1][c(1,4)])
model8_se <- freq_tab_3.1 |>
gnm(
Freq ~ R + C + L +
Rscore:Cscore +
Mult(1, R, C) +
Rscore:Lscore +
Mult(1, R, L) +
Cscore:Lscore +
Mult(1, C, L),
# constrain = c(13,16,17,20,23,26,27,30,33,36,37,40),
# constrainTo = con,
data = _,
family = poisson,
trace = F,
tolerance = 1e-15
)Initialising
Running start-up iterations..
Running main iterations.........................................................
................................................................................
................................................................................
................................................................................
................................................................................
................................................................................
...........................................
DoneWarning in gnmFit(modelTools = modelTools, y = y, constrain = constrain, : Fitting algorithm has either not converged or converged
to a non-solution of the likelihood equations.
Use exitInfo() for numerical details of last iteration.summary(model8_se)
Call:
gnm(formula = Freq ~ R + C + L + Rscore:Cscore + Mult(1, R, C) +
Rscore:Lscore + Mult(1, R, L) + Cscore:Lscore + Mult(1, C,
L), family = poisson, data = freq_tab_3.1, tolerance = 1e-15, trace = F)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.64348 -0.68031 -0.02712 0.48954 2.01814
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.594392 NA NA NA
R2 -0.977519 NA NA NA
R3 -0.761181 NA NA NA
R4 -1.587929 NA NA NA
C2 -1.613462 NA NA NA
C3 -2.038526 NA NA NA
C4 -3.427495 NA NA NA
L2 -1.439324 NA NA NA
L3 -2.287543 NA NA NA
L4 -3.815832 NA NA NA
Rscore:Cscore 0.129644 0.027279 4.752 2.01e-06 ***
Mult(., R, C). -0.338925 NA NA NA
Mult(1, ., C).R1 1.170488 NA NA NA
Mult(1, ., C).R2 -2.905396 NA NA NA
Mult(1, ., C).R3 -1.346926 NA NA NA
Mult(1, ., C).R4 1.507681 NA NA NA
Mult(1, R, .).C1 -0.106945 NA NA NA
Mult(1, R, .).C2 0.392924 NA NA NA
Mult(1, R, .).C3 -0.005986 NA NA NA
Mult(1, R, .).C4 -0.370150 NA NA NA
Rscore:Lscore 0.254496 0.028029 9.080 < 2e-16 ***
Mult(., R, L). -0.267323 NA NA NA
Mult(1, ., L).R1 -6.238651 NA NA NA
Mult(1, ., L).R2 6.266066 NA NA NA
Mult(1, ., L).R3 4.916803 NA NA NA
Mult(1, ., L).R4 -3.764631 NA NA NA
Mult(1, R, .).L1 0.011758 NA NA NA
Mult(1, R, .).L2 -0.268226 NA NA NA
Mult(1, R, .).L3 -0.061013 NA NA NA
Mult(1, R, .).L4 0.285593 NA NA NA
Cscore:Lscore 0.201725 0.025855 7.802 < 2e-16 ***
Mult(., C, L). -0.237549 NA NA NA
Mult(1, ., L).C1 -0.145041 NA NA NA
Mult(1, ., L).C2 0.341127 NA NA NA
Mult(1, ., L).C3 1.668373 NA NA NA
Mult(1, ., L).C4 -0.250748 NA NA NA
Mult(1, C, .).L1 0.161115 NA NA NA
Mult(1, C, .).L2 0.526771 NA NA NA
Mult(1, C, .).L3 -0.558873 NA NA NA
Mult(1, C, .).L4 0.101184 NA NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Std. Error is NA where coefficient has been constrained or is unidentified
Residual deviance: 45.854 on 36 degrees of freedom
AIC: 396.17
Number of iterations: 500options(contrasts = c(factor="contr.treatment",
ordered="contr.treatment"))
# Model 9 - Model 8 with consistent row scores
model9 <- freq_tab_3.1 |>
gnm(
Freq ~ R + C + L +
Rscore:Cscore +
Mult(1, R, C + L) +
Rscore:Lscore +
Cscore:Lscore +
Mult(1, C, L),
data = _,
family = poisson,
trace = F,
tolerance = 1e-12
)Initialising
Running start-up iterations..
Running main iterations.........................................................
Donesummary(model9)
Call:
gnm(formula = Freq ~ R + C + L + Rscore:Cscore + Mult(1, R, C +
L) + Rscore:Lscore + Cscore:Lscore + Mult(1, C, L), family = poisson,
data = freq_tab_3.1, tolerance = 1e-12, trace = F)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.56526 -0.67024 -0.02438 0.54037 1.98333
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.59143 NA NA NA
R2 -1.00021 NA NA NA
R3 -0.76239 NA NA NA
R4 -1.56089 NA NA NA
C2 -1.62160 NA NA NA
C3 -2.09063 NA NA NA
C4 -3.46410 NA NA NA
L2 -1.43836 NA NA NA
L3 -2.22931 NA NA NA
L4 -3.77811 NA NA NA
Rscore:Cscore 0.13353 0.02535 5.266 1.39e-07 ***
Mult(., R, C + L). 0.41490 NA NA NA
Mult(1, ., C + L).R1 -1.25113 NA NA NA
Mult(1, ., C + L).R2 1.54126 NA NA NA
Mult(1, ., C + L).R3 1.05354 NA NA NA
Mult(1, ., C + L).R4 -0.91803 NA NA NA
Mult(1, R, . + L).C1 -0.14160 NA NA NA
Mult(1, R, . + L).C2 0.44885 NA NA NA
Mult(1, R, . + L).C3 0.01394 NA NA NA
Mult(1, R, . + L).C4 -0.46676 NA NA NA
Mult(1, R, C + .).L2 0.80669 NA NA NA
Mult(1, R, C + .).L3 0.20043 NA NA NA
Mult(1, R, C + .).L4 -0.77812 NA NA NA
Rscore:Lscore 0.24848 0.02722 9.127 < 2e-16 ***
Cscore:Lscore 0.20270 0.02562 7.913 < 2e-16 ***
Mult(., C, L). -0.81663 NA NA NA
Mult(1, ., L).C1 0.11393 NA NA NA
Mult(1, ., L).C2 -0.04445 NA NA NA
Mult(1, ., L).C3 -0.64049 NA NA NA
Mult(1, ., L).C4 0.16394 NA NA NA
Mult(1, C, .).L1 -0.04492 NA NA NA
Mult(1, C, .).L2 -0.33272 NA NA NA
Mult(1, C, .).L3 0.42413 NA NA NA
Mult(1, C, .).L4 -0.03133 NA NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Std. Error is NA where coefficient has been constrained or is unidentified
Residual deviance: 47.134 on 38 degrees of freedom
AIC: 393.45
Number of iterations: 57mu <- getContrasts(
model9,
pickCoef(model9, "[.]R"),
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit"
)
nu2 <- getContrasts(
model9,
pickCoef(model9, "[.]C")[5:8],
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit"
)
eta2 <- getContrasts(
model9,
pickCoef(model9, "[.]L")[4:7],
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit"
)
model9.extended <-
freq_tab_3.1 |> gnm(
Freq ~ R + C + L +
Rscore:Cscore +
Mult(1, R, C) +
Rscore:Lscore +
Mult(1, R, L) +
Cscore:Lscore +
Mult(1, C, L),
constrain = c(13:16, 23:26),
constrainTo = c(mu$qvframe[, 1], mu$qvframe[, 1]),
data = _,
family = poisson,
trace = F,
tolerance = 1e-12
)Initialising
Running start-up iterations..
Running main iterations.........................................................
............
Donesummary(model9.extended)
Call:
gnm(formula = Freq ~ R + C + L + Rscore:Cscore + Mult(1, R, C) +
Rscore:Lscore + Mult(1, R, L) + Cscore:Lscore + Mult(1, C,
L), constrain = c(13:16, 23:26), constrainTo = c(mu$qvframe[,
1], mu$qvframe[, 1]), family = poisson, data = freq_tab_3.1,
tolerance = 1e-12, trace = F)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.56526 -0.67024 -0.02438 0.54037 1.98333
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.59300 NA NA NA
R2 -0.99919 NA NA NA
R3 -0.76155 NA NA NA
R4 -1.56077 NA NA NA
C2 -1.58393 NA NA NA
C3 -2.02848 NA NA NA
C4 -3.48212 NA NA NA
L2 -1.36713 NA NA NA
L3 -2.27850 NA NA NA
L4 -3.81414 NA NA NA
Rscore:Cscore 0.13353 0.02520 5.300 1.16e-07 ***
Mult(., R, C). -1.31709 NA NA NA
Mult(1, ., C).R1 -0.56136 NA NA NA
Mult(1, ., C).R2 0.59333 NA NA NA
Mult(1, ., C).R3 0.39165 NA NA NA
Mult(1, ., C).R4 -0.42362 NA NA NA
Mult(1, R, .).C1 0.08594 NA NA NA
Mult(1, R, .).C2 -0.36387 NA NA NA
Mult(1, R, .).C3 -0.03255 NA NA NA
Mult(1, R, .).C4 0.33364 NA NA NA
Rscore:Lscore 0.24848 0.02521 9.856 < 2e-16 ***
Mult(., R, L). 0.56373 NA NA NA
Mult(1, ., L).R1 -0.56136 NA NA NA
Mult(1, ., L).R2 0.59333 NA NA NA
Mult(1, ., L).R3 0.39165 NA NA NA
Mult(1, ., L).R4 -0.42362 NA NA NA
Mult(1, R, .).L1 -0.05283 NA NA NA
Mult(1, R, .).L2 1.38296 NA NA NA
Mult(1, R, .).L3 0.30391 NA NA NA
Mult(1, R, .).L4 -1.43778 NA NA NA
Cscore:Lscore 0.20270 0.02560 7.919 < 2e-16 ***
Mult(., C, L). 0.43424 NA NA NA
Mult(1, ., L).C1 0.05769 NA NA NA
Mult(1, ., L).C2 0.30073 NA NA NA
Mult(1, ., L).C3 1.21534 NA NA NA
Mult(1, ., L).C4 -0.01904 NA NA NA
Mult(1, C, .).L1 -0.16502 NA NA NA
Mult(1, C, .).L2 -0.51774 NA NA NA
Mult(1, C, .).L3 0.40983 NA NA NA
Mult(1, C, .).L4 -0.14836 NA NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Std. Error is NA where coefficient has been constrained or is unidentified
Residual deviance: 47.134 on 40 degrees of freedom
AIC: 389.45
Number of iterations: 69nu1 <- getContrasts(
model9.extended,
pickCoef(model9.extended, "[.]C")[1:4],
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit"
)
eta1 <- getContrasts(
model9.extended,
pickCoef(model9.extended, "[.]L")[1:4],
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit"
)
model9.phis <-
freq_tab_3.1 |> gnm(
Freq ~ R + C + L +
Rscore:Cscore +
Mult(1, R, C) +
Rscore:Lscore +
Mult(1, R, L) +
Cscore:Lscore +
Mult(1, C, L),
constrain = c(13:20, 23:30, 33:40),
constrainTo = c(
mu$qvframe[, 1],
nu1$qvframe[, 1],
mu$qvframe[, 1],
eta1$qvframe[, 1],
nu2$qvframe[, 1],
eta2$qvframe[, 1]
),
data = _,
family = poisson,
trace = F,
tolerance = 1e-12
)Initialising
Running main iterations.
Donesummary(model9.phis)
Call:
gnm(formula = Freq ~ R + C + L + Rscore:Cscore + Mult(1, R, C) +
Rscore:Lscore + Mult(1, R, L) + Cscore:Lscore + Mult(1, C,
L), constrain = c(13:20, 23:30, 33:40), constrainTo = c(mu$qvframe[,
1], nu1$qvframe[, 1], mu$qvframe[, 1], eta1$qvframe[, 1],
nu2$qvframe[, 1], eta2$qvframe[, 1]), family = poisson, data = freq_tab_3.1,
tolerance = 1e-12, trace = F)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.56526 -0.67024 -0.02438 0.54037 1.98333
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.56903 0.08137 43.862 < 2e-16 ***
R2 -0.97605 0.11881 -8.215 < 2e-16 ***
R3 -0.74245 0.14810 -5.013 5.36e-07 ***
R4 -1.55801 0.22175 -7.026 2.13e-12 ***
C2 -1.59504 0.10577 -15.081 < 2e-16 ***
C3 -2.08142 0.17240 -12.073 < 2e-16 ***
C4 -3.47861 0.26841 -12.960 < 2e-16 ***
L2 -1.42666 0.10739 -13.285 < 2e-16 ***
L3 -2.18148 0.17791 -12.262 < 2e-16 ***
L4 -3.81133 0.27937 -13.643 < 2e-16 ***
Rscore:Cscore 0.13353 0.02486 5.371 7.83e-08 ***
Mult(., R, C). -0.66121 0.14008 -4.720 2.36e-06 ***
Mult(1, ., C).R1 -0.56136 NA NA NA
Mult(1, ., C).R2 0.59333 NA NA NA
Mult(1, ., C).R3 0.39165 NA NA NA
Mult(1, ., C).R4 -0.42362 NA NA NA
Mult(1, R, .).C1 0.15965 NA NA NA
Mult(1, R, .).C2 -0.73633 NA NA NA
Mult(1, R, .).C3 -0.07638 NA NA NA
Mult(1, R, .).C4 0.65306 NA NA NA
Rscore:Lscore 0.24848 0.02515 9.879 < 2e-16 ***
Mult(., R, L). 1.13662 0.13964 8.140 < 2e-16 ***
Mult(1, ., L).R1 -0.56136 NA NA NA
Mult(1, ., L).R2 0.59333 NA NA NA
Mult(1, ., L).R3 0.39165 NA NA NA
Mult(1, ., L).R4 -0.42362 NA NA NA
Mult(1, R, .).L1 -0.05054 NA NA NA
Mult(1, R, .).L2 0.66157 NA NA NA
Mult(1, R, .).L3 0.12640 NA NA NA
Mult(1, R, .).L4 -0.73743 NA NA NA
Cscore:Lscore 0.20270 0.02315 8.756 < 2e-16 ***
Mult(., C, L). -0.28351 0.10068 -2.816 0.00486 **
Mult(1, ., L).C1 0.33661 NA NA NA
Mult(1, ., L).C2 0.08944 NA NA NA
Mult(1, ., L).C3 -0.84069 NA NA NA
Mult(1, ., L).C4 0.41464 NA NA NA
Mult(1, C, .).L1 -0.08991 NA NA NA
Mult(1, C, .).L2 -0.62111 NA NA NA
Mult(1, C, .).L3 0.77584 NA NA NA
Mult(1, C, .).L4 -0.06482 NA NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Std. Error is NA where coefficient has been constrained or is unidentified
Residual deviance: 47.134 on 48 degrees of freedom
AIC: 373.45
Number of iterations: 1mu Estimate Std. Error
Mult(1, ., C + L).R1 -0.5613577 0.07946201
Mult(1, ., C + L).R2 0.5933253 0.06976100
Mult(1, ., C + L).R3 0.3916506 0.08375491
Mult(1, ., C + L).R4 -0.4236182 0.08896023nu1 Estimate Std. Error
Mult(1, R, .).C1 0.15964986 0.14637516
Mult(1, R, .).C2 -0.73633460 0.08871449
Mult(1, R, .).C3 -0.07637713 0.15684115
Mult(1, R, .).C4 0.65306187 0.10868144eta1 Estimate Std. Error
Mult(1, R, .).L1 -0.05053742 0.08981862
Mult(1, R, .).L2 0.66156822 0.06038098
Mult(1, R, .).L3 0.12639555 0.08679365
Mult(1, R, .).L4 -0.73742635 0.04940126nu2 Estimate Std. Error
Mult(1, ., L).C1 0.33661004 0.30087703
Mult(1, ., L).C2 0.08944163 0.36058729
Mult(1, ., L).C3 -0.84069369 0.08936905
Mult(1, ., L).C4 0.41464201 0.32252708eta2 Estimate Std. Error
Mult(1, C, .).L1 -0.08990831 0.3610226
Mult(1, C, .).L2 -0.62111332 0.2399524
Mult(1, C, .).L3 0.77584334 0.1519881
Mult(1, C, .).L4 -0.06482171 0.3227342model.summary(model9)# A tibble: 1 × 6
`Model Description` df L2 BIC Delta p
<chr> <int> <dbl> <dbl> <dbl> <dbl>
1 model9 38 47.1 -231. 6.25 0.147# Model 10 - Model 8 with consistent col scores
model10 <-
freq_tab_3.1 |> gnm(
Freq ~ R + C + L +
Rscore:Cscore +
Mult(1, C, R + L) +
Rscore:Lscore +
Mult(1, R, L) +
Cscore:Lscore,
data = _,
family = poisson,
trace = F,
tolerance = 1e-12
)Initialising
Running start-up iterations..
Running main iterations.........................................................
..........................................................................
Donemodel.summary(model10)# A tibble: 1 × 6
`Model Description` df L2 BIC Delta p
<chr> <int> <dbl> <dbl> <dbl> <dbl>
1 model10 38 50.9 -227. 6.31 0.0788# Model 11 - Model 8 with consistent layer scores
model11 <-
freq_tab_3.1 |> gnm(
Freq ~ R + C + L +
Rscore:Cscore +
Mult(1, L, R + C) +
Rscore:Lscore +
Mult(1, R, C) +
Cscore:Lscore,
data = _,
family = poisson,
trace = F,
tolerance = 1e-12
)Initialising
Running start-up iterations..
Running main iterations.........................................................
................................
Donemodel.summary(model11)# A tibble: 1 × 6
`Model Description` df L2 BIC Delta p
<chr> <int> <dbl> <dbl> <dbl> <dbl>
1 model11 38 53.6 -225. 6.58 0.0477# Model 12 - Modek 9 + U(RL)=U(CL)
# To impose equality constraint on U(RL) & U(CL),
# a covariate is created in following manner:
freq_tab_3.1 <- freq_tab_3.1|>
mutate(cov_rcl = Rscore * Lscore + Cscore * Lscore)
model12 <-
freq_tab_3.1 |> gnm(
Freq ~ R + C + L +
Rscore:Cscore +
Mult(1, R, C + L) +
cov_rcl +
Mult(1, C, L),
data = _,
family = poisson,
trace = F,
tolerance = 1e-12
)Initialising
Running start-up iterations..
Running main iterations.........................................................
.......................................
Done#summary(model12)
model.summary(model12)# A tibble: 1 × 6
`Model Description` df L2 BIC Delta p
<chr> <int> <dbl> <dbl> <dbl> <dbl>
1 model12 39 48.6 -237. 6.57 0.140mu <- getContrasts(
model12,
pickCoef(model12, "[.]R"),
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit"
)
nu2 <-
getContrasts(
model12,
pickCoef(model12, "[.]C")[5:8],
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit"
)
eta2 <-
getContrasts(
model12,
pickCoef(model12, "[.]L")[4:7],
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit"
)
model12.extended <-
freq_tab_3.1 |>
gnm(Freq ~ R + C + L +
Rscore:Cscore +
Mult(1, R, C) +
cov_rcl +
Mult(1, R, L) +
Mult(1, C, L),
constrain = c(13:16, 23:26),
constrainTo = c(mu$qvframe[, 1], mu$qvframe[, 1]),
data = _,
family = poisson,
trace = F,
tolerance = 1e-12
)Initialising
Running start-up iterations..
Running main iterations.........................................................
................................................
Donesummary(model12.extended)
Call:
gnm(formula = Freq ~ R + C + L + Rscore:Cscore + Mult(1, R, C) +
cov_rcl + Mult(1, R, L) + Mult(1, C, L), constrain = c(13:16,
23:26), constrainTo = c(mu$qvframe[, 1], mu$qvframe[, 1]),
family = poisson, data = freq_tab_3.1, tolerance = 1e-12, trace = F)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.58844 -0.76809 0.02278 0.56330 1.92999
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.57371 NA NA NA
R2 -0.98116 NA NA NA
R3 -0.68473 NA NA NA
R4 -1.39730 NA NA NA
C2 -1.63413 NA NA NA
C3 -2.15384 NA NA NA
C4 -3.65708 NA NA NA
L2 -1.39412 NA NA NA
L3 -2.13584 NA NA NA
L4 -3.74767 NA NA NA
Rscore:Cscore 0.13336 0.02542 5.246 1.55e-07 ***
Mult(., R, C). 1.10155 NA NA NA
Mult(1, ., C).R1 -0.53827 NA NA NA
Mult(1, ., C).R2 0.59766 NA NA NA
Mult(1, ., C).R3 0.38941 NA NA NA
Mult(1, ., C).R4 -0.44880 NA NA NA
Mult(1, R, .).C1 -0.09323 NA NA NA
Mult(1, R, .).C2 0.44282 NA NA NA
Mult(1, R, .).C3 0.05887 NA NA NA
Mult(1, R, .).C4 -0.36049 NA NA NA
cov_rcl 0.22423 0.01709 13.123 < 2e-16 ***
Mult(., R, L). -1.84369 NA NA NA
Mult(1, ., L).R1 -0.53827 NA NA NA
Mult(1, ., L).R2 0.59766 NA NA NA
Mult(1, ., L).R3 0.38941 NA NA NA
Mult(1, ., L).R4 -0.44880 NA NA NA
Mult(1, R, .).L1 0.01040 NA NA NA
Mult(1, R, .).L2 -0.42684 NA NA NA
Mult(1, R, .).L3 -0.09320 NA NA NA
Mult(1, R, .).L4 0.43884 NA NA NA
Mult(., C, L). 0.63123 NA NA NA
Mult(1, ., L).C1 -0.25933 NA NA NA
Mult(1, ., L).C2 -0.12167 NA NA NA
Mult(1, ., L).C3 0.86218 NA NA NA
Mult(1, ., L).C4 -0.41136 NA NA NA
Mult(1, C, .).L1 -0.03042 NA NA NA
Mult(1, C, .).L2 -0.32432 NA NA NA
Mult(1, C, .).L3 0.29074 NA NA NA
Mult(1, C, .).L4 -0.08827 NA NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Std. Error is NA where coefficient has been constrained or is unidentified
Residual deviance: 48.569 on 41 degrees of freedom
AIC: 388.89
Number of iterations: 105nu1 <- getContrasts(
model12.extended,
pickCoef(model12.extended, "[.]C")[1:4],
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit"
)
eta1 <-
getContrasts(
model12.extended,
pickCoef(model12.extended, "[.]L")[1:4],
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit"
)
model12.phis <-
freq_tab_3.1 |> gnm(
Freq ~ R + C + L + Rscore:Cscore + Mult(1, R, C) +
cov_rcl + Mult(1, R, L) + Mult(1, C, L),
constrain = c(13:20, 23:30, 32:39),
constrainTo = c(mu$qvframe[, 1], nu1$qvframe[, 1],
mu$qvframe[, 1], eta1$qvframe[, 1],
nu2$qvframe[, 1], eta2$qvframe[, 1]),
data = _,
family = poisson,
trace = F,
tolerance = 1e-12)Initialising
Running main iterations.
Donesummary(model12.phis)
Call:
gnm(formula = Freq ~ R + C + L + Rscore:Cscore + Mult(1, R, C) +
cov_rcl + Mult(1, R, L) + Mult(1, C, L), constrain = c(13:20,
23:30, 32:39), constrainTo = c(mu$qvframe[, 1], nu1$qvframe[,
1], mu$qvframe[, 1], eta1$qvframe[, 1], nu2$qvframe[, 1],
eta2$qvframe[, 1]), family = poisson, data = freq_tab_3.1,
tolerance = 1e-12, trace = F)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.58844 -0.76809 0.02278 0.56330 1.92999
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.55535 0.08106 43.859 < 2e-16 ***
R2 -0.92909 0.11339 -8.194 < 2e-16 ***
R3 -0.64221 0.12609 -5.093 3.52e-07 ***
R4 -1.39320 0.17746 -7.851 < 2e-16 ***
C2 -1.63743 0.10043 -16.304 < 2e-16 ***
C3 -2.18079 0.15608 -13.973 < 2e-16 ***
C4 -3.65343 0.23846 -15.321 < 2e-16 ***
L2 -1.39736 0.10537 -13.261 < 2e-16 ***
L3 -2.13230 0.17235 -12.372 < 2e-16 ***
L4 -3.74831 0.27154 -13.804 < 2e-16 ***
Rscore:Cscore 0.13336 0.02490 5.357 8.48e-08 ***
Mult(., R, C). 0.64007 0.13725 4.663 3.11e-06 ***
Mult(1, ., C).R1 -0.53827 NA NA NA
Mult(1, ., C).R2 0.59766 NA NA NA
Mult(1, ., C).R3 0.38941 NA NA NA
Mult(1, ., C).R4 -0.44880 NA NA NA
Mult(1, R, .).C1 -0.18108 NA NA NA
Mult(1, R, .).C2 0.74146 NA NA NA
Mult(1, R, .).C3 0.08067 NA NA NA
Mult(1, R, .).C4 -0.64105 NA NA NA
cov_rcl 0.22423 0.01636 13.707 < 2e-16 ***
Mult(., R, L). -1.13998 0.13669 -8.340 < 2e-16 ***
Mult(1, ., L).R1 -0.53827 NA NA NA
Mult(1, ., L).R2 0.59766 NA NA NA
Mult(1, ., L).R3 0.38941 NA NA NA
Mult(1, ., L).R4 -0.44880 NA NA NA
Mult(1, R, .).L1 0.04544 NA NA NA
Mult(1, R, .).L2 -0.66170 NA NA NA
Mult(1, R, .).L3 -0.12210 NA NA NA
Mult(1, R, .).L4 0.73836 NA NA NA
Mult(., C, L). -0.27613 0.10158 -2.718 0.00656 **
Mult(1, ., L).C1 0.27770 NA NA NA
Mult(1, ., L).C2 0.13959 NA NA NA
Mult(1, ., L).C3 -0.84752 NA NA NA
Mult(1, ., L).C4 0.43024 NA NA NA
Mult(1, C, .).L1 0.01743 NA NA NA
Mult(1, C, .).L2 -0.65220 NA NA NA
Mult(1, C, .).L3 0.74916 NA NA NA
Mult(1, C, .).L4 -0.11439 NA NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Std. Error is NA where coefficient has been constrained or is unidentified
Residual deviance: 48.569 on 49 degrees of freedom
AIC: 372.89
Number of iterations: 1mu Estimate Std. Error
Mult(1, ., C + L).R1 -0.5382715 0.07940991
Mult(1, ., C + L).R2 0.5976649 0.06869502
Mult(1, ., C + L).R3 0.3894076 0.08390695
Mult(1, ., C + L).R4 -0.4488010 0.08495489nu1 Estimate Std. Error
Mult(1, R, .).C1 -0.18108336 0.1505749
Mult(1, R, .).C2 0.74145909 0.0888159
Mult(1, R, .).C3 0.08066998 0.1598862
Mult(1, R, .).C4 -0.64104571 0.1141041eta1 Estimate Std. Error
Mult(1, R, .).L1 0.0454423 0.08895571
Mult(1, R, .).L2 -0.6617033 0.05943373
Mult(1, R, .).L3 -0.1220991 0.08507409
Mult(1, R, .).L4 0.7383601 0.04862532nu2 Estimate Std. Error
Mult(1, ., L).C1 0.2776974 0.30818568
Mult(1, ., L).C2 0.1395889 0.36977790
Mult(1, ., L).C3 -0.8475231 0.08001496
Mult(1, ., L).C4 0.4302368 0.32326672eta2 Estimate Std. Error
Mult(1, C, .).L1 0.01742597 0.3681779
Mult(1, C, .).L2 -0.65219849 0.2296758
Mult(1, C, .).L3 0.74916478 0.1797929
Mult(1, C, .).L4 -0.11439226 0.3321123model.summary(model12)# A tibble: 1 × 6
`Model Description` df L2 BIC Delta p
<chr> <int> <dbl> <dbl> <dbl> <dbl>
1 model12 39 48.6 -237. 6.57 0.140# Model 15
freq_tab_3.1 <-
freq_tab_3.1 |> mutate(cov_rcl = Rscore * Lscore + Cscore * Lscore)
model15 <- freq_tab_3.1 |>
gnm(
Freq ~ R + C + L + Rscore:Cscore + Mult(1, R, C + L) + cov_rcl,
data = _,
family = poisson,
trace = F,
tolerance = 1e-12
)Initialising
Running start-up iterations..
Running main iterations....................................................
Donesummary(model15)
Call:
gnm(formula = Freq ~ R + C + L + Rscore:Cscore + Mult(1, R, C +
L) + cov_rcl, family = poisson, data = freq_tab_3.1, tolerance = 1e-12,
trace = F)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.98677 -0.75501 0.01756 0.52206 2.97799
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.577949 NA NA NA
R2 -0.959160 NA NA NA
R3 -0.676834 NA NA NA
R4 -1.422370 NA NA NA
C2 -1.700522 NA NA NA
C3 -2.184681 NA NA NA
C4 -3.668341 NA NA NA
L2 -1.487671 NA NA NA
L3 -2.177105 NA NA NA
L4 -3.748632 NA NA NA
Rscore:Cscore 0.133072 0.025518 5.215 1.84e-07 ***
Mult(., R, C + L). -4.496485 NA NA NA
Mult(1, ., C + L).R1 0.195482 NA NA NA
Mult(1, ., C + L).R2 -0.296443 NA NA NA
Mult(1, ., C + L).R3 -0.206130 NA NA NA
Mult(1, ., C + L).R4 0.151970 NA NA NA
Mult(1, R, . + L).C1 -0.075450 NA NA NA
Mult(1, R, . + L).C2 0.228093 NA NA NA
Mult(1, R, . + L).C3 0.003945 NA NA NA
Mult(1, R, . + L).C4 -0.230981 NA NA NA
Mult(1, R, C + .).L2 0.414836 NA NA NA
Mult(1, R, C + .).L3 0.104680 NA NA NA
Mult(1, R, C + .).L4 -0.406569 NA NA NA
cov_rcl 0.228615 0.017090 13.377 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Std. Error is NA where coefficient has been constrained or is unidentified
Residual deviance: 55.984 on 44 degrees of freedom
AIC: 390.3
Number of iterations: 52model.summary(model15)# A tibble: 1 × 6
`Model Description` df L2 BIC Delta p
<chr> <int> <dbl> <dbl> <dbl> <dbl>
1 model15 44 56.0 -266. 6.83 0.106mu <- getContrasts(
model15,
pickCoef(model15, "[.]R"),
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit"
)
nu <- getContrasts(
model15,
pickCoef(model15, "[.]C"),
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit"
)
model15.extended <-
freq_tab_3.1 |> gnm(
Freq ~ R + C + L +
Rscore:Cscore +
Mult(1, R, C) +
cov_rcl +
Mult(1, R, L),
constrain = c(13:16, 23:26),
constrainTo = c(mu$qvframe[, 1], mu$qvframe[, 1]),
data = _,
family = poisson,
trace = F,
tolerance = 1e-12
)Initialising
Running start-up iterations..
Running main iterations.......
Donesummary(model15.extended)
Call:
gnm(formula = Freq ~ R + C + L + Rscore:Cscore + Mult(1, R, C) +
cov_rcl + Mult(1, R, L), constrain = c(13:16, 23:26), constrainTo = c(mu$qvframe[,
1], mu$qvframe[, 1]), family = poisson, data = freq_tab_3.1,
tolerance = 1e-12, trace = F)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.98677 -0.75501 0.01756 0.52206 2.97799
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.581787 NA NA NA
R2 -0.994846 NA NA NA
R3 -0.705968 NA NA NA
R4 -1.425526 NA NA NA
C2 -1.647592 0.101523 -16.23 < 2e-16 ***
C3 -2.170836 0.158428 -13.70 < 2e-16 ***
C4 -3.695461 0.241287 -15.32 < 2e-16 ***
L2 -1.415334 0.106117 -13.34 < 2e-16 ***
L3 -2.158852 0.172660 -12.50 < 2e-16 ***
L4 -3.819528 0.272676 -14.01 < 2e-16 ***
Rscore:Cscore 0.133072 0.025394 5.24 1.6e-07 ***
Mult(., R, C). -0.777033 NA NA NA
Mult(1, ., C).R1 0.543679 NA NA NA
Mult(1, ., C).R2 -0.597987 NA NA NA
Mult(1, ., C).R3 -0.388388 NA NA NA
Mult(1, ., C).R4 0.442695 NA NA NA
Mult(1, R, .).C1 -0.115820 NA NA NA
Mult(1, R, .).C2 0.641040 NA NA NA
Mult(1, R, .).C3 0.082145 NA NA NA
Mult(1, R, .).C4 -0.503621 NA NA NA
cov_rcl 0.228615 0.016429 13.91 < 2e-16 ***
Mult(., R, L). -3.469981 NA NA NA
Mult(1, ., L).R1 0.543679 NA NA NA
Mult(1, ., L).R2 -0.597987 NA NA NA
Mult(1, ., L).R3 -0.388388 NA NA NA
Mult(1, ., L).R4 0.442695 NA NA NA
Mult(1, R, .).L1 -0.007184 NA NA NA
Mult(1, R, .).L2 0.224439 NA NA NA
Mult(1, R, .).L3 0.051264 NA NA NA
Mult(1, R, .).L4 -0.234192 NA NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Std. Error is NA where coefficient has been constrained or is unidentified
Residual deviance: 55.984 on 46 degrees of freedom
AIC: 386.3
Number of iterations: 7eta <-
getContrasts(
model15.extended,
pickCoef(model15.extended, "[.]L"),
ref = "mean",
scaleRef = "mean",
scaleWeights = "unit"
)
model15.phis <-
freq_tab_3.1 |> gnm(
Freq ~ R + C + L +
Rscore:Cscore +
Mult(1, R, C) +
cov_rcl +
Mult(1, R, L),
constrain = c(13:20, 23:30),
constrainTo = c(mu$qvframe[, 1],
nu$qvframe[, 1],
mu$qvframe[, 1],
eta$qvframe[, 1]),
data = _,
family = poisson,
trace = F,
tolerance = 1e-12
)Initialising
Running main iterations..
Donesummary(model15.phis)
Call:
gnm(formula = Freq ~ R + C + L + Rscore:Cscore + Mult(1, R, C) +
cov_rcl + Mult(1, R, L), constrain = c(13:20, 23:30), constrainTo = c(mu$qvframe[,
1], nu$qvframe[, 1], mu$qvframe[, 1], eta$qvframe[, 1]),
family = poisson, data = freq_tab_3.1, tolerance = 1e-12, trace = F)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.98677 -0.75501 0.01756 0.52206 2.97799
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.55464 0.08091 43.935 < 2e-16 ***
R2 -0.93784 0.11339 -8.271 < 2e-16 ***
R3 -0.65943 0.12613 -5.228 1.71e-07 ***
R4 -1.42048 0.17754 -8.001 < 2e-16 ***
C2 -1.64759 0.10044 -16.404 < 2e-16 ***
C3 -2.17084 0.15588 -13.927 < 2e-16 ***
C4 -3.69546 0.23861 -15.488 < 2e-16 ***
L2 -1.41533 0.10551 -13.414 < 2e-16 ***
L3 -2.15885 0.17245 -12.519 < 2e-16 ***
L4 -3.81953 0.27221 -14.032 < 2e-16 ***
Rscore:Cscore 0.13307 0.02492 5.341 9.24e-08 ***
Mult(., R, C). -0.64172 0.13829 -4.641 3.48e-06 ***
Mult(1, ., C).R1 0.54368 NA NA NA
Mult(1, ., C).R2 -0.59799 NA NA NA
Mult(1, ., C).R3 -0.38839 NA NA NA
Mult(1, ., C).R4 0.44270 NA NA NA
Mult(1, R, .).C1 -0.17165 NA NA NA
Mult(1, R, .).C2 0.74481 NA NA NA
Mult(1, R, .).C3 0.06806 NA NA NA
Mult(1, R, .).C4 -0.64122 NA NA NA
cov_rcl 0.22862 0.01640 13.939 < 2e-16 ***
Mult(., R, L). -1.13826 0.13726 -8.293 < 2e-16 ***
Mult(1, ., L).R1 0.54368 NA NA NA
Mult(1, ., L).R2 -0.59799 NA NA NA
Mult(1, ., L).R3 -0.38839 NA NA NA
Mult(1, ., L).R4 0.44270 NA NA NA
Mult(1, R, .).L1 -0.04806 NA NA NA
Mult(1, R, .).L2 0.65804 NA NA NA
Mult(1, R, .).L3 0.13012 NA NA NA
Mult(1, R, .).L4 -0.74010 NA NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Std. Error is NA where coefficient has been constrained or is unidentified
Residual deviance: 55.984 on 50 degrees of freedom
AIC: 378.3
Number of iterations: 2mu Estimate Std. Error
Mult(1, ., C + L).R1 0.5436794 0.07818382
Mult(1, ., C + L).R2 -0.5979869 0.06871974
Mult(1, ., C + L).R3 -0.3883879 0.08391754
Mult(1, ., C + L).R4 0.4426954 0.08454001nu Estimate Std. Error
Mult(1, R, . + L).C1 -0.17164747 0.14998137
Mult(1, R, . + L).C2 0.74480785 0.08583798
Mult(1, R, . + L).C3 0.06806188 0.15421719
Mult(1, R, . + L).C4 -0.64122226 0.11322236eta Estimate Std. Error
Mult(1, R, .).L1 -0.04806244 0.08932196
Mult(1, R, .).L2 0.65804171 0.05969167
Mult(1, R, .).L3 0.13011596 0.08427834
Mult(1, R, .).L4 -0.74009523 0.04837401model.summary(model15)# A tibble: 1 × 6
`Model Description` df L2 BIC Delta p
<chr> <int> <dbl> <dbl> <dbl> <dbl>
1 model15 44 56.0 -266. 6.83 0.106bind_rows(model.summary(model1, "Model 1"),
model.summary(model2, "Model 2"),
model.summary(model4, "Model 4"),
model.summary(model6, "Model 6"),
model.summary(model7, "Model 7"),
model.summary(model8, "Model 8"),
model.summary(model9, "Model 9"),
model.summary(model10, "Model 10"),
model.summary(model11, "Model 11"),
model.summary(model12, "Model 12"),
model.summary(model15, "Model 15"))# A tibble: 11 × 6
`Model Description` df L2 BIC Delta p
<chr> <int> <dbl> <dbl> <dbl> <dbl>
1 Model 1 54 544. 149. 23.7 3.44e-82
2 Model 2 39 109. -176. 9.81 1.42e- 8
3 Model 4 35 37.7 -218. 4.63 3.47e- 1
4 Model 6 27 29.0 -169. 4.86 3.61e- 1
5 Model 7 23 21.9 -146. 3.70 5.24e- 1
6 Model 8 36 45.9 -218. 6.25 1.26e- 1
7 Model 9 38 47.1 -231. 6.25 1.47e- 1
8 Model 10 38 50.9 -227. 6.31 7.88e- 2
9 Model 11 38 53.6 -225. 6.58 4.77e- 2
10 Model 12 39 48.6 -237. 6.57 1.40e- 1
11 Model 15 44 56.0 -266. 6.83 1.06e- 1練習問題
問題3.1 Goodman (1979) の表5A，表5B，表5Cの結果を再現しよう．使用するデータはmentalHealthのデータである．表5Aのモデルは，(1) Null Association，(2) Uniform Association，(3) Row-Effect Association，(4) Column-Effect Association，(5) Row and Column Effects (I)， (6) Row and Column Effects (II) の6つである．
- 問題3.1.1：Goodman (1979) の表5Bの分析（全体効果の分解）を再現する．表5Bではまずモデル (1) と (2) の比較から一般的な効果（General Effect）を求めている．次に，モデル (2) と (5) の比較から，行・列効果を求めている． 他の効果はモデル (5) で示され，全体効果はモデル (1) となる．
- 問題3.1.2：Goodman (1979) の表5Cの分析（行と列効果の分割）を再現する．モデル (2) と (3) から行効果が，さらにモデル (3) と (5) から行効果が，そして，モデル (2) と (5) から行・列効果の要素を取り出すことができる．
参考文献
Goodman, Leo A. 1979. “Simple Models for the Analysis of Association in Cross-Classifications Having Ordered Categories.” Journal of the American Statistical Association 74 (367): 537–52.