4  第4章

第4章では3元表の条件付き連関モデルを扱っている．

library(tidyverse)
library(magrittr)
library(gnm)
library(broom)
library(knitr)
library(logmult)

モデル適合度を表示するための関数を準備する．すべて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)     #p<-ifelse(p<0.001,"<0.001",p)
  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)
}

4.1 表4.3

# Table 4.3
Freq <- c(201,  29,   8,  13,   5,     152,  29,   2,   8,   0,
           18,   6,   3,   6,   0,      17,  12,   0,   3,   0,
          109,  74, 164,  89,  16,     101, 336,   9, 134,   2,
            7,   6,  45,  30,   6,       7,  41,   7,  63,   0,
          247,  58,  20,  23,   2,     288,  51,   1,  17,   3,
           48,  11,  16,  13,   1,      47,  38,   2,  18,   0,
          157,  68, 178, 116,  27,     165, 321,  27, 168,   1,
            7,   7,  50,  42,   5,      12,  25,   5,  29,   6)

Educ <- gl(4, 10, 4 * 5 * 2 * 2)
Occ <- gl(5, 1, 4 * 5 * 2 * 2)
Sex <- gl(2, 5, 4 * 5 * 2 * 2)
Year <- gl(2, 40, 4 * 5 * 2 * 2)
L <- Year:Sex
levels(L) <- factor(1:4)

Rscore <- as.numeric(Educ)
Cscore <- as.numeric(Occ)

freq_tab_3.1 <- tibble(Freq, Educ, Occ, Sex, Year, L, Rscore, Cscore) 
freq_tab_3.1

# A tibble: 80 × 8
    Freq Educ  Occ   Sex   Year  L     Rscore Cscore
   <dbl> <fct> <fct> <fct> <fct> <fct>  <dbl>  <dbl>
 1   201 1     1     1     1     1          1      1
 2    29 1     2     1     1     1          1      2
 3     8 1     3     1     1     1          1      3
 4    13 1     4     1     1     1          1      4
 5     5 1     5     1     1     1          1      5
 6   152 1     1     2     1     2          1      1
 7    29 1     2     2     1     2          1      2
 8     2 1     3     2     1     2          1      3
 9     8 1     4     2     1     2          1      4
10     0 1     5     2     1     2          1      5
# ℹ 70 more rows

# Model 1 - RC(0)-L(homogeneous)
model1 <- freq_tab_3.1 |> 
  gnm(
    Freq ~ Educ + Occ + L + Educ:L + Occ:L,
    family = poisson,
    data = _,
    trace = F,
    tolerance = 1e-12
  )
summary(model1)

Call:
gnm(formula = Freq ~ Educ + Occ + L + Educ:L + Occ:L, family = poisson, 
    data = freq_tab_3.1, tolerance = 1e-12, trace = F)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-9.2080  -2.2181  -0.4393   1.7534  10.3493  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  4.63188    0.07546  61.385  < 2e-16 ***
Educ2       -2.04867    0.18496 -11.076  < 2e-16 ***
Educ3        0.56850    0.07822   7.268 3.65e-13 ***
Educ4       -1.00188    0.12060  -8.307  < 2e-16 ***
Occ2        -1.06920    0.10808  -9.893  < 2e-16 ***
Occ3        -0.42050    0.08678  -4.846 1.26e-06 ***
Occ4        -0.88688    0.10115  -8.768  < 2e-16 ***
Occ5        -2.51829    0.20006 -12.588  < 2e-16 ***
L2          -0.58321    0.11600  -5.028 4.96e-07 ***
L3           0.35568    0.09908   3.590 0.000331 ***
L4           0.38267    0.09802   3.904 9.46e-05 ***
Educ2:L2     0.26213    0.26589   0.986 0.324189    
Educ3:L2     0.54569    0.11433   4.773 1.82e-06 ***
Educ4:L2     0.52029    0.16809   3.095 0.001966 ** 
Educ2:L3     0.67937    0.21978   3.091 0.001994 ** 
Educ3:L3    -0.12382    0.10396  -1.191 0.233635    
Educ4:L3    -0.14652    0.16251  -0.902 0.367276    
Educ2:L4     0.81653    0.21566   3.786 0.000153 ***
Educ3:L4     0.07042    0.10180   0.692 0.489078    
Educ4:L4    -0.54042    0.17410  -3.104 0.001908 ** 
Occ2:L2      1.48066    0.13298  11.135  < 2e-16 ***
Occ3:L2     -2.31314    0.25826  -8.957  < 2e-16 ***
Occ4:L2      0.60040    0.13656   4.397 1.10e-05 ***
Occ5:L2     -2.41258    0.73731  -3.272 0.001067 ** 
Occ2:L3     -0.09004    0.14424  -0.624 0.532465    
Occ3:L3     -0.13260    0.11618  -1.141 0.253725    
Occ4:L3      0.02568    0.13253   0.194 0.846332    
Occ5:L3     -0.05541    0.26603  -0.208 0.835011    
Occ2:L4      0.90622    0.12622   7.179 7.00e-13 ***
Occ3:L4     -2.26247    0.19508 -11.598  < 2e-16 ***
Occ4:L4      0.09529    0.12843   0.742 0.458122    
Occ5:L4     -1.41745    0.37680  -3.762 0.000169 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

