library(tidyverse)
library(gnm)
library(broom)
library(vcdExtra)先行研究の再現
訳者のGitHubのスクリプトのフォルダ にRのスクリプトファイルを保存しているのでそれも活用してください．
Duncan (1979)
Duncan (1979) の結果を再現する．まずはデータを作成する．
# Data (p. 795)
Freq_duncan_1979 <- c(50, 19, 26, 8, 7, 11, 6, 2,
16, 40, 34, 18, 11, 20, 8, 3,
12, 35, 65, 66, 35, 88, 23, 21,
11, 20, 58, 110, 40, 183, 64, 32,
2, 8, 12, 23, 25, 46, 28, 12,
12, 28, 102, 162, 90, 554, 230, 177,
0, 6, 19, 40, 21, 158, 143, 71,
0, 3, 14, 32, 15, 126, 91, 106
)
d_duncan <- tibble(Freq = Freq_duncan_1979,
O = gl(n = 8, k = 8, length = 8 * 8),
D = gl(n = 8, k = 1, length = 8 * 8))
# データの完成
d_duncan# A tibble: 64 × 3
Freq O D
<dbl> <fct> <fct>
1 50 1 1
2 19 1 2
3 26 1 3
4 8 1 4
5 7 1 5
6 11 1 6
7 6 1 7
8 2 1 8
9 16 2 1
10 40 2 2
# ℹ 54 more rows# Row and column scores
d_duncan <- d_duncan |>
mutate(U = as.numeric(O),
V = as.numeric(D))
# 確認
d_duncan |> count(O, U)# A tibble: 8 × 3
O U n
<fct> <dbl> <int>
1 1 1 8
2 2 2 8
3 3 3 8
4 4 4 8
5 5 5 8
6 6 6 8
7 7 7 8
8 8 8 8d_duncan |> count(D, V)# A tibble: 8 × 3
D V n
<fct> <dbl> <int>
1 1 1 8
2 2 2 8
3 3 3 8
4 4 4 8
5 5 5 8
6 6 6 8
7 7 7 8
8 8 8 8# Diagonal parameters
d_duncan <- d_duncan |>
mutate(Diag = case_when(U == V ~ U,
.default = 0))
# 確認
d_duncan |> xtabs(Diag ~ O + D, data = _) D
O 1 2 3 4 5 6 7 8
1 1 0 0 0 0 0 0 0
2 0 2 0 0 0 0 0 0
3 0 0 3 0 0 0 0 0
4 0 0 0 4 0 0 0 0
5 0 0 0 0 5 0 0 0
6 0 0 0 0 0 6 0 0
7 0 0 0 0 0 0 7 0
8 0 0 0 0 0 0 0 8# factorにする
d_duncan <- d_duncan |> mutate(Diag = factor(Diag))
# weight
d_duncan <- d_duncan |> mutate(W = ifelse(Diag != 0, 0, 1))
# weight変数がどのようになっているのかを確認
d_duncan |> xtabs(W ~ O + D, data = _) D
O 1 2 3 4 5 6 7 8
1 0 1 1 1 1 1 1 1
2 1 0 1 1 1 1 1 1
3 1 1 0 1 1 1 1 1
4 1 1 1 0 1 1 1 1
5 1 1 1 1 0 1 1 1
6 1 1 1 1 1 0 1 1
7 1 1 1 1 1 1 0 1
8 1 1 1 1 1 1 1 0d_duncan |> uncount(weights = W) |>
xtabs(Freq ~ O + D, data = _) D
O 1 2 3 4 5 6 7 8
1 0 19 26 8 7 11 6 2
2 16 0 34 18 11 20 8 3
3 12 35 0 66 35 88 23 21
4 11 20 58 0 40 183 64 32
5 2 8 12 23 0 46 28 12
6 12 28 102 162 90 0 230 177
7 0 6 19 40 21 158 0 71
8 0 3 14 32 15 126 91 0# データの完成
d_duncan# A tibble: 64 × 7
Freq O D U V Diag W
<dbl> <fct> <fct> <dbl> <dbl> <fct> <dbl>
1 50 1 1 1 1 1 0
2 19 1 2 1 2 0 1
3 26 1 3 1 3 0 1
4 8 1 4 1 4 0 1
5 7 1 5 1 5 0 1
6 11 1 6 1 6 0 1
7 6 1 7 1 7 0 1
8 2 1 8 1 8 0 1
9 16 2 1 2 1 0 1
10 40 2 2 2 2 2 0
# ℹ 54 more rowsこのデータについてgnmを用いたモデルの推定を行う．glmでもよい．結果はsummaryでみることができるが，tidyやglanceでも確認でき，これらを使ったほうが表の整理やその後のggplotによる図の作成が簡単である．tidyはgnmをサポートの対象外としているが，glmとほぼ同じ構造なのではじめに警告がでるものの，結果は問題なく出力される．
# (1) Independence
fit_d_1 <- d_duncan |>
gnm(Freq ~ O + D, data = _, family = poisson)
summary(fit_d_1)
Call:
gnm(formula = Freq ~ O + D, family = poisson, data = d_duncan)
Deviance Residuals:
Min 1Q Median 3Q Max
-5.9252 -3.6322 -0.7127 2.0150 12.8584
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.33460 0.13099 10.189 < 2e-16 ***
O2 0.15082 0.12007 1.256 0.209053
O3 0.98373 0.10319 9.533 < 2e-16 ***
O4 1.39016 0.09839 14.130 < 2e-16 ***
O5 0.19004 0.11899 1.597 0.110244
O6 2.35174 0.09213 25.527 < 2e-16 ***
O7 1.26706 0.09966 12.713 < 2e-16 ***
O8 1.09861 0.10165 10.808 < 2e-16 ***
D2 0.43417 0.12642 3.434 0.000594 ***
D3 1.16436 0.11280 10.322 < 2e-16 ***
D4 1.49432 0.10896 13.714 < 2e-16 ***
D5 0.86244 0.11744 7.344 2.08e-13 ***
D6 2.44361 0.10265 23.805 < 2e-16 ***
D7 1.75046 0.10668 16.409 < 2e-16 ***
D8 1.41500 0.10978 12.889 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Residual deviance: 954.49 on 49 degrees of freedom
AIC: 1306.8
Number of iterations: 5tidy(fit_d_1)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.# A tibble: 15 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 1.33 0.131 10.2 2.22e- 24
2 O2 0.151 0.120 1.26 2.09e- 1
3 O3 0.984 0.103 9.53 1.52e- 21
4 O4 1.39 0.0984 14.1 2.49e- 45
5 O5 0.190 0.119 1.60 1.10e- 1
6 O6 2.35 0.0921 25.5 9.84e-144
7 O7 1.27 0.0997 12.7 4.98e- 37
8 O8 1.10 0.102 10.8 3.17e- 27
9 D2 0.434 0.126 3.43 5.94e- 4
10 D3 1.16 0.113 10.3 5.59e- 25
11 D4 1.49 0.109 13.7 8.33e- 43
12 D5 0.862 0.117 7.34 2.08e- 13
13 D6 2.44 0.103 23.8 2.93e-125
14 D7 1.75 0.107 16.4 1.65e- 60
15 D8 1.42 0.110 12.9 5.18e- 38glance(fit_d_1)# 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 -638. 1307. 1339. 954. 49 64# (2) Row effects
fit_d_2 <- d_duncan |>
gnm(Freq ~ O + D + O:V, data = _, family = poisson)
summary(fit_d_2)
Call:
gnm(formula = Freq ~ O + D + O:V, family = poisson, data = d_duncan)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.1378 -1.0958 -0.3250 0.8183 3.3582
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.78811 0.19161 24.989 < 2e-16 ***
O2 -0.69056 0.25946 -2.662 0.00778 **
O3 -1.15391 0.24063 -4.795 1.62e-06 ***
O4 -1.63391 0.24085 -6.784 1.17e-11 ***
O5 -3.24022 0.33365 -9.711 < 2e-16 ***
O6 -1.89930 0.23156 -8.202 < 2e-16 ***
O7 -4.21681 0.30489 -13.831 < 2e-16 ***
O8 -5.13641 0.34152 -15.040 < 2e-16 ***
D2 1.24434 0.13842 8.989 < 2e-16 ***
D3 2.64279 0.15304 17.269 < 2e-16 ***
D4 3.51911 0.18143 19.396 < 2e-16 ***
D5 3.33853 0.21802 15.313 < 2e-16 ***
D6 5.29848 0.24258 21.842 < 2e-16 ***
D7 4.92699 0.27452 17.947 < 2e-16 ***
D8 4.86597 0.30300 16.059 < 2e-16 ***
O1:V -1.32268 0.07546 -17.529 < 2e-16 ***
O2:V -1.05277 0.06472 -16.267 < 2e-16 ***
O3:V -0.73047 0.05137 -14.221 < 2e-16 ***
O4:V -0.54801 0.04808 -11.397 < 2e-16 ***
O5:V -0.47086 0.06225 -7.564 3.92e-14 ***
O6:V -0.32377 0.04351 -7.441 1.00e-13 ***
O7:V -0.11858 0.05167 -2.295 0.02175 *
O8:V 0.00000 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: 130.08 on 42 degrees of freedom
AIC: 496.42
Number of iterations: 4tidy(fit_d_2)# A tibble: 23 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 4.79 0.192 25.0 8.03e-138
2 O2 -0.691 0.259 -2.66 7.78e- 3
3 O3 -1.15 0.241 -4.80 1.62e- 6
4 O4 -1.63 0.241 -6.78 1.17e- 11
5 O5 -3.24 0.334 -9.71 2.69e- 22
6 O6 -1.90 0.232 -8.20 2.36e- 16
7 O7 -4.22 0.305 -13.8 1.67e- 43
8 O8 -5.14 0.342 -15.0 4.03e- 51
9 D2 1.24 0.138 8.99 2.49e- 19
10 D3 2.64 0.153 17.3 8.10e- 67
# ℹ 13 more rowsglance(fit_d_2)# 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 -226. 496. 544. 130. 42 64# (3) Quasi independence, diagonal omitted
fit_d_3 <- d_duncan |> gnm(Freq ~ O + D, data = _, weights = W, family = poisson)
summary(fit_d_3)
Call:
gnm(formula = Freq ~ O + D, family = poisson, data = d_duncan,
weights = W)
Deviance Residuals:
Min 1Q Median 3Q Max
-4.567 -2.318 0.000 1.018 6.152
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.3635 0.1789 2.031 0.0422 *
O2 0.3545 0.1475 2.402 0.0163 *
O3 1.3492 0.1276 10.575 < 2e-16 ***
O4 1.7696 0.1232 14.362 < 2e-16 ***
O5 0.5661 0.1426 3.971 7.16e-05 ***
O6 2.7261 0.1193 22.845 < 2e-16 ***
O7 1.5500 0.1262 12.283 < 2e-16 ***
O8 1.3749 0.1276 10.777 < 2e-16 ***
D2 0.8204 0.1652 4.966 6.83e-07 ***
D3 1.6901 0.1506 11.219 < 2e-16 ***
D4 2.0265 0.1477 13.719 < 2e-16 ***
D5 1.4397 0.1532 9.399 < 2e-16 ***
D6 2.9686 0.1444 20.559 < 2e-16 ***
D7 2.2451 0.1455 15.434 < 2e-16 ***
D8 1.8753 0.1486 12.624 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Residual deviance: 446.84 on 41 degrees of freedom
AIC: 748.82
Number of iterations: 5tidy(fit_d_3)# A tibble: 15 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 0.363 0.179 2.03 4.22e- 2
2 O2 0.354 0.148 2.40 1.63e- 2
3 O3 1.35 0.128 10.6 3.88e- 26
4 O4 1.77 0.123 14.4 9.01e- 47
5 O5 0.566 0.143 3.97 7.16e- 5
6 O6 2.73 0.119 22.8 1.65e-115
7 O7 1.55 0.126 12.3 1.12e- 34
8 O8 1.37 0.128 10.8 4.42e- 27
9 D2 0.820 0.165 4.97 6.83e- 7
10 D3 1.69 0.151 11.2 3.30e- 29
11 D4 2.03 0.148 13.7 7.82e- 43
12 D5 1.44 0.153 9.40 5.50e- 21
13 D6 2.97 0.144 20.6 6.33e- 94
14 D7 2.25 0.145 15.4 9.66e- 54
15 D8 1.88 0.149 12.6 1.56e- 36glance(fit_d_3)# 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 -359. 749. 781. 447. 41 56# (3) Quasi independence, diagonal fitted
fit_d_3b <- d_duncan |> gnm(Freq ~ O + D + Diag, data = _, family = poisson)
summary(fit_d_3b)
Call:
gnm(formula = Freq ~ O + D + Diag, family = poisson, data = d_duncan)
Deviance Residuals:
Min 1Q Median 3Q Max
-4.567 -2.318 0.000 1.018 6.152
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.36346 0.17892 2.031 0.042211 *
O2 0.35445 0.14755 2.402 0.016292 *
O3 1.34923 0.12758 10.575 < 2e-16 ***
O4 1.76957 0.12322 14.362 < 2e-16 ***
O5 0.56615 0.14257 3.971 7.16e-05 ***
O6 2.72612 0.11933 22.845 < 2e-16 ***
O7 1.54999 0.12619 12.283 < 2e-16 ***
O8 1.37488 0.12757 10.777 < 2e-16 ***
D2 0.82044 0.16521 4.966 6.83e-07 ***
D3 1.69007 0.15065 11.219 < 2e-16 ***
D4 2.02648 0.14771 13.719 < 2e-16 ***
D5 1.43966 0.15317 9.399 < 2e-16 ***
D6 2.96856 0.14439 20.559 < 2e-16 ***
D7 2.24512 0.14547 15.434 < 2e-16 ***
D8 1.87534 0.14855 12.624 < 2e-16 ***
Diag1 3.54857 0.22806 15.560 < 2e-16 ***
Diag2 2.15053 0.20737 10.370 < 2e-16 ***
Diag3 0.77163 0.15275 5.052 4.38e-07 ***
Diag4 0.54098 0.12287 4.403 1.07e-05 ***
Diag5 0.84961 0.22965 3.700 0.000216 ***
Diag6 0.25903 0.07531 3.439 0.000583 ***
Diag7 0.80428 0.11438 7.031 2.04e-12 ***
Diag8 1.04976 0.12946 8.109 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Residual deviance: 446.84 on 41 degrees of freedom
AIC: 815.18
Number of iterations: 5tidy(fit_d_3b)# A tibble: 23 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 0.363 0.179 2.03 4.22e- 2
2 O2 0.354 0.148 2.40 1.63e- 2
3 O3 1.35 0.128 10.6 3.88e- 26
4 O4 1.77 0.123 14.4 9.01e- 47
5 O5 0.566 0.143 3.97 7.16e- 5
6 O6 2.73 0.119 22.8 1.65e-115
7 O7 1.55 0.126 12.3 1.12e- 34
8 O8 1.37 0.128 10.8 4.42e- 27
9 D2 0.820 0.165 4.97 6.83e- 7
10 D3 1.69 0.151 11.2 3.30e- 29
# ℹ 13 more rowsglance(fit_d_3b)# 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 -385. 815. 865. 447. 41 64# (4) Uniform association, diagonal omitted
fit_d_4 <- d_duncan |> gnm(Freq ~ O + D + U:V, data = _, weights = W, family = poisson)
summary(fit_d_4)
Call:
gnm(formula = Freq ~ O + D + U:V, family = poisson, data = d_duncan,
weights = W)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.6522 -0.6267 0.0000 0.5913 2.0964
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.109554 0.188406 11.197 < 2e-16 ***
O2 -0.185071 0.150342 -1.231 0.218324
O3 0.229092 0.136939 1.673 0.094337 .
O4 -0.060423 0.147813 -0.409 0.682699
O5 -2.024528 0.187084 -10.822 < 2e-16 ***
O6 -0.590189 0.198922 -2.967 0.003008 **
O7 -2.570773 0.244589 -10.511 < 2e-16 ***
O8 -3.518248 0.285046 -12.343 < 2e-16 ***
D2 0.256818 0.168036 1.528 0.126426
D3 0.581222 0.158565 3.666 0.000247 ***
D4 0.233361 0.167699 1.392 0.164059
D5 -1.099156 0.193262 -5.687 1.29e-08 ***
D6 -0.265555 0.210615 -1.261 0.207362
D7 -1.794390 0.250606 -7.160 8.06e-13 ***
D8 -2.916108 0.289546 -10.071 < 2e-16 ***
U:V 0.136974 0.007489 18.289 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Residual deviance: 58.436 on 40 degrees of freedom
AIC: 362.41
Number of iterations: 4tidy(fit_d_4)# A tibble: 16 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 2.11 0.188 11.2 4.22e-29
2 O2 -0.185 0.150 -1.23 2.18e- 1
3 O3 0.229 0.137 1.67 9.43e- 2
4 O4 -0.0604 0.148 -0.409 6.83e- 1
5 O5 -2.02 0.187 -10.8 2.72e-27
6 O6 -0.590 0.199 -2.97 3.01e- 3
7 O7 -2.57 0.245 -10.5 7.72e-26
8 O8 -3.52 0.285 -12.3 5.33e-35
9 D2 0.257 0.168 1.53 1.26e- 1
10 D3 0.581 0.159 3.67 2.47e- 4
11 D4 0.233 0.168 1.39 1.64e- 1
12 D5 -1.10 0.193 -5.69 1.29e- 8
13 D6 -0.266 0.211 -1.26 2.07e- 1
14 D7 -1.79 0.251 -7.16 8.06e-13
15 D8 -2.92 0.290 -10.1 7.40e-24
16 U:V 0.137 0.00749 18.3 1.01e-74glance(fit_d_4)# 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 -165. 362. 397. 58.4 40 56# log-odds b_hat
fit_d_4$coefficients["U:V"] |> exp() U:V
1.146799 fit_d_4b <- d_duncan |> gnm(Freq ~ O + D + U:V + Diag, data = _, family = poisson)
summary(fit_d_4b)
Call:
gnm(formula = Freq ~ O + D + U:V + Diag, family = poisson, data = d_duncan)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.652e+00 -6.267e-01 -3.932e-08 5.913e-01 2.096e+00
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.109554 0.188406 11.197 < 2e-16 ***
O2 -0.185071 0.150342 -1.231 0.218324
O3 0.229092 0.136939 1.673 0.094337 .
O4 -0.060423 0.147813 -0.409 0.682699
O5 -2.024528 0.187084 -10.822 < 2e-16 ***
O6 -0.590189 0.198922 -2.967 0.003008 **
O7 -2.570773 0.244589 -10.511 < 2e-16 ***
O8 -3.518248 0.285046 -12.343 < 2e-16 ***
D2 0.256818 0.168036 1.528 0.126426
D3 0.581222 0.158565 3.666 0.000247 ***
D4 0.233361 0.167699 1.392 0.164059
D5 -1.099156 0.193262 -5.687 1.29e-08 ***
D6 -0.265555 0.210615 -1.261 0.207362
D7 -1.794390 0.250606 -7.160 8.06e-13 ***
D8 -2.916108 0.289546 -10.071 < 2e-16 ***
U:V 0.136974 0.007489 18.289 < 2e-16 ***
Diag1 1.665495 0.237383 7.016 2.28e-12 ***
Diag2 0.959681 0.212122 4.524 6.06e-06 ***
Diag3 0.021750 0.156530 0.139 0.889490
Diag4 0.226399 0.124364 1.820 0.068689 .
Diag5 0.808646 0.229754 3.520 0.000432 ***
Diag6 0.132277 0.077191 1.714 0.086597 .
Diag7 0.506709 0.115936 4.371 1.24e-05 ***
Diag8 0.221880 0.134803 1.646 0.099771 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Residual deviance: 58.436 on 40 degrees of freedom
AIC: 428.78
Number of iterations: 4tidy(fit_d_4b)# A tibble: 24 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 2.11 0.188 11.2 4.22e-29
2 O2 -0.185 0.150 -1.23 2.18e- 1
3 O3 0.229 0.137 1.67 9.43e- 2
4 O4 -0.0604 0.148 -0.409 6.83e- 1
5 O5 -2.02 0.187 -10.8 2.72e-27
6 O6 -0.590 0.199 -2.97 3.01e- 3
7 O7 -2.57 0.245 -10.5 7.72e-26
8 O8 -3.52 0.285 -12.3 5.33e-35
9 D2 0.257 0.168 1.53 1.26e- 1
10 D3 0.581 0.159 3.67 2.47e- 4
# ℹ 14 more rowsglance(fit_d_4b)# 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 -190. 429. 481. 58.4 40 64# log-odds b_hat
fit_d_4b$coefficients["U:V"] |> exp() U:V
1.146799 # (5) Row effects, diagonal omitted*
fit_d_5 <- d_duncan |> gnm(Freq ~ O + D + O:V, data = _, weights = W, family = poisson)
summary(fit_d_5)
Call:
gnm(formula = Freq ~ O + D + O:V, family = poisson, data = d_duncan,
weights = W)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.5731 -0.5561 0.0000 0.6020 1.4669
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.30716 0.35290 9.371 < 2e-16 ***
O2 0.18160 0.40668 0.447 0.6552
O3 0.07048 0.35788 0.197 0.8439
O4 -0.59546 0.34779 -1.712 0.0869 .
O5 -2.17219 0.42137 -5.155 2.54e-07 ***
O6 -0.68517 0.33031 -2.074 0.0380 *
O7 -2.75820 0.38757 -7.117 1.11e-12 ***
O8 -3.51465 0.43967 -7.994 < 2e-16 ***
D2 1.38525 0.18314 7.564 3.91e-14 ***
D3 2.75348 0.19480 14.135 < 2e-16 ***
D4 3.45305 0.22703 15.210 < 2e-16 ***
D5 3.20664 0.26857 11.940 < 2e-16 ***
D6 5.10309 0.30368 16.804 < 2e-16 ***
D7 4.64707 0.34112 13.623 < 2e-16 ***
D8 4.58416 0.40156 11.416 < 2e-16 ***
O1:V -0.97261 0.09074 -10.719 < 2e-16 ***
O2:V -0.91747 0.07796 -11.768 < 2e-16 ***
O3:V -0.65692 0.06243 -10.522 < 2e-16 ***
O4:V -0.45024 0.05914 -7.613 2.68e-14 ***
O5:V -0.38508 0.07183 -5.361 8.29e-08 ***
O6:V -0.25744 0.05495 -4.685 2.80e-06 ***
O7:V -0.10400 0.06481 -1.605 0.1085
O8:V 0.00000 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.906 on 34 degrees of freedom
AIC: 361.88
Number of iterations: 4tidy(fit_d_5)# A tibble: 23 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 3.31 0.353 9.37 7.16e-21
2 O2 0.182 0.407 0.447 6.55e- 1
3 O3 0.0705 0.358 0.197 8.44e- 1
4 O4 -0.595 0.348 -1.71 8.69e- 2
5 O5 -2.17 0.421 -5.16 2.54e- 7
6 O6 -0.685 0.330 -2.07 3.80e- 2
7 O7 -2.76 0.388 -7.12 1.11e-12
8 O8 -3.51 0.440 -7.99 1.31e-15
9 D2 1.39 0.183 7.56 3.91e-14
10 D3 2.75 0.195 14.1 2.32e-45
# ℹ 13 more rowsglance(fit_d_5)# 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 -159. 362. 409. 45.9 34 56# (5) Row effects, diagonal fitted
fit_d_5b <- d_duncan |> gnm(Freq ~ O + D + O:V + Diag, data = _, family = poisson)
summary(fit_d_5b)
Call:
gnm(formula = Freq ~ O + D + O:V + Diag, family = poisson, data = d_duncan)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.573e+00 -5.561e-01 3.495e-08 6.020e-01 1.467e+00
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.30716 0.35290 9.371 < 2e-16 ***
O2 0.18160 0.40668 0.447 0.655198
O3 0.07048 0.35788 0.197 0.843875
O4 -0.59546 0.34779 -1.712 0.086876 .
O5 -2.17219 0.42137 -5.155 2.54e-07 ***
O6 -0.68517 0.33031 -2.074 0.038049 *
O7 -2.75820 0.38757 -7.117 1.11e-12 ***
O8 -3.51465 0.43967 -7.994 < 2e-16 ***
D2 1.38525 0.18314 7.564 3.91e-14 ***
D3 2.75348 0.19480 14.135 < 2e-16 ***
D4 3.45305 0.22703 15.210 < 2e-16 ***
D5 3.20664 0.26857 11.940 < 2e-16 ***
D6 5.10309 0.30368 16.804 < 2e-16 ***
D7 4.64707 0.34112 13.623 < 2e-16 ***
D8 4.58416 0.40156 11.416 < 2e-16 ***
O1:V -0.97261 0.09074 -10.719 < 2e-16 ***
O2:V -0.91747 0.07796 -11.768 < 2e-16 ***
O3:V -0.65692 0.06243 -10.522 < 2e-16 ***
O4:V -0.45024 0.05914 -7.613 2.68e-14 ***
O5:V -0.38508 0.07183 -5.361 8.29e-08 ***
O6:V -0.25744 0.05495 -4.685 2.80e-06 ***
O7:V -0.10400 0.06481 -1.605 0.108538
O8:V 0.00000 NA NA NA
Diag1 1.57747 0.32240 4.893 9.94e-07 ***
Diag2 0.64980 0.24519 2.650 0.008045 **
Diag3 0.01404 0.17040 0.082 0.934327
Diag4 0.33668 0.13137 2.563 0.010381 *
Diag5 0.80264 0.23091 3.476 0.000509 ***
Diag6 0.13674 0.07784 1.757 0.078973 .
Diag7 0.49482 0.12541 3.946 7.96e-05 ***
Diag8 0.28677 0.17735 1.617 0.105894
---
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.906 on 34 degrees of freedom
AIC: 428.25
Number of iterations: 4tidy(fit_d_5b)# A tibble: 31 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 3.31 0.353 9.37 7.16e-21
2 O2 0.182 0.407 0.447 6.55e- 1
3 O3 0.0705 0.358 0.197 8.44e- 1
4 O4 -0.595 0.348 -1.71 8.69e- 2
5 O5 -2.17 0.421 -5.16 2.54e- 7
6 O6 -0.685 0.330 -2.07 3.80e- 2
7 O7 -2.76 0.388 -7.12 1.11e-12
8 O8 -3.51 0.440 -7.99 1.31e-15
9 D2 1.39 0.183 7.56 3.91e-14
10 D3 2.75 0.195 14.1 2.32e-45
# ℹ 21 more rowsglance(fit_d_5b)# 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 -184. 428. 493. 45.9 34 64# Tabel 1
d_duncan_fitted <-
d_duncan |>
mutate(fitted_1 = predict(fit_d_1, type = "response"),
fitted_2 = predict(fit_d_2, type = "response"),
fitted_3 = predict(fit_d_3, type = "response"),
fitted_4 = predict(fit_d_4, type = "response"),
fitted_5 = predict(fit_d_5, type = "response"))
tab_obs <- d_duncan_fitted |> xtabs(Freq ~ O + D, data = _)
tab_fitted_1 <- d_duncan_fitted |> xtabs(fitted_1 ~ O + D, data = _)
tab_fitted_2 <- d_duncan_fitted |> xtabs(fitted_2 ~ O + D, data = _)
tab_fitted_3 <- d_duncan_fitted |> xtabs(fitted_3 ~ O + D, data = _)
tab_fitted_4 <- d_duncan_fitted |> xtabs(fitted_4 ~ O + D, data = _)
tab_fitted_5 <- d_duncan_fitted |> xtabs(fitted_5 ~ O + D, data = _)
# Observed Counts
tab_obs D
O 1 2 3 4 5 6 7 8
1 50 19 26 8 7 11 6 2
2 16 40 34 18 11 20 8 3
3 12 35 65 66 35 88 23 21
4 11 20 58 110 40 183 64 32
5 2 8 12 23 25 46 28 12
6 12 28 102 162 90 554 230 177
7 0 6 19 40 21 158 143 71
8 0 3 14 32 15 126 91 106# Fitted Counts, Model (4)
tab_fitted_4 |> round(1) D
O 1 2 3 4 5 6 7 8
1 9.5 14.0 22.2 18.0 5.4 14.4 3.6 1.3
2 9.0 15.3 27.9 25.9 9.0 27.2 7.8 3.3
3 15.6 30.5 63.6 67.7 27.0 93.6 30.6 15.0
4 13.4 30.0 71.8 87.7 40.0 159.3 59.7 33.7
5 2.2 5.5 15.2 21.3 11.1 50.8 21.9 14.1
6 10.4 30.6 96.2 154.5 92.7 485.4 239.3 177.3
7 1.6 5.5 20.0 36.9 25.4 152.3 86.2 73.2
8 0.7 2.8 11.7 24.7 19.5 134.4 87.1 84.9# Fitted Counts, Model (5)
tab_fitted_5 |> round(1) D
O 1 2 3 4 5 6 7 8
1 10.3 15.6 23.2 17.6 5.2 13.1 3.1 1.1
2 13.1 20.9 32.8 26.4 8.2 21.9 5.5 2.1
3 15.2 31.5 64.1 66.9 27.1 93.6 30.8 15.0
4 9.6 24.4 61.2 78.6 39.1 166.2 67.2 40.2
5 2.1 5.8 15.4 21.1 11.2 50.8 21.9 14.0
6 10.6 32.9 99.8 155.3 93.8 483.2 236.7 171.8
7 1.6 5.6 19.9 36.1 25.4 152.6 87.2 73.8
8 0.8 3.2 12.8 25.7 20.1 133.7 84.7 79.6d_1 <- (tab_obs/tab_fitted_1) |> diag()
d_2 <- (tab_obs/tab_fitted_2) |> diag()
d_3 <- (tab_obs/tab_fitted_3) |> diag()
d_4 <- (tab_obs/tab_fitted_4) |> diag()
d_5 <- (tab_obs/tab_fitted_5) |> diag()
# Diagonal Ratios
occ <- c("I","II","III","IV","Va","Vb","VI","VII")
tibble(occ, d_1, d_2, d_3, d_4, d_5)# A tibble: 8 × 6
occ d_1 d_2 d_3 d_4 d_5
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 I 13.2 1.56 34.8 5.29 4.84
2 II 5.87 1.57 8.59 2.61 1.92
3 III 2.00 1.09 2.16 1.02 1.01
4 IV 1.62 1.25 1.72 1.25 1.40
5 Va 2.30 1.99 2.34 2.24 2.23
6 Vb 1.21 1.08 1.30 1.14 1.15
7 VI 1.84 1.34 2.24 1.66 1.64
8 VII 2.26 1.16 2.86 1.25 1.33# Table 2
Model <- c("(1) Independence",
"(2) Row effects",
"(3) Quasi independence, diagonal omitted",
"(4) Uniform association, diagonal omitted",
"(5) Row effects, diagonal omitted")
# 複数のモデルの適合度をglanceで示してまとめる
list(fit_d_1, fit_d_2, fit_d_3, fit_d_4, fit_d_5) |>
map_dfr(glance) |>
mutate(Model, .before = 1) |>
dplyr::select(Model, L2 = deviance, df = df.residual)# A tibble: 5 × 3
Model L2 df
<chr> <dbl> <int>
1 (1) Independence 954. 49
2 (2) Row effects 130. 42
3 (3) Quasi independence, diagonal omitted 447. 41
4 (4) Uniform association, diagonal omitted 58.4 40
5 (5) Row effects, diagonal omitted 45.9 34# Figure 2
b_4 <- exp(fit_d_4$coefficients["U:V"] * (7:0))
b_5 <- fit_d_5$coefficients[grep(":V", names(fit_d_5$coefficients))]
b_5["O8:V"] <- 0
b_5 <- exp(b_5 * (-1))
d_duncan_b_4 <- tibble(y = b_4, Model = "Model 4", x = gl(8, 1))
d_duncan_b_5 <- tibble(y = b_5, Model = "Model 5", x = gl(8, 1))
d <- bind_rows(d_duncan_b_4, d_duncan_b_5)
d |>
mutate(x = factor(x, levels = 8:1)) |>
ggplot(aes(x = x, y = y,
group = Model,
linetype = Model,
color = Model)) +
geom_line() +
ylim(1,3) +
scale_linetype_manual(values = c("longdash","solid")) +
labs(x = "FATHER'S OCCUPATION (i)", y = "b_i/b_8") +
theme_classic() +
ggthemes::scale_color_colorblind()
# mosaic plot
mosaic(fit_d_1, ~ O + D, residuals_type = "rstandard")
mosaic(fit_d_2, ~ O + D, residuals_type = "rstandard")
mosaic(fit_d_3b, ~ O + D, residuals_type = "rstandard")
mosaic(fit_d_4b, ~ O + D, residuals_type = "rstandard")
mosaic(fit_d_5b, ~ O + D, residuals_type = "rstandard")
Yamaguchi (1987)
Yamaguchi (1987) の結果を再現する．なおデータはvcdExtraパッケージのYamaguchi87からも入手できる．
# Data
Freq_yamaguchi_1987 <- c(1275, 364, 274, 272, 17,
1055, 597, 394, 443, 31,
1043, 587, 1045, 951, 47,
1159, 791, 1323, 2046, 52,
666, 496, 1031, 1632, 646,
474, 129, 87, 124, 11,
300, 218, 171, 220, 8,
438, 254, 669, 703, 16,
601, 388, 932, 1789, 37,
76, 56, 125, 295, 191,
127, 101, 24, 30, 12,
86, 207, 64, 61, 13,
43, 73, 122, 60, 13,
35, 51, 62, 66, 11,
109, 206, 184, 253, 325)
d_yamaguchi <- tibble(O = gl(5, k = 5, length = 3 * 5 * 5),
D = gl(5, k = 1, length = 3 * 5 * 5),
C = gl(3, k = 5 * 5, length = 3 * 5 * 5),
Freq = Freq_yamaguchi_1987)
d_yamaguchi# A tibble: 75 × 4
O D C Freq
<fct> <fct> <fct> <dbl>
1 1 1 1 1275
2 1 2 1 364
3 1 3 1 274
4 1 4 1 272
5 1 5 1 17
6 2 1 1 1055
7 2 2 1 597
8 2 3 1 394
9 2 4 1 443
10 2 5 1 31
# ℹ 65 more rows# Diagonal parameters
d_yamaguchi <- d_yamaguchi |>
mutate(Diag = case_when(O == D ~ O, .default = "0") |> factor())
# Row and column scores
d_yamaguchi <- d_yamaguchi |>
mutate(U = as.numeric(O),
V = as.numeric(D))
# Country dummy
d_yamaguchi <- d_yamaguchi |>
mutate(C1 = case_match(C,
c("2","3") ~ 0,
"1" ~ 1),
C2 = case_match(C,
c("1","3") ~ 0,
"2" ~ 1),
C3 = case_match(C,
c("1","2") ~ 0,
"3" ~ 1))
# Off diagonal cells (Diag*C)
# No association between R and C, given L
fit_y_1 <- d_yamaguchi |> glm(Freq ~ O*C + D*C + Diag*C, data = _, family = poisson)
glance(fit_y_1)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -929. 1943. 2040. 1336. 33 75# U0-Cross-nationally homogeneous uniform association
fit_y_2 <- d_yamaguchi |> glm(Freq ~ O*C + D*C + Diag*C + U:V, data = _, family = poisson)
glance(fit_y_2)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -349. 784. 884. 175. 32 75b_2 <- fit_y_2$coefficients[grep("U:V", names(fit_y_2$coefficients))]
b_2 U:V
0.1850182 # U1-Cross-nationally uniform uniform association
fit_y_3 <- d_yamaguchi |> glm(Freq ~ O*C + D*C + Diag*C + U:V*C, data = _, family = poisson)
glance(fit_y_3)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -344. 779. 883. 166. 30 75b_3 <- fit_y_3$coefficients[grep("U:V", names(fit_y_3$coefficients))]
b_3 U:V C2:U:V C3:U:V
0.18852859 0.01069127 -0.04990850 # U2-Cross-nationally uniform uniform association
fit_y_4 <- d_yamaguchi |>
glm(Freq ~ O*C + D*C + Diag*C + U:V + U:V:C3, data = _, family = poisson)
glance(fit_y_4)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -345. 777. 879. 167. 31 75b_4 <- fit_y_4$coefficients[grep("U:V", names(fit_y_4$coefficients))]
b_4 U:V U:V:C3
0.19052983 -0.05190974 # R0-Cross-nationally homogeneous row effect association
fit_y_5 <- d_yamaguchi |>
glm(Freq ~ O*C + D*C + Diag*C + O:V, data = _, family = poisson)
glance(fit_y_5)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -339. 771. 877. 156. 29 75# R1-Cross-nationally uniform row effect association
fit_y_6 <- d_yamaguchi |>
glm(Freq ~ O*C + D*C + Diag*C + O:V + U:V*C, data = _, family = poisson)
glance(fit_y_6)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -335. 766. 878. 148. 27 75b_6 <- fit_y_6$coefficients[grep("V:U", names(fit_y_6$coefficients))]
b_6 V:U C2:V:U C3:V:U
NA 0.00369588 -0.05053513 # R2-Cross-nationally uniform row effect association
fit_y_7 <- d_yamaguchi |>
glm(Freq ~ O*C + D*C + Diag*C + O:V + U:V:C3, data = _, family = poisson)
glance(fit_y_7)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -335. 764. 873. 148. 28 75b_7 <- fit_y_7$coefficients[grep("V:U", names(fit_y_7$coefficients))]
b_7 V:U:C3
-0.05122437 # C0-Cross-nationally homogeneous column effect association
fit_y_8 <- d_yamaguchi |>
glm(Freq ~ O*C + D*C + Diag*C + U:D, data = _, family = poisson)
glance(fit_y_8)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -295. 682. 789. 67.7 29 75# C1-Cross-nationally uniform column effect association
fit_y_9 <- d_yamaguchi |>
glm(Freq ~ O*C + D*C + Diag*C + U:D + U:V*C, data = _, family = poisson)
glance(fit_y_9)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -291. 679. 790. 60.3 27 75b_9 <- fit_y_9$coefficients[grep("U:V", names(fit_y_9$coefficients))]
b_9 U:V C2:U:V C3:U:V
NA 0.01491735 -0.04242365 # C2-Cross-nationally uniform column effect association
fit_y_10 <- d_yamaguchi |>
glm(Freq ~ O*C + D*C + Diag*C + U:D + U:V:C3, data = _, family = poisson)
glance(fit_y_10)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -292. 678. 787. 61.3 28 75b_10 <- fit_y_10$coefficients[grep("U:V", names(fit_y_10$coefficients))]
b_10 U:V:C3
-0.04530956 # H0-Cross-nationally homogeneous homogeneous row & column effect association
fit_y_11 <- d_yamaguchi |>
gnm(Freq ~ O*C + D*C + Diag*C + MultHomog(O, D),
data = _, family = poisson)Initialising
Running start-up iterations..
Running main iterations.....
Doneglance(fit_y_11)# 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 -315. 721. 828. 107. 29 75# H1-Cross-nationally uniform homogeneous row & column effect association
fit_y_12 <- d_yamaguchi |>
gnm(Freq ~ O*C + D*C + Diag*C + MultHomog(O, D) + U:V:C2 + U:V:C3,
data = _, family = poisson)Initialising
Running start-up iterations..
Running main iterations.....
Doneglance(fit_y_12)# 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 -310. 717. 828. 98.2 27 75b_12 <- fit_y_12$coefficients[grep("U:V", names(fit_y_12$coefficients))]
b_12 U:V:C2 U:V:C3
-0.002548059 -0.052199354 # H2-Cross-nationally uniform homogeneous row & column effect association
fit_y_13 <- d_yamaguchi |> gnm(Freq ~ O*C + D*C + Diag*C + MultHomog(O, D) + U:V:C3,
data = _, family = poisson)Initialising
Running start-up iterations..
Running main iterations.....
Doneglance(fit_y_13)# 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 -310. 715. 824. 98.2 28 75b_13 <- fit_y_13$coefficients[grep("U:V", names(fit_y_13$coefficients))]
b_13 U:V:C3
-0.05173675 # R+C0 -Cross-nationally uniform row & column effect association
fit_y_14 <- d_yamaguchi |>
glm(Freq ~ O*C + D*C + Diag*C + U:D + O:V, data = _, family = poisson)
glance(fit_y_14)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -281. 659. 773. 38.8 26 75# R+C1 -Cross-nationally uniform row & column effect association
fit_y_15 <- d_yamaguchi |>
glm(Freq ~ O*C + D*C + Diag*C + U:D + O:V + U:V:C2 + U:V:C3, data = _, family = poisson)
glance(fit_y_15)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -278. 658. 776. 33.3 24 75b_15 <- fit_y_15$coefficients[grep("U:V", names(fit_y_15$coefficients))]
b_15 U:V:C2 U:V:C3
0.003287747 -0.041213885 # R+C2-Cross-nationally uniform row & column effect association
fit_y_16 <- d_yamaguchi |>
glm(Freq ~ O*C + D*C + Diag*C + U:D + O:V + U:V:C3, data = _, family = poisson)
glance(fit_y_16)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -278. 656. 772. 33.3 25 75b_16 <- fit_y_16$coefficients[grep("U:V", names(fit_y_16$coefficients))]
b_16 U:V:C3
-0.0418251 d_yamaguchi <- d_yamaguchi |>
mutate(SYM = c(rep(c(0, 1, 2, 3, 4,
1, 0, 5, 6, 7,
2, 5, 0, 8, 9,
3, 6, 8, 0, 10,
4, 7, 9, 10, 0), 3)) |> factor())
d_yamaguchi |> dplyr::select(O, D, C, SYM) |> unique()# A tibble: 75 × 4
O D C SYM
<fct> <fct> <fct> <fct>
1 1 1 1 0
2 1 2 1 1
3 1 3 1 2
4 1 4 1 3
5 1 5 1 4
6 2 1 1 1
7 2 2 1 0
8 2 3 1 5
9 2 4 1 6
10 2 5 1 7
# ℹ 65 more rows# QS0-Cross-nationally homogeneous quasi-symmetry
fit_y_17 <- d_yamaguchi |>
glm(Freq ~ O*C + D*C + Diag*C + SYM, data = _, family = poisson)
glance(fit_y_17)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -315. 723. 832. 107. 28 75tidy(fit_y_17) |> print(n = Inf)# A tibble: 52 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 4.47 0.0870 51.4 0
2 O2 0.688 0.0658 10.4 1.59e- 25
3 O3 1.64 0.0602 27.2 1.06e-162
4 O4 2.34 0.0575 40.7 0
5 O5 2.09 0.0650 32.1 1.13e-226
6 C2 -1.23 0.0742 -16.6 6.96e- 62
7 C3 -2.41 0.115 -21.0 3.48e- 98
8 D2 -0.304 0.0519 -5.86 4.74e- 9
9 D3 0.359 0.0447 8.04 9.13e- 16
10 D4 0.811 0.0414 19.6 2.04e- 85
11 D5 -2.40 0.0947 -25.3 1.06e-141
12 Diag1 2.68 0.0914 29.3 3.33e-189
13 Diag2 1.54 0.0992 15.5 3.40e- 54
14 Diag3 0.487 0.0535 9.10 9.21e- 20
15 Diag4 -0.000521 0.0479 -0.0109 9.91e- 1
16 Diag5 2.31 0.0964 24.0 6.95e-127
17 SYM1 1.78 0.0722 24.7 1.90e-134
18 SYM2 0.811 0.0517 15.7 1.73e- 55
19 SYM3 0.256 0.0532 4.80 1.56e- 6
20 SYM4 NA NA NA NA
21 SYM5 0.541 0.0553 9.78 1.35e- 22
22 SYM6 0.132 0.0563 2.34 1.92e- 2
23 SYM7 NA NA NA NA
24 SYM8 NA NA NA NA
25 SYM9 NA NA NA NA
26 SYM10 NA NA NA NA
27 O2:C2 0.0124 0.0798 0.155 8.77e- 1
28 O3:C2 0.373 0.0734 5.09 3.57e- 7
29 O4:C2 0.574 0.0714 8.03 9.65e- 16
30 O5:C2 -1.02 0.0788 -12.9 2.39e- 38
31 O2:C3 0.0805 0.120 0.673 5.01e- 1
32 O3:C3 -0.688 0.116 -5.92 3.20e- 9
33 O4:C3 -1.19 0.120 -9.92 3.45e- 23
34 O5:C3 0.409 0.0967 4.23 2.29e- 5
35 C2:D2 0.0286 0.0541 0.529 5.97e- 1
36 C3:D2 1.05 0.0901 11.6 2.65e- 31
37 C2:D3 0.285 0.0489 5.83 5.64e- 9
38 C3:D3 0.336 0.0893 3.77 1.65e- 4
39 C2:D4 0.516 0.0498 10.4 3.12e- 25
40 C3:D4 0.215 0.0862 2.50 1.25e- 2
41 C2:D5 0.0726 0.148 0.491 6.23e- 1
42 C3:D5 2.05 0.179 11.4 4.08e- 30
43 C2:Diag1 0.243 0.0917 2.65 8.12e- 3
44 C3:Diag1 0.106 0.148 0.720 4.72e- 1
45 C2:Diag2 0.184 0.102 1.80 7.26e- 2
46 C3:Diag2 0.225 0.128 1.76 7.91e- 2
47 C2:Diag3 0.128 0.0743 1.72 8.60e- 2
48 C3:Diag3 0.617 0.140 4.41 1.05e- 5
49 C2:Diag4 0.00793 0.0624 0.127 8.99e- 1
50 C3:Diag4 -0.0490 0.163 -0.300 7.64e- 1
51 C2:Diag5 0.961 0.174 5.53 3.26e- 8
52 C3:Diag5 -0.730 0.189 -3.86 1.12e- 4# QS1-Cross-nationally uniform quasi-symmetry
fit_y_18 <- d_yamaguchi |>
glm(Freq ~ O*C + D*C + Diag*C + SYM + U:V:C2 + U:V:C3, data = _, family = poisson)
glance(fit_y_18)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -310. 719. 832. 98.2 26 75b_18 <- fit_y_18$coefficients[grep("U:V", names(fit_y_18$coefficients))]
b_18 U:V:C2 U:V:C3
-0.002538939 -0.052099153 # QS2-Cross-nationally uniform quasi-symmetry
fit_y_19 <- d_yamaguchi |>
glm(Freq ~ O*C + D*C + Diag*C + SYM + U:V:C3, data = _, family = poisson)
glance(fit_y_19)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -310. 717. 828. 98.2 27 75b_19 <- fit_y_19$coefficients[grep("U:V", names(fit_y_19$coefficients))]
b_19 U:V:C3
-0.05163812 # FI0-Cross-nationallyhomogeneous full two-way interaction
fit_y_20 <- d_yamaguchi |>
glm(Freq ~ O*C + D*C + Diag*C + O:D, data = _, family = poisson)
glance(fit_y_20)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -279. 665. 788. 36.2 22 75# FI1-Cross-nationally uniform full two-way interaction
fit_y_21 <- d_yamaguchi |>
glm(Freq ~ O*C + D*C + Diag*C + O:D + U:V:C2 + U:V:C3, data = _, family = poisson)
glance(fit_y_21)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -277. 663. 791. 30.7 20 75b_21 <- fit_y_21$coefficients[grep("U:V", names(fit_y_21$coefficients))]
b_21 U:V:C2 U:V:C3
0.003476173 -0.041081646 # FI2-Cross-nationally uniform full two-way interaction
fit_y_22 <- d_yamaguchi |>
glm(Freq ~ O*C + D*C + Diag*C + O:D + U:V:C3, data = _, family = poisson)
glance(fit_y_22)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -277. 661. 787. 30.8 21 75b_22 <- fit_y_22$coefficients[grep("U:V", names(fit_y_22$coefficients))]
b_22 U:V:C3
-0.04172632 Model <- c("No association between R and C, given L",
"U0-Cross-nationally homogeneous uniform association",
"U1-Cross-nationally uniform uniform association",
"U2-Cross-nationally uniform uniform association",
"R0-Cross-nationally homogeneous row effect association",
"R1-Cross-nationally uniform row effect association",
"R2-Cross-nationally uniform row effect association",
"C0-Cross-nationally homogeneous column effect association",
"C1-Cross-nationally uniform column effect association",
"C2-Cross-nationally uniform column effect association",
"H0-Cross-nationally homogeneous homogeneous row & column effect association",
"H1-Cross-nationally uniform homogeneousrow & column effect association",
"H2-Cross-nationally uniform homogeneousrow & column effect association",
"R+C0 -Cross-nationally uniform row & column effect association",
"R+C1 -Cross-nationally uniform row & column effect association",
"R+C2-Cross-nationally uniform row & column effect association",
"QS0-Cross-nationally homogeneous quasi-symmetry",
"QS1-Cross-nationally uniform quasi-symmetry",
"QS2-Cross-nationally uniform quasi-symmetry",
"FI0-Cross-nationallyhomogeneous full two-way interaction",
"FI1-Cross-nationally uniform full two-way interaction",
"FI2-Cross-nationally uniform full two-way interaction")
list(fit_y_1, fit_y_2, fit_y_3, fit_y_4, fit_y_5, fit_y_6, fit_y_7, fit_y_8, fit_y_9, fit_y_10,
fit_y_11, fit_y_12, fit_y_13, fit_y_14, fit_y_15, fit_y_16, fit_y_17, fit_y_18, fit_y_19,
fit_y_20, fit_y_21, fit_y_22) |>
map_dfr(glance) |>
mutate(Model = Model, .before = 1) |>
dplyr::select(-c(2:3)) |>
print(n = Inf)# A tibble: 22 × 7
Model logLik AIC BIC deviance df.residual nobs
<chr> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 No association between R and C… -929. 1943. 2040. 1336. 33 75
2 U0-Cross-nationally homogeneou… -349. 784. 884. 175. 32 75
3 U1-Cross-nationally uniform un… -344. 779. 883. 166. 30 75
4 U2-Cross-nationally uniform un… -345. 777. 879. 167. 31 75
5 R0-Cross-nationally homogeneou… -339. 771. 877. 156. 29 75
6 R1-Cross-nationally uniform ro… -335. 766. 878. 148. 27 75
7 R2-Cross-nationally uniform ro… -335. 764. 873. 148. 28 75
8 C0-Cross-nationally homogeneou… -295. 682. 789. 67.7 29 75
9 C1-Cross-nationally uniform co… -291. 679. 790. 60.3 27 75
10 C2-Cross-nationally uniform co… -292. 678. 787. 61.3 28 75
11 H0-Cross-nationally homogeneou… -315. 721. 828. 107. 29 75
12 H1-Cross-nationally uniform ho… -310. 717. 828. 98.2 27 75
13 H2-Cross-nationally uniform ho… -310. 715. 824. 98.2 28 75
14 R+C0 -Cross-nationally uniform… -281. 659. 773. 38.8 26 75
15 R+C1 -Cross-nationally uniform… -278. 658. 776. 33.3 24 75
16 R+C2-Cross-nationally uniform … -278. 656. 772. 33.3 25 75
17 QS0-Cross-nationally homogeneo… -315. 723. 832. 107. 28 75
18 QS1-Cross-nationally uniform q… -310. 719. 832. 98.2 26 75
19 QS2-Cross-nationally uniform q… -310. 717. 828. 98.2 27 75
20 FI0-Cross-nationallyhomogeneou… -279. 665. 788. 36.2 22 75
21 FI1-Cross-nationally uniform f… -277. 663. 791. 30.7 20 75
22 FI2-Cross-nationally uniform f… -277. 661. 787. 30.8 21 75H1とH2が微妙に異なるが誤差の範囲か．また，FI1の\(\beta_1\)については符号が Yamaguchi (1987) のものとは異なっているがこちらが正しいか．
# R+C2: No uniform difference between any two
fit_y_23 <- d_yamaguchi |>
gnm(Freq ~ O*C + D*C + Diag*C + U:D + O:V + U:V:C3,
data = _, family = poisson)
glance(fit_y_23)# 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 -278. 656. 772. 33.3 25 75tidy(fit_y_23) |> print(n = Inf)# A tibble: 53 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 7.07 0.177 39.9 0
2 O2 -0.197 0.253 -0.778 4.37e- 1
3 O3 -0.536 0.376 -1.43 1.54e- 1
4 O4 -0.722 0.549 -1.32 1.88e- 1
5 O5 -1.43 0.656 -2.18 2.96e- 2
6 C2 -1.24 0.0742 -16.6 3.30e-62
7 C3 -2.60 0.130 -19.9 1.61e-88
8 O2:C2 0.0105 0.0799 0.131 8.96e- 1
9 O3:C2 0.375 0.0734 5.12 3.13e- 7
10 O4:C2 0.576 0.0716 8.04 8.77e-16
11 O5:C2 -1.03 0.0790 -13.1 4.41e-39
12 O2:C3 0.244 0.131 1.87 6.19e- 2
13 O3:C3 -0.405 0.154 -2.63 8.56e- 3
14 O4:C3 -0.777 0.189 -4.11 4.00e- 5
15 O5:C3 0.911 0.219 4.17 3.08e- 5
16 D2 -0.670 0.161 -4.16 3.18e- 5
17 D3 -0.922 0.240 -3.84 1.25e- 4
18 D4 -0.625 0.351 -1.78 7.49e- 2
19 D5 -2.26 0.545 -4.15 3.31e- 5
20 C2:D2 0.0189 0.0536 0.354 7.24e- 1
21 C3:D2 1.21 0.109 11.1 7.17e-29
22 C2:D3 0.285 0.0491 5.80 6.73e- 9
23 C3:D3 0.640 0.156 4.10 4.20e- 5
24 C2:D4 0.539 0.0498 10.8 2.77e-27
25 C3:D4 0.672 0.216 3.11 1.87e- 3
26 C2:D5 0.157 0.148 1.06 2.90e- 1
27 C3:D5 2.31 0.275 8.38 5.22e-17
28 Diag1 0.247 0.135 1.83 6.76e- 2
29 Diag2 0.249 0.0715 3.49 4.90e- 4
30 Diag3 0.361 0.0460 7.85 4.25e-15
31 Diag4 -0.00158 0.0564 -0.0280 9.78e- 1
32 Diag5 3.09 0.177 17.5 1.15e-68
33 C2:Diag1 0.246 0.0917 2.69 7.22e- 3
34 C3:Diag1 0.332 0.168 1.98 4.80e- 2
35 C2:Diag2 0.199 0.102 1.95 5.13e- 2
36 C3:Diag2 0.246 0.128 1.93 5.38e- 2
37 C2:Diag3 0.130 0.0742 1.75 8.06e- 2
38 C3:Diag3 0.590 0.140 4.22 2.49e- 5
39 C2:Diag4 -0.0134 0.0625 -0.214 8.31e- 1
40 C3:Diag4 -0.0629 0.163 -0.387 6.99e- 1
41 C2:Diag5 0.893 0.174 5.13 2.89e- 7
42 C3:Diag5 -0.264 0.219 -1.21 2.28e- 1
43 D1:U 0.177 0.136 1.30 1.93e- 1
44 D2:U 0.245 0.118 2.08 3.73e- 2
45 D3:U 0.442 0.0941 4.70 2.60e- 6
46 D4:U 0.476 0.0784 6.07 1.28e- 9
47 D5:U 0 NA NA NA
48 O1:V -0.342 0.105 -3.27 1.08e- 3
49 O2:V -0.275 0.0678 -4.05 5.05e- 5
50 O3:V -0.116 0.0440 -2.63 8.55e- 3
51 O4:V 0 NA NA NA
52 O5:V 0 NA NA NA
53 U:V:C3 -0.0418 0.0178 -2.35 1.88e- 2# R+C3: Uniform difference between G.B. and U.S.Xie (1992)
Xie (1992) の結果を再現する．
# Data
d_xie <- d_yamaguchi
# 確認
d_xie# A tibble: 75 × 11
O D C Freq Diag U V C1 C2 C3 SYM
<fct> <fct> <fct> <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <fct>
1 1 1 1 1275 1 1 1 1 0 0 0
2 1 2 1 364 0 1 2 1 0 0 1
3 1 3 1 274 0 1 3 1 0 0 2
4 1 4 1 272 0 1 4 1 0 0 3
5 1 5 1 17 0 1 5 1 0 0 4
6 2 1 1 1055 0 2 1 1 0 0 1
7 2 2 1 597 2 2 2 1 0 0 0
8 2 3 1 394 0 2 3 1 0 0 5
9 2 4 1 443 0 2 4 1 0 0 6
10 2 5 1 31 0 2 5 1 0 0 7
# ℹ 65 more rows# Off diagonal cells (Diag*C)
# Null association between R and C, given L
fit_x_1 <- d_xie |> glm(Freq ~ O*C + D*C + Diag*C, data = _, family = poisson)
glance(fit_x_1)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -929. 1943. 2040. 1336. 33 75# R0 Cross-nationally homogeneous row effect association
fit_x_2 <- d_xie |> glm(Freq ~ O*C + D*C + Diag*C + O:V, data = _, family = poisson)
glance(fit_x_2)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -339. 771. 877. 156. 29 75# R1 Cross-nationally uniform row effect association
fit_x_3 <- d_xie |> glm(Freq ~ O*C + D*C + Diag*C + O:V + U:V*C, data = _, family = poisson)
glance(fit_x_3)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -335. 766. 878. 148. 27 75# Rx Cross-nationally log-multiplicative row effect association
fit_x_4 <- d_xie |> gnm(Freq ~ O*C + D*C + Diag*C + Mult(C,O,V), data = _, family = poisson)Initialising
Running start-up iterations..
Running main iterations.......
Doneglance(fit_x_4)# 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 -335. 766. 877. 148. 27 75# C0-Cross-nationally homogeneous column effect association
fit_x_5 <- d_xie |> glm(Freq ~ O*C + D*C + Diag*C + U:D, data = _, family = poisson)
glance(fit_x_5)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -295. 682. 789. 67.7 29 75# C1-Cross-nationally uniform column effect association
fit_x_6 <- d_xie |> glm(Freq ~ O*C + D*C + Diag*C + U:D + U:V*C, data = _, family = poisson)
glance(fit_x_6)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 34313. 74 -291. 679. 790. 60.3 27 75# Cx-Cross-nationally log-multiplicative column effect association
fit_x_7 <- d_xie |> gnm(Freq ~ O*C + D*C + Diag*C + Mult(C,U,D), data = _, family = poisson)Initialising
Running start-up iterations..
Running main iterations.......
Doneglance(fit_x_7)# 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 -291. 677. 789. 58.8 27 75# (R+C) 0 - Cross-nationally homogeneous row and column effects association I
fit_x_8 <- d_xie |> gnm(Freq ~ O*C + D*C + Diag*C + O:V + U:D, data = _, family = poisson)
glance(fit_x_8)# 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 -281. 659. 773. 38.8 26 75# (R+C) u - Cross-nationally uniform row and column effects association I
fit_x_9 <- d_xie |> gnm(Freq ~ O*C + D*C + Diag*C + O:V + U:D + U:V*C, data = _, family = poisson)
glance(fit_x_9)# 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 -278. 658. 776. 33.3 24 75# (R+C) x - Cross-nationally log-multiplicative row and column effects association I
fit_x_10 <- d_xie |> gnm(Freq ~ O*C + D*C + Diag*C + Mult(C,O:V+U:D), data = _, family = poisson)Initialising
Running start-up iterations..
Running main iterations........
Doneglance(fit_x_10)# 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 -278. 658. 776. 33.0 24 75# RC0 Cross-nationally homogeneous row and column effects association
fit_x_11 <- d_xie |> gnm(Freq ~ O*C + D*C + Diag*C + Mult(O,D), data = _, family = poisson)Initialising
Running start-up iterations..
Running main iterations.......
Doneglance(fit_x_11)# 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 -280. 658. 772. 37.7 26 75# RCx Cross-nationally log-multiplicative row and column effects association
fit_x_12 <- d_xie |> gnm(Freq ~ O*C + D*C + Diag*C + Mult(Exp(C),O,D), data = _, family = poisson)Initialising
Running start-up iterations..
Running main iterations........
Doneglance(fit_x_12)# 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 -277. 657. 775. 32.1 24 75# FI0 Cross-nationally homogeneous full two-way R and C interaction
fit_x_13 <- d_xie |> gnm(Freq ~ O*C + D*C + Diag*C + Mult(1,O*D), data = _, family = poisson)Initialising
Running start-up iterations..
Running main iterations....
Doneglance(fit_x_13)# 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 -279. 665. 788. 36.2 22 75# FIu Cross-nationally uniform full two-way R and C interaction
fit_x_14 <- d_xie |> gnm(Freq ~ O*C + D*C + Diag*C + Mult(1,O*D) + U:V*C, data = _, family = poisson)Initialising
Running start-up iterations..
Running main iterations.....
Doneglance(fit_x_14)# 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 -277. 663. 791. 30.7 20 75# FIx Cross-nationally log-multiplicative full two-way R and C interaction
fit_x_15 <- d_xie |> gnm(Freq ~ O*C + D*C + Diag*C + Mult(Exp(C),O*D), data = _, family = poisson)Initialising
Running start-up iterations..
Running main iterations........
Doneglance(fit_x_15)# 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 -277. 664. 791. 30.9 20 75Model <- c("Null association between R and C, given L",
"R0 Cross-nationally homogeneous row effect association",
"R1 Cross-nationally uniform row effect association",
"Rx Cross-nationally log-multiplicative row effect association",
"C0-Cross-nationally homogeneous column effect association",
"C1-Cross-nationally uniform column effect association",
"Cx-Cross-nationally log-multiplicative column effect association",
"(R+C) 0 - Cross-nationally homogeneous row and column effects association I",
"(R+C) u - Cross-nationally uniform row and column effects association I",
"(R+C) x - Cross-nationally log-multiplicative row and column effects association I",
"RC0 Cross-nationally homogeneous row and column effects association",
"RCx Cross-nationally log-multiplicative row and column effects association",
"FI0 Cross-nationally homogeneous full two-way R and C interaction",
"FIu Cross-nationally uniform full two-way R and C interaction",
"FIx Cross-nationally log-multiplicative full two-way R and C interaction")
list(fit_x_1, fit_x_2, fit_x_3, fit_x_4, fit_x_5,
fit_x_6, fit_x_7, fit_x_8, fit_x_9, fit_x_10,
fit_x_11, fit_x_12, fit_x_13, fit_x_14, fit_x_15) |>
map_dfr(glance) |>
mutate(Model = Model, .before = 1) |>
dplyr::select(-c(2:3)) |>
print(n = Inf)# A tibble: 15 × 7
Model logLik AIC BIC deviance df.residual nobs
<chr> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 Null association between R and… -929. 1943. 2040. 1336. 33 75
2 R0 Cross-nationally homogeneou… -339. 771. 877. 156. 29 75
3 R1 Cross-nationally uniform ro… -335. 766. 878. 148. 27 75
4 Rx Cross-nationally log-multip… -335. 766. 877. 148. 27 75
5 C0-Cross-nationally homogeneou… -295. 682. 789. 67.7 29 75
6 C1-Cross-nationally uniform co… -291. 679. 790. 60.3 27 75
7 Cx-Cross-nationally log-multip… -291. 677. 789. 58.8 27 75
8 (R+C) 0 - Cross-nationally hom… -281. 659. 773. 38.8 26 75
9 (R+C) u - Cross-nationally uni… -278. 658. 776. 33.3 24 75
10 (R+C) x - Cross-nationally log… -278. 658. 776. 33.0 24 75
11 RC0 Cross-nationally homogeneo… -280. 658. 772. 37.7 26 75
12 RCx Cross-nationally log-multi… -277. 657. 775. 32.1 24 75
13 FI0 Cross-nationally homogeneo… -279. 665. 788. 36.2 22 75
14 FIu Cross-nationally uniform f… -277. 663. 791. 30.7 20 75
15 FIx Cross-nationally log-multi… -277. 664. 791. 30.9 20 75参考文献
Duncan, Otis Dudley. 1979. “How Destination Depends on Origin in the Occupational Mobility Table.” The American Journal of Sociology 84 (4): 793–803. https://www.jstor.org/stable/2778024.
Xie, Yu. 1992. “The Log-Multiplicative Layer Effect Model for Comparing Mobility Tables.” American Sociological Review 57 (3): 380–95. https://doi.org/10.2307/2096242.
Yamaguchi, Kazuo. 1987. “Models for Comparing Mobility Tables: Toward Parsimony and Substance.” American Sociological Review 52 (4): 482–94. https://doi.org/10.2307/2095293.