R Under development (unstable) (2025-08-24 r88696 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2025 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> ## van der Heijden & Mooijaart (1995), Table 2c, p. 23
> # See also ?hmskew
> 
> library(logmult)
Loading required package: gnm

Attaching package: 'logmult'

The following object is masked from 'package:gnm':

    se

> data(ocg1973)
> 
> # 5:1 is here to take "Farmers" as reference category (angle 0)
> model <- hmskew(ocg1973[5:1, 5:1], weighting="uniform", start=NA)
Running base model to find starting values...
Running real model...
Initialising
Running main iterations...........
Done
Warning message:
Using make.unique() to make default parameter labels unique 
> model
Call:
hmskew(tab = ocg1973[5:1, 5:1], weighting = "uniform", start = NA)

Intrinsic skew association coefficients:
  Dim1   Dim2 
0.0976 0.0976 

Normalized skew association scores:
    Dim1   Dim2  
F    1.967  0.000
LM  -0.347 -1.639
UM  -0.742 -0.423
LNM -0.630  0.970
UNM -0.249  1.092

Normalization weights: uniform
Deviance:              0.08265023
Pearson chi-squared:   0.08267867
Dissimilarity index:   0.05559037%
Residual df:           1
BIC:                   -9.728558
AIC:                   -1.91735
> ass <- model$assoc.hmskew
> 
> # First column of the table
> round(ass$row[,,1] * sqrt(ass$phi[1,1]), d=2)[5:1,]
     Dim1  Dim2
UNM -0.08  0.34
LNM -0.20  0.30
UM  -0.23 -0.13
LM  -0.11 -0.51
F    0.61  0.00
> 
> summary(model)
Call:
hmskew(tab = ocg1973[5:1, 5:1], weighting = "uniform", start = NA)

Deviance Residuals:
     D
O              F         LM         UM        LNM        UNM
  F     0.00e+00   4.66e-03  -1.17e-02   1.96e-02  -9.68e-03
  LM   -2.61e-02  -4.71e-08   6.86e-03  -5.04e-02   3.99e-02
  UM    5.47e-02  -8.09e-03   0.00e+00   9.26e-02  -7.32e-02
  LNM  -7.79e-02   6.74e-02  -1.13e-01  -1.97e-07   3.88e-02
  UNM   6.09e-02  -8.21e-02   1.43e-01  -6.59e-02  -3.47e-07

Association coefficients:
             Normalized  Adjusted
Dim1           0.097622        NA
Dim2           0.097622        NA
Dim1:O|DF      1.967392    0.6147
Dim1:O|DLM    -0.346896   -0.1084
Dim1:O|DUM    -0.741997   -0.2318
Dim1:O|DLNM   -0.629767   -0.1968
Dim1:O|DUNM   -0.248732   -0.0777
Dim2:O|DF      0.000000    0.0000
Dim2:O|DLM    -1.639087   -0.5121
Dim2:O|DUM    -0.423322   -0.1323
Dim2:O|DLNM    0.970278    0.3032
Dim2:O|DUNM    1.092132    0.3412

No standard errors available
(jackknife and bootstrap disabled, or weighting changed since fitting).

Normalization weights: uniform
Deviance:              0.08265023
Pearson chi-squared:   0.08267867
Dissimilarity index:   0.05559037%
Residual df:           1
BIC:                   -9.728558
AIC:                   -1.91735
> 
> # First score for Farmers is slightly different from the original article
> stopifnot(isTRUE(all.equal(round(ass$row[,,1] * sqrt(ass$phi[1,1]), d=2)[5:1,],
+                                  matrix(c(-0.08, -0.2, -0.23, -0.11, 0.61,
+                                            0.34,  0.3, -0.13, -0.51, 0), 5, 2),
+                           check.attributes=FALSE)))
> 
> # Right part of the table
> round(ass$phi[1] * (ass$row[,2,1] %o% ass$row[,1,1] - ass$row[,1,1] %o% ass$row[,2,1]), d=3)[5:1, 5:1]
       UNM    LNM     UM     LM      F
UNM  0.000 -0.044 -0.089 -0.077  0.210
LNM  0.044  0.000 -0.096 -0.134  0.186
UM   0.089  0.096  0.000 -0.104 -0.081
LM   0.077  0.134  0.104  0.000 -0.315
F   -0.210 -0.186  0.081  0.315  0.000
> 
> # Plot
> plot(model, coords="cartesian")
> 
> # Test anova
> indep <- gnm(Freq ~ O + D, data=ocg1973, family=poisson)
> symm <- gnm(Freq ~ O + D + Symm(O, D), data=ocg1973, family=poisson)
> anova(indep, symm, model, test="LR")
Analysis of Deviance Table

Model 1: Freq ~ O + D
Model 2: Freq ~ O + D + Symm(O, D)
Model 3: Freq ~ O + D + Symm(O, D) + HMSkew(O, D)
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1        16     3278.7                          
2         6       46.2 10   3232.5 < 2.2e-16 ***
3         1        0.1  5     46.2 8.462e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> proc.time()
   user  system elapsed 
   1.07    0.17    1.23