Residual deviance: 1370.9 on 48 degrees of freedom
AIC: 1796.7

Number of iterations: 6

model.summary(model1)

# A tibble: 1 × 6
  `Model Description`    df    L2   BIC Delta         p
  <chr>               <int> <dbl> <dbl> <dbl>     <dbl>
1 model1                 48 1371.  972.  23.8 1.37e-255

# Model 2 - RC(1)-L(homogeneous)
model2 <- freq_tab_3.1 |> 
  gnm(
  Freq ~ Educ + Occ + L + Educ:L + Occ:L + Mult(1, Educ, Occ),
  family = poisson,
  data = _,
  trace = F,
  tolerance = 1e-12
)

Initialising
Running start-up iterations..
Running main iterations..................................
Done

model.summary(model2)

# A tibble: 1 × 6
  `Model Description`    df    L2   BIC Delta        p
  <chr>               <int> <dbl> <dbl> <dbl>    <dbl>
1 model2                 42  144. -205.  6.39 4.54e-13

# Model 3 - RC(1)-L(heterogeneous)
model3 <- freq_tab_3.1 |> 
  gnm(
  Freq ~ Educ + Occ + L + Educ:L + Occ:L + Mult(L, L:Educ, L:Occ),
  family = poisson,
  data = _,
  trace = F,
  tolerance = 1e-12
)

Initialising
Running start-up iterations..
Running main iterations.........................................................
........................
Done

model.summary(model3)

# A tibble: 1 × 6
  `Model Description`    df    L2   BIC Delta          p
  <chr>               <int> <dbl> <dbl> <dbl>      <dbl>
1 model3                 24  69.1 -130.  3.04 0.00000303

# Model 4 - RC(2)-L(homogeneous)
model4 <- freq_tab_3.1 |> 
  gnm(
    Freq ~ Educ + Occ + L + Educ:L + Occ:L + instances(Mult(1, Educ, Occ), 2),
    family = poisson,
    data = _, 
    trace = F,
    tolerance = 1e-8
  )

Initialising
Running start-up iterations..
Running main iterations................
Done

model.summary(model4)

# A tibble: 1 × 6
  `Model Description`    df    L2   BIC Delta        p
  <chr>               <int> <dbl> <dbl> <dbl>    <dbl>
1 model4                 38  117. -199.  5.24 4.93e-10

# Model 4 - RC(2)-L(heterogeneous)
model5 <- freq_tab_3.1 |> 
  gnm(
    Freq ~ Educ + Occ + L + Educ:L + Occ:L + 
      instances(Mult(L, L:Educ, L:Occ), 2),
    family = poisson,
    data = _,
    trace = F,
    tolerance = 1e-10,
    iterStart = 5,
    iterMax = 1e6,
    verbose = F
  )
model.summary(model5)

# A tibble: 1 × 6
  `Model Description`    df    L2   BIC Delta     p
  <chr>               <int> <dbl> <dbl> <dbl> <dbl>
1 model5                  8  5.83 -60.7 0.437 0.666

# Model 6 - RC(3)-L(homogeneous)
model6 <- freq_tab_3.1 |> 
  gnm(
    Freq ~ Educ + Occ + L + Educ:L + Occ:L + 
      instances(Mult(1, Educ, Occ), 3),
    family = poisson,
    data = _,
    trace = F,
    tolerance = 1e-12
  )

Initialising
Running start-up iterations..
Running main iterations........
Done

model.summary(model6)

# A tibble: 1 × 6
  `Model Description`    df    L2   BIC Delta        p
  <chr>               <int> <dbl> <dbl> <dbl>    <dbl>
1 model6                 36  113. -186.  5.19 6.60e-10

# Model 7 - RC(3)-L(heterogeneous)
model7 <- freq_tab_3.1 |> 
  gnm(
    Freq ~ Educ + Occ + L + Educ:L + Occ:L + 
      instances(Mult(L, L:Educ, L:Occ), 3),
    family = poisson,
    data = _,
    trace = F,
    tolerance = 1e-12
  )

Initialising
Running start-up iterations..
Running main iterations.........................................................
................................................................................
................................................................................
................................................................................
................................................................................
................................................................................
...........................................
Done

