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.44889506 1
polviews2 polviews2 1.11173082 2
polviews3 polviews3 1.07062463 3
polviews4 polviews4 2.38620551 4
polviews5 polviews5 1.37103186 5
polviews6 polviews6 1.30913698 6
polviews7 polviews7 -0.07836634 7
fefam2 fefam2 0.98038359 8
fefam3 fefam3 0.43295515 9
fefam4 fefam4 -0.87285503 10
natfare2 natfare2 0.35356150 11
natfare3 natfare3 0.29974988 12
Mult(., polviews, fefam). Mult(., polviews, fefam). 2.72741175 13
Mult(1, ., fefam).polviews1 Mult(1, ., fefam).polviews1 0.21153027 14
Mult(1, ., fefam).polviews2 Mult(1, ., fefam).polviews2 0.57366789 15
Mult(1, ., fefam).polviews3 Mult(1, ., fefam).polviews3 0.45026306 16
Mult(1, ., fefam).polviews4 Mult(1, ., fefam).polviews4 -0.02697483 17
Mult(1, ., fefam).polviews5 Mult(1, ., fefam).polviews5 0.04038332 18
Mult(1, ., fefam).polviews6 Mult(1, ., fefam).polviews6 -0.61521139 19
Mult(1, ., fefam).polviews7 Mult(1, ., fefam).polviews7 -0.94507437 20
Mult(1, polviews, .).fefam1 Mult(1, polviews, .).fefam1 0.33829591 21
Mult(1, polviews, .).fefam2 Mult(1, polviews, .).fefam2 0.07300986 22
Mult(1, polviews, .).fefam3 Mult(1, polviews, .).fefam3 -0.19982036 23
Mult(1, polviews, .).fefam4 Mult(1, polviews, .).fefam4 -0.36614113 24
Mult(., polviews, natfare). Mult(., polviews, natfare). 3.04182815 25
Mult(1, ., natfare).polviews1 Mult(1, ., natfare).polviews1 -1.17242040 26
Mult(1, ., natfare).polviews2 Mult(1, ., natfare).polviews2 -0.56839055 27
Mult(1, ., natfare).polviews3 Mult(1, ., natfare).polviews3 -0.35469807 28
Mult(1, ., natfare).polviews4 Mult(1, ., natfare).polviews4 0.05233282 29
Mult(1, ., natfare).polviews5 Mult(1, ., natfare).polviews5 -0.15794806 30
Mult(1, ., natfare).polviews6 Mult(1, ., natfare).polviews6 0.86877820 31
Mult(1, ., natfare).polviews7 Mult(1, ., natfare).polviews7 0.87257574 32
Mult(1, polviews, .).natfare1 Mult(1, polviews, .).natfare1 -0.19244649 33
Mult(1, polviews, .).natfare2 Mult(1, polviews, .).natfare2 -0.06881662 34
Mult(1, polviews, .).natfare3 Mult(1, polviews, .).natfare3 0.19771269 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)
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.53383944 1
polviews2 1.05910145 2
polviews3 1.00850170 3
polviews4 2.28433681 4
polviews5 1.28035534 5
polviews6 1.15025088 6
polviews7 -0.24284087 7
fefam2 0.99266224 8
fefam3 0.45226874 9
fefam4 -0.84193457 10
natfare2 0.34921745 11
natfare3 0.30640340 12
Mult(., polviews, fefam + natfare). -1.05643325 13
Mult(1, ., fefam + natfare).polviews1 -1.36835469 14
Mult(1, ., fefam + natfare).polviews2 -1.28039686 15
Mult(1, ., fefam + natfare).polviews3 -0.90102278 16
Mult(1, ., fefam + natfare).polviews4 0.06170142 17
Mult(1, ., fefam + natfare).polviews5 -0.20795228 18
Mult(1, ., fefam + natfare).polviews6 1.57835547 19
Mult(1, ., fefam + natfare).polviews7 1.92456047 20
Mult(1, polviews, . + natfare).fefam1 0.61653553 21
Mult(1, polviews, . + natfare).fefam2 0.30072953 22
Mult(1, polviews, . + natfare).fefam3 0.03156322 23
Mult(1, polviews, . + natfare).fefam4 -0.13260419 24
Mult(1, polviews, fefam + .).natfare2 -0.15612825 25
Mult(1, polviews, fefam + .).natfare3 -0.57283168 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 = 151.608771
Running start-up iterations
Start-up iteration 1. Deviance = 74.992568
Start-up iteration 2. Deviance = 72.948285
Running main iterations
Iteration 1. Deviance = 72.773956
Iteration 2. Deviance = 72.773834
Iteration 3. 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.497 NA NA NA
2 polviews2 1.05 0.184 5.70 1.20e- 8
3 polviews3 1.01 NA NA NA
4 polviews4 2.32 NA NA NA
5 polviews5 1.32 NA NA NA
6 polviews6 1.25 NA NA NA
7 polviews7 -0.131 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.5843655 NA NA NA
R2 -0.9844527 NA NA NA
R3 -0.7522372 NA NA NA
R4 -1.5738349 NA NA NA
C2 -1.6097190 NA NA NA
C3 -2.0783249 NA NA NA
C4 -3.4330190 NA NA NA
L2 -1.4738150 NA NA NA
L3 -2.2418273 NA NA NA
L4 -3.8020714 NA NA NA
Rscore:Cscore 0.1296440 0.0272792 4.752 2.01e-06 ***
Mult(., R, C). 0.3755915 NA NA NA
Mult(1, ., C).R1 0.1452548 NA NA NA
Mult(1, ., C).R2 -0.3244586 NA NA NA
Mult(1, ., C).R3 -0.1448572 NA NA NA
Mult(1, ., C).R4 0.1841136 NA NA NA
Mult(1, R, .).C1 0.5249563 NA NA NA
Mult(1, R, .).C2 -3.3891552 NA NA NA
Mult(1, R, .).C3 -0.2655787 NA NA NA
Mult(1, R, .).C4 2.5859285 NA NA NA
Rscore:Lscore 0.2544964 0.0280285 9.080 < 2e-16 ***
Mult(., R, L). 3.2515305 NA NA NA
Mult(1, ., L).R1 0.2882521 NA NA NA
Mult(1, ., L).R2 -0.3021096 NA NA NA
Mult(1, ., L).R3 -0.2384094 NA NA NA
Mult(1, ., L).R4 0.1714509 NA NA NA
Mult(1, R, .).L1 0.0455793 NA NA NA
Mult(1, R, .).L2 -0.4419901 NA NA NA
Mult(1, R, .).L3 -0.0811452 NA NA NA
Mult(1, R, .).L4 0.5224424 NA NA NA
Cscore:Lscore 0.2017249 0.0258549 7.802 < 2e-16 ***
Mult(., C, L). -0.8754712 NA NA NA
Mult(1, ., L).C1 0.3002275 NA NA NA
Mult(1, ., L).C2 -0.0412263 NA NA NA
Mult(1, ., L).C3 -0.9734008 NA NA NA
Mult(1, ., L).C4 0.3744697 NA NA NA
Mult(1, C, .).L1 -0.0237909 NA NA NA
Mult(1, C, .).L2 -0.1650570 NA NA NA
Mult(1, C, .).L3 0.2543660 NA NA NA
Mult(1, C, .).L4 -0.0006375 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: 71pickCoef(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.5843654594 1
R2 R2 -0.9844527452 2
R3 R3 -0.7522372217 3
R4 R4 -1.5738348698 4
C2 C2 -1.6097190368 5
C3 C3 -2.0783248718 6
C4 C4 -3.4330189812 7
L2 L2 -1.4738149842 8
L3 L3 -2.2418272853 9
L4 L4 -3.8020714474 10
Rscore:Cscore Rscore:Cscore 0.1296439973 11
Mult(., R, C). Mult(., R, C). 0.3755915111 12
Mult(1, ., C).R1 Mult(1, ., C).R1 0.1452548156 13
Mult(1, ., C).R2 Mult(1, ., C).R2 -0.3244585766 14
Mult(1, ., C).R3 Mult(1, ., C).R3 -0.1448572169 15
Mult(1, ., C).R4 Mult(1, ., C).R4 0.1841135671 16
Mult(1, R, .).C1 Mult(1, R, .).C1 0.5249563024 17
Mult(1, R, .).C2 Mult(1, R, .).C2 -3.3891551966 18
Mult(1, R, .).C3 Mult(1, R, .).C3 -0.2655787314 19
Mult(1, R, .).C4 Mult(1, R, .).C4 2.5859285336 20
Rscore:Lscore Rscore:Lscore 0.2544963888 21
Mult(., R, L). Mult(., R, L). 3.2515305041 22
Mult(1, ., L).R1 Mult(1, ., L).R1 0.2882521361 23
Mult(1, ., L).R2 Mult(1, ., L).R2 -0.3021096268 24
Mult(1, ., L).R3 Mult(1, ., L).R3 -0.2384094176 25
Mult(1, ., L).R4 Mult(1, ., L).R4 0.1714508503 26
Mult(1, R, .).L1 Mult(1, R, .).L1 0.0455792508 27
Mult(1, R, .).L2 Mult(1, R, .).L2 -0.4419900743 28
Mult(1, R, .).L3 Mult(1, R, .).L3 -0.0811452356 29
Mult(1, R, .).L4 Mult(1, R, .).L4 0.5224423967 30
Cscore:Lscore Cscore:Lscore 0.2017248677 31
Mult(., C, L). Mult(., C, L). -0.8754712159 32
Mult(1, ., L).C1 Mult(1, ., L).C1 0.3002274500 33
Mult(1, ., L).C2 Mult(1, ., L).C2 -0.0412263278 34
Mult(1, ., L).C3 Mult(1, ., L).C3 -0.9734007547 35
Mult(1, ., L).C4 Mult(1, ., L).C4 0.3744696897 36
Mult(1, C, .).L1 Mult(1, C, .).L1 -0.0237909080 37
Mult(1, C, .).L2 Mult(1, C, .).L2 -0.1650570164 38
Mult(1, C, .).L3 Mult(1, C, .).L3 0.2543660164 39
Mult(1, C, .).L4 Mult(1, C, .).L4 -0.0006375429 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.........................................................
................................................................................
................................................................................
................................................................................
................................................................................
................................................................................
...........................................
Donesummary(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.578416 NA NA NA
R2 -0.944753 NA NA NA
R3 -0.738628 NA NA NA
R4 -1.588272 NA NA NA
C2 -1.595823 NA NA NA
C3 -2.098590 NA NA NA
C4 -3.442358 NA NA NA
L2 -1.486233 NA NA NA
L3 -2.173240 NA NA NA
L4 -3.816157 NA NA NA
Rscore:Cscore 0.129644 0.027279 4.752 2.01e-06 ***
Mult(., R, C). 0.259905 NA NA NA
Mult(1, ., C).R1 -1.608369 NA NA NA
Mult(1, ., C).R2 3.138381 NA NA NA
Mult(1, ., C).R3 1.323396 NA NA NA
Mult(1, ., C).R4 -2.001061 NA NA NA
Mult(1, R, .).C1 -0.139469 NA NA NA
Mult(1, R, .).C2 0.420250 NA NA NA
Mult(1, R, .).C3 -0.026422 NA NA NA
Mult(1, R, .).C4 -0.434189 NA NA NA
Rscore:Lscore 0.254496 0.028029 9.080 < 2e-16 ***
Mult(., R, L). -1.396314 NA NA NA
Mult(1, ., L).R1 -0.559839 NA NA NA
Mult(1, ., L).R2 0.539134 NA NA NA
Mult(1, ., L).R3 0.420554 NA NA NA
Mult(1, ., L).R4 -0.342411 NA NA NA
Mult(1, R, .).L1 0.031112 NA NA NA
Mult(1, R, .).L2 -0.578807 NA NA NA
Mult(1, R, .).L3 -0.127413 NA NA NA
Mult(1, R, .).L4 0.627638 NA NA NA
Cscore:Lscore 0.201725 0.025855 7.802 < 2e-16 ***
Mult(., C, L). -0.860457 NA NA NA
Mult(1, ., L).C1 -0.362025 NA NA NA
Mult(1, ., L).C2 -0.141578 NA NA NA
Mult(1, ., L).C3 0.460249 NA NA NA
Mult(1, ., L).C4 -0.409957 NA NA NA
Mult(1, C, .).L1 0.003017 NA NA NA
Mult(1, C, .).L2 0.225644 NA NA NA
Mult(1, C, .).L3 -0.435340 NA NA NA
Mult(1, C, .).L4 -0.033471 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.57990 NA NA NA
R2 -0.94951 NA NA NA
R3 -0.72054 NA NA NA
R4 -1.55484 NA NA NA
C2 -1.62289 NA NA NA
C3 -2.02529 NA NA NA
C4 -3.45997 NA NA NA
L2 -1.42960 NA NA NA
L3 -2.28435 NA NA NA
L4 -3.75862 NA NA NA
Rscore:Cscore 0.13353 0.02535 5.266 1.39e-07 ***
Mult(., R, C + L). 2.92403 NA NA NA
Mult(1, ., C + L).R1 -0.20409 NA NA NA
Mult(1, ., C + L).R2 0.27631 NA NA NA
Mult(1, ., C + L).R3 0.19241 NA NA NA
Mult(1, ., C + L).R4 -0.14679 NA NA NA
Mult(1, R, . + L).C1 -0.15289 NA NA NA
Mult(1, R, . + L).C2 0.33409 NA NA NA
Mult(1, R, . + L).C3 -0.02460 NA NA NA
Mult(1, R, . + L).C4 -0.42106 NA NA NA
Mult(1, R, C + .).L2 0.66532 NA NA NA
Mult(1, R, C + .).L3 0.16531 NA NA NA
Mult(1, R, C + .).L4 -0.64176 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). -2.02786 NA NA NA
Mult(1, ., L).C1 -0.00891 NA NA NA
Mult(1, ., L).C2 -0.09883 NA NA NA
Mult(1, ., L).C3 -0.43722 NA NA NA
Mult(1, ., L).C4 0.01948 NA NA NA
Mult(1, C, .).L1 -0.11190 NA NA NA
Mult(1, C, .).L2 -0.31604 NA NA NA
Mult(1, C, .).L3 0.22080 NA NA NA
Mult(1, C, .).L4 -0.10226 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: 92mu <- 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.58684 NA NA NA
R2 -1.00837 NA NA NA
R3 -0.76912 NA NA NA
R4 -1.56186 NA NA NA
C2 -1.58868 NA NA NA
C3 -2.05112 NA NA NA
C4 -3.48062 NA NA NA
L2 -1.39505 NA NA NA
L3 -2.23301 NA NA NA
L4 -3.81283 NA NA NA
Rscore:Cscore 0.13353 0.02520 5.300 1.16e-07 ***
Mult(., R, C). 0.18009 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.51254 NA NA NA
Mult(1, R, .).C2 2.77718 NA NA NA
Mult(1, R, .).C3 0.35407 NA NA NA
Mult(1, R, .).C4 -2.32416 NA NA NA
Rscore:Lscore 0.24848 0.02521 9.856 < 2e-16 ***
Mult(., R, L). -0.14415 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.29632 NA NA NA
Mult(1, R, .).L2 -5.31844 NA NA NA
Mult(1, R, .).L3 -1.09875 NA NA NA
Mult(1, R, .).L4 5.71226 NA NA NA
Cscore:Lscore 0.20270 0.02560 7.919 < 2e-16 ***
Mult(., C, L). -0.56592 NA NA NA
Mult(1, ., L).C1 0.12534 NA NA NA
Mult(1, ., L).C2 -0.11925 NA NA NA
Mult(1, ., L).C3 -1.03970 NA NA NA
Mult(1, ., L).C4 0.20256 NA NA NA
Mult(1, C, .).L1 -0.09148 NA NA NA
Mult(1, C, .).L2 -0.36040 NA NA NA
Mult(1, C, .).L3 0.34681 NA NA NA
Mult(1, C, .).L4 -0.07878 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: 72nu1 <- 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.59464 NA NA NA
R2 -0.99884 NA NA NA
R3 -0.69917 NA NA NA
R4 -1.39870 NA NA NA
C2 -1.63854 NA NA NA
C3 -2.18983 NA NA NA
C4 -3.65220 NA NA NA
L2 -1.45492 NA NA NA
L3 -2.06941 NA NA NA
L4 -3.75963 NA NA NA
Rscore:Cscore 0.13336 0.02542 5.246 1.55e-07 ***
Mult(., R, C). 0.72034 NA NA NA
Mult(1, ., C).R1 -0.53827 NA NA NA
Mult(1, ., C).R2 0.59767 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.14477 NA NA NA
Mult(1, R, .).C2 0.67497 NA NA NA
Mult(1, R, .).C3 0.08782 NA NA NA
Mult(1, R, .).C4 -0.55348 NA NA NA
cov_rcl 0.22423 0.01709 13.123 < 2e-16 ***
Mult(., R, L). 1.15302 NA NA NA
Mult(1, ., L).R1 -0.53827 NA NA NA
Mult(1, ., L).R2 0.59767 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.00175 NA NA NA
Mult(1, R, .).L2 0.69740 NA NA NA
Mult(1, R, .).L3 0.16390 NA NA NA
Mult(1, R, .).L4 -0.68683 NA NA NA
Mult(., C, L). -0.51691 NA NA NA
Mult(1, ., L).C1 -0.74057 NA NA NA
Mult(1, ., L).C2 -0.56691 NA NA NA
Mult(1, ., L).C3 0.67430 NA NA NA
Mult(1, ., L).C4 -0.93238 NA NA NA
Mult(1, C, .).L1 -0.01977 NA NA NA
Mult(1, C, .).L2 0.26471 NA NA NA
Mult(1, C, .).L3 -0.33064 NA NA NA
Mult(1, C, .).L4 0.03623 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: 90nu1 <- 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: 8mu 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.58914 NA NA NA
R2 -0.98648 NA NA NA
R3 -0.69914 NA NA NA
R4 -1.42479 NA NA NA
C2 -1.70214 NA NA NA
C3 -2.18510 NA NA NA
C4 -3.66751 NA NA NA
L2 -1.48988 NA NA NA
L3 -2.17766 NA NA NA
L4 -3.74647 NA NA NA
Rscore:Cscore 0.13307 0.02552 5.215 1.84e-07 ***
Mult(., R, C + L). -0.52349 NA NA NA
Mult(1, ., C + L).R1 -1.46764 NA NA NA
Mult(1, ., C + L).R2 2.24813 NA NA NA
Mult(1, ., C + L).R3 1.56595 NA NA NA
Mult(1, ., C + L).R4 -1.13897 NA NA NA
Mult(1, R, . + L).C1 0.07175 NA NA NA
Mult(1, R, . + L).C2 -0.27342 NA NA NA
Mult(1, R, . + L).C3 -0.01853 NA NA NA
Mult(1, R, . + L).C4 0.24861 NA NA NA
Mult(1, R, C + .).L2 -0.47172 NA NA NA
Mult(1, R, C + .).L3 -0.11903 NA NA NA
Mult(1, R, C + .).L4 0.46232 NA NA NA
cov_rcl 0.22862 0.01709 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: 48model.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.58597 NA NA NA
R2 -1.00362 NA NA NA
R3 -0.71313 NA NA NA
R4 -1.42630 NA NA NA
C2 -1.64759 0.10152 -16.23 < 2e-16 ***
C3 -2.17084 0.15843 -13.70 < 2e-16 ***
C4 -3.69546 0.24129 -15.32 < 2e-16 ***
L2 -1.41533 0.10612 -13.34 < 2e-16 ***
L3 -2.15885 0.17266 -12.50 < 2e-16 ***
L4 -3.81953 0.27268 -14.01 < 2e-16 ***
Rscore:Cscore 0.13307 0.02539 5.24 1.6e-07 ***
Mult(., R, C). -1.10096 NA NA NA
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.08925 NA NA NA
Mult(1, R, .).C2 -0.44493 NA NA NA
Mult(1, R, .).C3 -0.05047 NA NA NA
Mult(1, R, .).C4 0.36295 NA NA NA
cov_rcl 0.22862 0.01643 13.91 < 2e-16 ***
Mult(., R, L). 0.34715 NA NA NA
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.02587 NA NA NA
Mult(1, R, .).L2 2.28937 NA NA NA
Mult(1, R, .).L3 0.55836 NA NA NA
Mult(1, R, .).L4 -2.29498 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: 4eta <-
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.