Warning in gnmFit(modelTools = modelTools, y = y, constrain = constrain, : Fitting algorithm has either not converged or converged
to a non-solution of the likelihood equations.
Use exitInfo() for numerical details of last iteration.

model.summary(model7)

# A tibble: 1 × 6
  `Model Description`    df       L2      BIC         Delta     p
  <chr>               <int>    <dbl>    <dbl>         <dbl> <dbl>
1 model7                  0 4.32e-11 4.32e-11 0.00000000657     0

# Model 8- U+RC (homogeneous)
model8 <- freq_tab_3.1 |> 
  gnm(
    Freq ~ Educ + Occ + L + Educ:L + Occ:L + 
      Rscore:Cscore + 
      Mult(1, Educ, Occ),
    family = poisson,
    data = _,
    trace = F,
    tolerance = 1e-6,
    iterStart = 20,
    iterMax = 100000
  )

Initialising
Running start-up iterations....................
Running main iterations..............
Done

summary(model8)

Call:
gnm(formula = Freq ~ Educ + Occ + L + Educ:L + Occ:L + Rscore:Cscore + 
    Mult(1, Educ, Occ), family = poisson, data = freq_tab_3.1, 
    tolerance = 1e-06, iterStart = 20, iterMax = 1e+05, trace = F)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-3.09104  -0.98816  -0.05028   0.65984   4.48366  

Coefficients:
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)            4.39711         NA      NA       NA    
Educ2                 -3.14283         NA      NA       NA    
Educ3                 -2.04278         NA      NA       NA    
Educ4                 -5.76600         NA      NA       NA    
Occ2                  -2.85536         NA      NA       NA    
Occ3                  -4.44772         NA      NA       NA    
Occ4                  -6.38812         NA      NA       NA    
Occ5                  -9.68116         NA      NA       NA    
L2                    -0.45256    0.10138  -4.464 8.04e-06 ***
L3                     0.30784    0.08546   3.602 0.000315 ***
L4                     0.27495    0.08559   3.212 0.001317 ** 
Educ2:L2               0.21438    0.26790   0.800 0.423574    
Educ3:L2               0.60680    0.13272   4.572 4.83e-06 ***
Educ4:L2               0.84863    0.19679   4.312 1.62e-05 ***
Educ2:L3               0.67195    0.22185   3.029 0.002455 ** 
Educ3:L3              -0.13777    0.12302  -1.120 0.262766    
Educ4:L3              -0.17785    0.18557  -0.958 0.337863    
Educ2:L4               0.83889    0.21766   3.854 0.000116 ***
Educ3:L4               0.25621    0.11970   2.140 0.032320 *  
Educ4:L4              -0.10081    0.19889  -0.507 0.612248    
Occ2:L2                1.22119    0.14272   8.557  < 2e-16 ***
Occ3:L2               -2.67793    0.26863  -9.969  < 2e-16 ***
Occ4:L2                0.25078    0.15354   1.633 0.102387    
Occ5:L2               -2.76550    0.74142  -3.730 0.000191 ***
Occ2:L3               -0.02264    0.15282  -0.148 0.882219    
Occ3:L3               -0.01696    0.13679  -0.124 0.901323    
Occ4:L3                0.12779    0.14886   0.858 0.390653    
Occ5:L3                0.04374    0.27494   0.159 0.873605    
Occ2:L4                0.83473    0.13514   6.177 6.54e-10 ***
Occ3:L4               -2.31754    0.20662 -11.217  < 2e-16 ***
Occ4:L4                0.05314    0.14388   0.369 0.711886    
Occ5:L4               -1.43905    0.38269  -3.760 0.000170 ***
Rscore:Cscore          0.80506    0.23183   3.473 0.000515 ***
Mult(., Educ, Occ).    0.69031         NA      NA       NA    
Mult(1, ., Occ).Educ1 -0.58680         NA      NA       NA    
Mult(1, ., Occ).Educ2 -0.05428         NA      NA       NA    
Mult(1, ., Occ).Educ3  0.48246         NA      NA       NA    
Mult(1, ., Occ).Educ4  0.74448         NA      NA       NA    
Mult(1, Educ, .).Occ1 -0.30028         NA      NA       NA    
Mult(1, Educ, .).Occ2  0.07222         NA      NA       NA    
Mult(1, Educ, .).Occ3 -0.20754         NA      NA       NA    
Mult(1, Educ, .).Occ4 -3.36324         NA      NA       NA    
Mult(1, Educ, .).Occ5 -5.98616         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: 124.75 on 41 degrees of freedom
AIC: 564.56

Number of iterations: 14

model.summary(model8)

# A tibble: 1 × 6
  `Model Description`    df    L2   BIC Delta        p
  <chr>               <int> <dbl> <dbl> <dbl>    <dbl>
1 model8                 41  125. -216.  5.58 2.18e-10

# Model 9 - U+RC (heterogeneous)
model9 <- freq_tab_3.1 |> 
  gnm(
    Freq ~ Educ + Occ + L + Educ:L + Occ:L + 
      Rscore:Cscore:L + 
      Mult(L, L:Educ, L:Occ),
    family = poisson,
    data = _,
    trace = F,
    tolerance = 1e-6,
    iterStart = 20,
    iterMax = 100000,
    verbose = F
  )
model.summary(model9)

# A tibble: 1 × 6
  `Model Description`    df    L2   BIC Delta     p
  <chr>               <int> <dbl> <dbl> <dbl> <dbl>
1 model9                 20  27.5 -139.  1.50 0.122

4.2 表4.4

models <- list()
models[[1]] <- model.summary(model1)
models[[2]] <- model.summary(model2)
models[[3]] <- model.summary(model3)
models[[4]] <- model.summary(model4)
models[[5]] <- model.summary(model5)
models[[6]] <- model.summary(model6)
models[[7]] <- model.summary(model7)
models[[8]] <- model.summary(model8)
models[[9]] <- model.summary(model9)
models |> bind_rows() |> kable(digit = 3)

Model Description	df	L2	BIC	Delta	p
model1	48	1370.934	971.892	23.750	0.000
model2	42	143.827	-205.334	6.388	0.000
model3	24	69.062	-130.459	3.039	0.000
model4	38	117.384	-198.524	5.237	0.000
model5	8	5.834	-60.673	0.437	0.666
model6	36	113.180	-186.101	5.187	0.000
model7	0	0.000	0.000	0.000	0.000
model8	41	124.751	-216.097	5.581	0.000
model9	20	27.472	-138.795	1.496	0.122

4.3 表4.5

model_comparison <- function(x, y = 0) {
  models |> 
    bind_rows() |>
    dplyr::summarise(`Model Used` = 
                ifelse(y == 0,
                       paste0(x),
                       paste0(x,"-",y)),
              df = ifelse(y == 0, 
                          df[x], 
                          df[x] - df[y]),
              L2 = ifelse(y == 0, 
                          L2[x], 
                          L2[x] - L2[y]))
}

表4.5Aは次のように再現できる．

# Table 4.5 Panel A
bind_rows(model_comparison(1,2),
          model_comparison(2,4),
          model_comparison(4,6),
          model_comparison(6),
          model_comparison(1)) |> kable(digit = 3)

Model Used	df	L2
1-2	6	1227.106
2-4	4	26.444
4-6	2	4.204
6	36	113.180
1	48	1370.934

# Table 4.5 Panel B
bind_rows(model_comparison(1,3),
          model_comparison(3,5),
          model_comparison(5,7),
          model_comparison(1)) |> kable(digit = 3)

Model Used	df	L2
1-3	24	1301.872
3-5	16	63.228
5-7	8	5.834
1	48	1370.934

4.4 表4.6

# Table 4.6
Freq <- c(201, 29, 8,13, 5, 152,29, 2, 8, 0,
          18, 6, 3, 6, 0,17,12, 0, 3, 0,
          109,74, 164,89,16, 101, 336, 9, 134, 2,
          7, 6,45,30, 6, 7,41, 7,63, 0,
          247,58,20,23, 2, 288,51, 1,17, 3,
          48,11,16,13, 1,47,38, 2,18, 0,
          157,68, 178, 116,27, 165, 321,27, 168, 1,
          7, 7,50,42, 5,12,25, 5,29, 6)

Educ <- gl(4, 10, 4 * 5 * 2 * 2)
Occ <- gl(5, 1, 4 * 5 * 2 * 2)
Sex <- gl(2, 5, 4 * 5 * 2 * 2)
Year <- gl(2, 40, 4 * 5 * 2 * 2)
L <- Year:Sex
levels(L) <- factor(1:4)
Rscore <- as.numeric(Educ)
Cscore <- as.numeric(Occ)

# Model 1
model1.un <- freq_tab_3.1 |> gnm(
  Freq ~ Educ + Occ + L + Educ:L + Occ:L +
    Mult(L, Educ, L:Occ, inst = 1) +
    Mult(L, Educ, L:Occ, inst = 2),
  family = poisson,
  data = _,
  trace = F,
  tolerance = 1e-8,
  iterStart = 20,
  iterMax = 100000,
  verbose = F)
model.summary(model1.un)

# A tibble: 1 × 6
  `Model Description`    df    L2   BIC Delta     p
  <chr>               <int> <dbl> <dbl> <dbl> <dbl>
1 model1.un              14  7.97 -108. 0.687 0.891

# Model 2
model2.un <- freq_tab_3.1 |> gnm(
  Freq ~ Educ + Occ + L + Educ:L + Occ:L + 
    Mult(L, L:Educ, Occ, inst = 1) + 
    Mult(L, L:Educ, Occ, inst = 2),
    family = poisson,
    data = _,
    trace = F,
    tolerance = 1e-8,
    iterStart = 20,
    iterMax = 100000,
    verbose = F
  )
model.summary(model2.un)

# A tibble: 1 × 6
  `Model Description`    df    L2   BIC Delta      p
  <chr>               <int> <dbl> <dbl> <dbl>  <dbl>
1 model2.un              20  29.2 -137.  1.82 0.0830

# Model 3
# it might not be fitted with gnm!

model2.un <- freq_tab_3.1 |> gnm(
  Freq ~ Educ + Occ + L + Educ:L + Occ:L + 
    Mult(L, Educ, Occ, inst = 1) + 
    Mult(L, Educ, Occ, inst = 2),
    family = poisson,
    data = _,
    trace = F,
    tolerance = 1e-8,
    iterStart = 20,
    iterMax = 100000,
    verbose = F
  )
model.summary(model2.un)

# A tibble: 1 × 6
  `Model Description`    df    L2   BIC Delta     p
  <chr>               <int> <dbl> <dbl> <dbl> <dbl>
1 model2.un              30  38.5 -211.  2.12 0.138

# Model 4
# it might not be fitted with gnm!

# Model 5
# it might not be fitted with gnm!

# Model 6
# it might not be fitted with gnm!

# Model 7
# it might not be fitted with gnm!

# Model 8
# it might not be fitted with gnm!

# Model 9
model9.un <- gnm(
  Freq ~ Educ + Occ + L + Educ:L + Occ:L + 
    Mult(L, Educ, Occ, inst = 1) + 
    Mult(L, Educ, Occ, inst = 2),
  family = poisson,
  trace = F,
  tolerance = 1e-12,
  iterStart = 20,
  iterMax = 100000
)

Initialising
Running start-up iterations....................
Running main iterations.........................................................
......................................................
Done

model.summary(model9.un)

# A tibble: 1 × 6
  `Model Description`    df    L2   BIC Delta     p
  <chr>               <int> <dbl> <dbl> <dbl> <dbl>
1 model9.un              30  38.5 -211.  2.12 0.138

mu1 <- getContrasts(
  model9.un,
  pickCoef(model9.un, "[.]Educ")[1:4],
  ref = "mean",
  scaleRef = "mean",
  scaleWeights = "unit"
)

nu1 <- getContrasts(
  model9.un,
  pickCoef(model9.un, "[.]Occ")[1:5],
  ref = "mean",
  scaleRef = "mean",
  scaleWeights = "unit"
)

mu2 <- getContrasts(
  model9.un,
  pickCoef(model9.un, "[.]Educ")[5:8],
  ref = "mean",
  scaleRef = "mean",
  scaleWeights = "unit"
)

nu2 <- getContrasts(
  model9.un,
  pickCoef(model9.un, "[.]Occ")[6:10],
  ref = "mean",
  scaleRef = "mean",
  scaleWeights = "unit"
)

con <- c(mu1$qvframe[, 1][c(1, 4)], 
         nu1$qvframe[, 1][c(1, 5)],
         mu2$qvframe[, 1][c(1, 4)], 
         nu2$qvframe[, 1][c(1, 5)])


set.seed(123456)
model9 <- gnm(
  Freq ~ Educ + Occ + L + Educ:L + Occ:L +
    Mult(L, Educ, Occ, inst = 1) + 
    Mult(L, Educ, Occ, inst = 2),
  constrain = c(37, 40, 41, 45, 50, 53, 54, 58),
  constrainTo = con,
  family = poisson,
  tolerance = 1e-12,
  iterStart = 20,
  iterMax = 100000
)

Initialising
Running start-up iterations....................
Running main iterations.........................................................
................................

Deviance is not finite

Warning: Algorithm failed - no model could be estimated

summary(model9)

Length  Class   Mode 
     0   NULL   NULL 

mu1

                                  Estimate Std. Error
Mult(L, ., Occ, inst = 1).Educ1 -0.6402222 0.02882211
Mult(L, ., Occ, inst = 1).Educ2 -0.2386164 0.03949527
Mult(L, ., Occ, inst = 1).Educ3  0.1683115 0.03139203
Mult(L, ., Occ, inst = 1).Educ4  0.7105272 0.02071164

nu1

                                  Estimate Std. Error
Mult(L, Educ, ., inst = 1).Occ1 -0.7654028 0.02524700
Mult(L, Educ, ., inst = 1).Occ2 -0.2501417 0.04589839
Mult(L, Educ, ., inst = 1).Occ3  0.3979180 0.05170619
Mult(L, Educ, ., inst = 1).Occ4  0.2733327 0.04831936
Mult(L, Educ, ., inst = 1).Occ5  0.3442938 0.08637688

mu2

                                  Estimate Std. Error
Mult(L, ., Occ, inst = 2).Educ1 -0.7309442 0.07862794
Mult(L, ., Occ, inst = 2).Educ2  0.2168538 0.14103742
Mult(L, ., Occ, inst = 2).Educ3  0.6355624 0.08287906
Mult(L, ., Occ, inst = 2).Educ4 -0.1214720 0.17120754

nu2

                                   Estimate Std. Error
Mult(L, Educ, ., inst = 2).Occ1  0.07059438 0.10141736
Mult(L, Educ, ., inst = 2).Occ2 -0.48030206 0.07332532
Mult(L, Educ, ., inst = 2).Occ3 -0.19772154 0.10545108
Mult(L, Educ, ., inst = 2).Occ4 -0.21626011 0.06608227
Mult(L, Educ, ., inst = 2).Occ5  0.82368933 0.03866984

model.summary(model9.un)

# A tibble: 1 × 6
  `Model Description`    df    L2   BIC Delta     p
  <chr>               <int> <dbl> <dbl> <dbl> <dbl>
1 model9.un              30  38.5 -211.  2.12 0.138

# Model 10
Lc <- L
levels(Lc) <- c(1, 1, 2, 2)
model10.un <- gnm(
  Freq ~ Educ + Occ + L + Educ:L + Occ:L + 
    Mult(Lc, Educ, Occ) + 
    Mult(L, Educ, Occ),
  family = poisson,
  trace = F,
  tolerance = 1e-8,
  iterStart = 20,
  iterMax = 100000,
  verbose = F
)
model.summary(model10.un)

# A tibble: 1 × 6
  `Model Description`    df    L2   BIC Delta     p
  <chr>               <int> <dbl> <dbl> <dbl> <dbl>
1 model10.un             32  40.4 -226.  2.17 0.146

# Model 11
set.seed(1234)
model11.un <- gnm(
  Freq ~ Educ + Occ + L + Educ:L + Occ:L + 
    Mult(1, Educ, Occ) + 
    Mult(L, Educ, Occ),
  family = poisson,
  trace = F,
  tolerance = 1e-8,
  iterStart = 20,
  iterMax = 100000,
  verbose = F
)
model.summary(model11.un)

# A tibble: 1 × 6
  `Model Description`    df    L2   BIC Delta     p
  <chr>               <int> <dbl> <dbl> <dbl> <dbl>
1 model11.un             33  43.4 -231.  2.31 0.106

# Model 12
L.c <- L
levels(L.c) <- factor(c(1, 2, 3, 2))

set.seed(1234)
model12.un <-
  gnm(Freq ~ Educ + Occ + L + Educ:L + Occ:L + 
        Mult(1, Educ, Occ) + 
        Mult(L.c, Educ, Occ),
    family = poisson,
    trace = F,
    tolerance = 1e-12,
    iterStart = 20,
    iterMax = 100000
  )

Initialising
Running start-up iterations....................
Running main iterations.........................................................
................................................................
Done

model.summary(model12.un)

# A tibble: 1 × 6
  `Model Description`    df    L2   BIC Delta     p
  <chr>               <int> <dbl> <dbl> <dbl> <dbl>
1 model12.un             34  44.6 -238.  2.70 0.105

mu1 <- getContrasts(
  model12.un,
  pickCoef(model12.un, "[.]Educ")[1:4],
  ref = "mean",
  scaleRef = "mean",
  scaleWeights = "unit"
)

nu1 <- getContrasts(
  model12.un,
  pickCoef(model12.un, "[.]Occ")[1:5],
  ref = "mean",
  scaleRef = "mean",
  scaleWeights = "unit"
)

mu2 <- getContrasts(
  model12.un,
  pickCoef(model12.un, "[.]Educ")[5:8],
  ref = "mean",
  scaleRef = "mean",
  scaleWeights = "unit"
)

nu2 <- getContrasts(
  model12.un,
  pickCoef(model12.un, "[.]Occ")[6:10],
  ref = "mean",
  scaleRef = "mean",
  scaleWeights = "unit"
)

con <- c(mu1$qvframe[, 1][c(1, 4)], nu1$qvframe[, 1][c(1, 5)],
         mu2$qvframe[, 1][c(1, 4)], nu2$qvframe[, 1][c(1, 5)])

model12 <-
  gnm(Freq ~ Educ + Occ + L + Educ:L + Occ:L + 
        Mult(1, Educ, Occ) + 
        Mult(L.c, Educ, Occ),
    constrain = c(34, 37, 38, 42, 46, 49, 50, 54),
    constrainTo = con,
    family = poisson,
    trace = F,
    tolerance = 1e-12,
    iterStart = 20,
    iterMax = 100000
  )

Initialising
Running start-up iterations....................
Running main iterations.........................................................
...................

Deviance is not finite

Warning: Algorithm failed - no model could be estimated

summary(model12)

Length  Class   Mode 
     0   NULL   NULL 

mu1

                        Estimate Std. Error
Mult(1, ., Occ).Educ1 -0.6496856 0.02689904
Mult(1, ., Occ).Educ2 -0.2266856 0.03916088
Mult(1, ., Occ).Educ3  0.1712502 0.02789996
Mult(1, ., Occ).Educ4  0.7051210 0.01869796

nu1

                        Estimate Std. Error
Mult(1, Educ, .).Occ1  0.7477472 0.03046619
Mult(1, Educ, .).Occ2  0.2798375 0.04956910
Mult(1, Educ, .).Occ3 -0.3926778 0.05281437
Mult(1, Educ, .).Occ4 -0.2590761 0.04731457
Mult(1, Educ, .).Occ5 -0.3758307 0.08654401

mu2

                          Estimate Std. Error
Mult(L.c, ., Occ).Educ1  0.7697514 0.05982573
Mult(L.c, ., Occ).Educ2 -0.1649990 0.14337525
Mult(L.c, ., Occ).Educ3 -0.6165380 0.08295724
Mult(L.c, ., Occ).Educ4  0.0117856 0.16680372

nu2

                           Estimate Std. Error
Mult(L.c, Educ, .).Occ1  0.04308284 0.07694408
Mult(L.c, Educ, .).Occ2 -0.47712634 0.07295173
Mult(L.c, Educ, .).Occ3 -0.22737747 0.10213888
Mult(L.c, Educ, .).Occ4 -0.16931700 0.06314427
Mult(L.c, Educ, .).Occ5  0.83073797 0.03112819

model.summary(model12.un)

# A tibble: 1 × 6
  `Model Description`    df    L2   BIC Delta     p
  <chr>               <int> <dbl> <dbl> <dbl> <dbl>
1 model12.un             34  44.6 -238.  2.70 0.105

# Model 13
model13.un <-
  gnm(Freq ~ Educ + Occ + L + Educ:L + Occ:L + 
        Mult(1, Educ, Occ, inst = 1) + 
        Mult(1, Educ, Occ, inst = 2),
    family = poisson,
    trace = F,
    tolerance = 1e-8,
    iterStart = 20,
    iterMax = 100000,
    verbose = F
  )
model.summary(model13.un)

# A tibble: 1 × 6
  `Model Description`    df    L2   BIC Delta        p
  <chr>               <int> <dbl> <dbl> <dbl>    <dbl>
1 model13.un             38  117. -199.  5.24 4.93e-10

library(logmult)
d <- tibble(Freq, Educ, Occ, Sex, Year, L, Rscore, Cscore)
tab <- xtabs(Freq ~ Educ + Occ + L, data = d)

# Model 2
fit <- rcL(tab, nd = 1,
    layer.effect = "none")

Initialising
Running start-up iterations..
Running main iterations................
Done

fit

Call:
rcL(tab = tab, nd = 1, layer.effect = "none")

Intrinsic association coefficients:
     Dim1 
[1,] 0.657

Normalized row scores (Educ) for all layers:
     1      2      3      4 
 1.428  0.416 -0.477 -1.701 

Normalized column scores (Occ) for all layers:
     1      2      3      4      5 
 1.171 -0.236 -1.272 -1.076 -1.053 

Normalization weights: marginal
Deviance:              143.8274
Pearson chi-squared:   152.4751
Dissimilarity index:   6.388286%
Residual df:           42
BIC:                   -205.3338
AIC:                   59.82743

# Model 3
fit <- rcL(tab, nd = 1,
    layer.effect = "heterogeneous")

Initialising
Running start-up iterations..
Running main iterations.........................................................
...................
Done

fit

Call:
rcL(tab = tab, nd = 1, layer.effect = "heterogeneous")

Intrinsic association coefficients:
  Dim1 
1 0.644
2 0.870
3 0.562
4 0.683

Normalized row scores (Educ):
, , 1

  Dim1  
1  1.457
2  0.520
3 -0.542
4 -1.488

, , 2

  Dim1  
1  1.395
2  0.588
3 -0.470
4 -1.759

, , 3

  Dim1  
1  1.403
2  0.447
3 -0.452
4 -1.793

, , 4

  Dim1  
1  1.534
2  0.159
3 -0.598
4 -1.161


Normalized column scores (Occ):
, , 1

  Dim1  
1  1.093
2 -0.024
3 -1.537
4 -1.069
5 -0.715

, , 2

  Dim1  
1  1.195
2 -0.341
3 -1.118
4 -1.141
5 -0.413

, , 3

  Dim1  
1  0.958
2  0.312
3 -1.415
4 -1.294
5 -1.417

, , 4

  Dim1  
1  1.185
2 -0.383
3 -1.466
4 -0.827
5 -0.337


Normalization weights: marginal
Deviance:              69.06197
Pearson chi-squared:   88.05995
Dissimilarity index:   3.038593%
Residual df:           24
BIC:                   -130.4587
AIC:                   21.06197

# Model 4
fit <- rcL(tab, nd = 2,
    layer.effect = "none")

Initialising
Running start-up iterations..
Running main iterations.................
Done

fit

Call:
rcL(tab = tab, nd = 2, layer.effect = "none")

Intrinsic association coefficients:
     Dim1   Dim2  
[1,] 0.6277 0.0863

Normalized row scores (Educ) for all layers:
  Dim1   Dim2  
1  1.458  0.480
2  0.425 -0.103
3 -0.520 -0.695
4 -1.549  2.607

Normalized column scores (Occ) for all layers:
  Dim1    Dim2   
1  1.1703  0.3450
2 -0.2552 -1.4735
3 -1.3890  0.0154
4 -0.9827  1.1830
5 -0.8680  2.3096

Normalization weights: marginal
Deviance:              117.3838
Pearson chi-squared:   122.1628
Dissimilarity index:   5.236715%
Residual df:           38
BIC:                   -198.524
AIC:                   41.38378

# Model 5
fit <- rcL(tab, nd = 2,
    layer.effect = "heterogeneous")

Initialising
Running start-up iterations..
Running main iterations.........................................................
................................................................................
................................................................................
................................................................................
................................................................................
................................................................................
................................................................................
................................................................................
................................................................................
................................................................................
...........................
Done

fit

Call:
rcL(tab = tab, nd = 2, layer.effect = "heterogeneous")

Intrinsic association coefficients:
  Dim1   Dim2  
1 0.6502 0.0943
2 4.2413 0.4108
3 0.5696 0.0757
4 0.6765 0.1914

Normalized row scores (Educ):
, , 1

  Dim1     Dim2    
1  1.45352  0.45162
2  0.51234 -0.02060
3 -0.53379 -0.69282
4 -1.51746  2.62495

, , 2

  Dim1     Dim2    
1  1.52661  0.24250
2 -0.00802  1.36857
3 -0.78150  0.23483
4  0.00888 -2.91553

, , 3

  Dim1     Dim2    
1  1.38340  0.71281
2  0.44660 -0.57412
3 -0.42823 -0.71739
4 -1.86901  2.36677

, , 4

  Dim1     Dim2    
1  1.54808  0.21241
2  0.13984 -0.58450
3 -0.63135 -0.54666
4 -0.99808  2.85541


Normalized column scores (Occ):
, , 1

  Dim1     Dim2    
1  1.06815  0.26848
2 -0.00474 -1.18599
3 -1.70563 -0.75487
4 -0.94238  1.42995
5 -0.56999  2.63870

, , 2

  Dim1     Dim2    
1  0.30521  1.09050
2  0.01328 -0.11121
3  0.10725 -1.32567
4 -0.02157 -1.23729
5 -7.28179  0.87125

, , 3

  Dim1     Dim2    
1  0.98305 -0.55082
2  0.24927  0.94783
3 -1.48810 -1.21427
4 -1.16872  0.97409
5 -1.78377 -3.81045

, , 4

  Dim1     Dim2    
1  1.15659 -0.16207
2 -0.42550 -0.63682
3 -1.48072  0.27185
4 -0.81835  0.40919
5  0.93490  6.79479


Normalization weights: marginal
Deviance:              5.833864
Pearson chi-squared:   4.354907
Dissimilarity index:   0.436722%
Residual df:           8
BIC:                   -60.67303
AIC:                   -10.16614

# Model 6
fit <- rcL(tab, nd = 3,
    layer.effect = "none")

Initialising
Running start-up iterations..
Running main iterations.....
Done

fit

Call:
rcL(tab = tab, nd = 3, layer.effect = "none")

Intrinsic association coefficients:
     Dim1   Dim2   Dim3  
[1,] 0.6287 0.0867 0.0499

Normalized row scores (Educ) for all layers:
  Dim1   Dim2   Dim3  
1  1.456  0.533  0.346
2  0.437 -0.631 -3.762
3 -0.521 -0.657  0.315
4 -1.549  2.584 -0.348

Normalized column scores (Occ) for all layers:
  Dim1    Dim2    Dim3   
1  1.1701  0.3410 -0.1004
2 -0.2525 -1.4473  0.3582
3 -1.3804  0.0413  0.1242
4 -0.9949  1.0643 -1.0180
5 -0.8408  3.0518  6.4839

Normalization weights: marginal
Deviance:              113.1797
Pearson chi-squared:   117.7575
Dissimilarity index:   5.186927%
Residual df:           36
BIC:                   -186.1013
AIC:                   41.17975

# Model 7
fit <- rcL(tab, nd = 3,
    layer.effect = "heterogeneous")

Initialising
Running start-up iterations..
Running main iterations.

Deviance is not finite

Warning: Algorithm failed - no model could be estimated

fit

NULL