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.

> # Avoid random failures to converge
> set.seed(1)
> 
> # See Becker-Clogg (1989) test
> 
> library(logmult)
Loading required package: gnm

Attaching package: 'logmult'

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

    se

> data(color)
> 
> rcm <- rc(color[,,1], 2, weighting="marginal", start=NA)
Running base model to find starting values...
Running real model...
Initialising
Running main iterations......................
Done
> rcu <- rc(color[,,1], 2, weighting="uniform", start=NA)
Running base model to find starting values...
Running real model...
Initialising
Running main iterations......................
Done
> rcn <- rc(color[,,1], 2, weighting="none", start=NA)
Running base model to find starting values...
Running real model...
Initialising
Running main iterations......................
Done
> 
> stopifnot(all.equal(fitted(rcm), fitted(rcu)))
> stopifnot(all.equal(fitted(rcm), fitted(rcn)))
> 
> phim <- iac(fitted(rcm), weighting="marginal")
> phiu <- iac(fitted(rcu), weighting="uniform")
> phin <- iac(fitted(rcn), weighting="none")
> 
> cphim <- iac(fitted(rcm), TRUE, weighting="marginal")
> cphiu <- iac(fitted(rcm), TRUE, weighting="uniform")
> cphin <- iac(fitted(rcm), TRUE, weighting="none")
> 
> stopifnot(all.equal(phim, sqrt(sum(rcm$assoc$phi^2))))
> stopifnot(all.equal(phiu, sqrt(sum(rcu$assoc$phi^2))))
> stopifnot(all.equal(phin, sqrt(sum(rcn$assoc$phi^2))))
> 
> stopifnot(all.equal(phim, sqrt(sum(cphim))))
> stopifnot(all.equal(phiu, sqrt(sum(cphiu))))
> stopifnot(all.equal(phin, sqrt(sum(cphin))))
> 
> 
> # Test on perfectly symmetric association
> data(ocg1973)
> 
> rcm <- rc(ocg1973, 2, symmetric=TRUE, weighting="marginal", start=NA)
Running base model to find starting values...
Running real model...
Initialising
Running main iterations...................................
Done
> rcu <- rc(ocg1973, 2, symmetric=TRUE, weighting="uniform", start=NA)
Running base model to find starting values...
Running real model...
Initialising
Running main iterations...................................
Done
> rcn <- rc(ocg1973, 2, symmetric=TRUE, weighting="none", start=NA)
Running base model to find starting values...
Running real model...
Initialising
Running main iterations...................................
Done
> 
> stopifnot(all.equal(fitted(rcm), fitted(rcu)))
> stopifnot(all.equal(fitted(rcm), fitted(rcn)))
> 
> w <- (rcm$assoc$row.weights + rcm$assoc$col.weights)[,1]/2
> phim <- iac(fitted(rcm), weighting="marginal", row.weights=w, col.weights=w)
> phiu <- iac(fitted(rcu), weighting="uniform")
> phin <- iac(fitted(rcn), weighting="none")
> 
> sphim <- iac(fitted(rcm), component="symmetric", weighting="marginal",
+               row.weights=w, col.weights=w)
> sphiu <- iac(fitted(rcu), component="symmetric", weighting="uniform")
> sphin <- iac(fitted(rcn), component="symmetric", weighting="none")
> 
> aphim <- iac(fitted(rcm), component="antisymmetric", weighting="marginal",
+              row.weights=w, col.weights=w)
> aphiu <- iac(fitted(rcu), component="antisymmetric", weighting="uniform")
> aphin <- iac(fitted(rcn), component="antisymmetric", weighting="none")
> 
> cphim <- iac(fitted(rcm), TRUE, weighting="marginal",
+              row.weights=w, col.weights=w)
> cphiu <- iac(fitted(rcu), TRUE, weighting="uniform")
> cphin <- iac(fitted(rcn), TRUE, weighting="none")
> 
> stopifnot(all.equal(phim, sphim))
> stopifnot(all.equal(phiu, sphiu))
> stopifnot(all.equal(phin, sphin))
> stopifnot(all.equal(phim, sqrt(sum(rcm$assoc$phi^2))))
> stopifnot(all.equal(phiu, sqrt(sum(rcu$assoc$phi^2))))
> stopifnot(all.equal(phin, sqrt(sum(rcn$assoc$phi^2))))
> 
> stopifnot(all(c(aphim, aphiu, aphin) < 1e-10))
> 
> stopifnot(all.equal(phim, sqrt(sum(cphim))))
> stopifnot(all.equal(phiu, sqrt(sum(cphiu))))
> stopifnot(all.equal(phin, sqrt(sum(cphin))))
> 
> 
> # Test on perfectly anti-symmetric association
> hmm <- hmskew(ocg1973, nd.symm=0, weighting="marginal", 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 
> hmu <- hmskew(ocg1973, nd.symm=0, 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 
> hmn <- hmskew(ocg1973, nd.symm=0, weighting="none", 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 
> 
> stopifnot(all.equal(fitted(hmm), fitted(hmu)))
> stopifnot(all.equal(fitted(hmm), fitted(hmn)))
> 
> w <- (hmm$assoc$row.weights + hmm$assoc$col.weights)[,1]/2
> phim <- iac(fitted(hmm), weighting="marginal", row.weights=w, col.weights=w)
> phiu <- iac(fitted(hmu), weighting="uniform")
> phin <- iac(fitted(hmn), weighting="none")
> 
> sphim <- iac(fitted(hmm), component="symmetric", weighting="marginal",
+               row.weights=w, col.weights=w)
> sphiu <- iac(fitted(hmu), component="symmetric", weighting="uniform")
> sphin <- iac(fitted(hmn), component="symmetric", weighting="none")
> 
> aphim <- iac(fitted(hmm), component="antisymmetric", weighting="marginal",
+               row.weights=w, col.weights=w)
> aphiu <- iac(fitted(hmu), component="antisymmetric", weighting="uniform")
> aphin <- iac(fitted(hmn), component="antisymmetric", weighting="none")
> 
> cphim <- iac(fitted(hmm), TRUE, weighting="marginal",
+               row.weights=w, col.weights=w)
> cphiu <- iac(fitted(hmu), TRUE, weighting="uniform")
> cphin <- iac(fitted(hmn), TRUE, weighting="none")
> 
> stopifnot(all.equal(phim, aphim))
> stopifnot(all.equal(phiu, aphiu))
> stopifnot(all.equal(phin, aphin))
> stopifnot(all.equal(phim, sqrt(sum(hmm$assoc$phi^2))))
> stopifnot(all.equal(phiu, sqrt(sum(hmu$assoc$phi^2))))
> stopifnot(all.equal(phin, sqrt(sum(hmn$assoc$phi^2))))
> 
> stopifnot(all(c(sphim, sphiu, sphin) < 1e-10))
> 
> stopifnot(all.equal(phim, sqrt(sum(cphim))))
> stopifnot(all.equal(phiu, sqrt(sum(cphiu))))
> stopifnot(all.equal(phin, sqrt(sum(cphin))))
> 
> 
> # Test on symmetric and anti-symmetric association
> # Without starting values, model too often converges to wrong solution
> start <- c(6.540,  0.106,  0.407, 0.666, 1.006, -0.581, -0.261, 0.060,
+           -4.411, -0.567, -0.310, 0.264, 0.652, -1.794, -1.610, -1.627,
+           -0.743, -0.012,  6.311, 0.295, 0.198, -0.015, -0.167, 0.010)
> hmm <- hmskew(ocg1973, nd.symm=1, weighting="marginal", start=start)
Running main iterations................
Done
Warning message:
Using make.unique() to make default parameter labels unique 
> hmu <- hmskew(ocg1973, nd.symm=1, weighting="uniform", start=start)
Running main iterations................
Done
Warning message:
Using make.unique() to make default parameter labels unique 
> hmn <- hmskew(ocg1973, nd.symm=1, weighting="none", start=start)
Running main iterations................
Done
Warning message:
Using make.unique() to make default parameter labels unique 
> 
> stopifnot(all.equal(fitted(hmm), fitted(hmu)))
> stopifnot(all.equal(fitted(hmm), fitted(hmn)))
> 
> w <- (hmm$assoc$row.weights + hmm$assoc$col.weights)[,1]/2
> phim <- iac(fitted(hmm), weighting="marginal", row.weights=w, col.weights=w)
> phiu <- iac(fitted(hmu), weighting="uniform")
> phin <- iac(fitted(hmn), weighting="none")
> 
> sphim <- iac(fitted(hmm), component="symmetric", weighting="marginal",
+               row.weights=w, col.weights=w)
> sphiu <- iac(fitted(hmu), component="symmetric", weighting="uniform")
> sphin <- iac(fitted(hmn), component="symmetric", weighting="none")
> 
> aphim <- iac(fitted(hmm), component="antisymmetric", weighting="marginal",
+              row.weights=w, col.weights=w)
> aphiu <- iac(fitted(hmu), component="antisymmetric", weighting="uniform")
> aphin <- iac(fitted(hmn), component="antisymmetric", weighting="none")
> 
> cphim <- iac(fitted(hmm), TRUE, weighting="marginal",
+              row.weights=w, col.weights=w)
> cphiu <- iac(fitted(hmu), TRUE, weighting="uniform")
> cphin <- iac(fitted(hmn), TRUE, weighting="none")
> 
> stopifnot(all.equal(phim, sqrt(sphim^2 + aphim^2)))
> stopifnot(all.equal(phiu, sqrt(sphiu^2 + aphiu^2)))
> stopifnot(all.equal(phin, sqrt(sphin^2 + aphin^2)))
> stopifnot(all.equal(sphim, sqrt(sum(hmm$assoc$phi^2))))
> stopifnot(all.equal(sphiu, sqrt(sum(hmu$assoc$phi^2))))
> stopifnot(all.equal(sphin, sqrt(sum(hmn$assoc$phi^2))))
> stopifnot(all.equal(aphim, sqrt(sum(hmm$assoc.hmskew$phi^2))))
> stopifnot(all.equal(aphiu, sqrt(sum(hmu$assoc.hmskew$phi^2))))
> stopifnot(all.equal(aphin, sqrt(sum(hmn$assoc.hmskew$phi^2))))
> 
> stopifnot(all.equal(phim, sqrt(sum(cphim))))
> stopifnot(all.equal(phiu, sqrt(sum(cphiu))))
> stopifnot(all.equal(phin, sqrt(sum(cphin))))
> 
> 
> # Test for phi computed from UNIDIFF two-way interaction coefficients
> data(yaish)
> tab <- aperm(yaish[,,-7], 3:1)
> 
> u2m <- unidiff(tab, weighting="marginal")
Initialising
Running start-up iterations..
Running main iterations.........................................................
........
Done
> stopifnot(all.equal(u2m$unidiff$phi,
+                     iac(fitted(u2m), weighting="marginal")))
> stopifnot(all.equal(u2m$unidiff$phi[1] * exp(u2m$unidiff$layer$qvframe$estimate),
+                     iac(fitted(u2m), weighting="marginal"), check.attributes=FALSE))
> 
> u2u <- unidiff(tab, weighting="uniform")
Initialising
Running start-up iterations..
Running main iterations.........................................................
........
Done
> stopifnot(all.equal(u2u$unidiff$phi,
+                     iac(fitted(u2u), weighting="uniform")))
> stopifnot(all.equal(u2u$unidiff$phi[1] * exp(u2u$unidiff$layer$qvframe$estimate),
+                     iac(fitted(u2u), weighting="uniform"), check.attributes=FALSE))
> 
> u2n <- unidiff(tab, weighting="none", norm=2)
Initialising
Running start-up iterations..
Running main iterations.........................................................
....................
Done
> stopifnot(all.equal(u2n$unidiff$phi,
+                     iac(fitted(u2n), weighting="none")))
> stopifnot(all.equal(u2n$unidiff$phi[1] * exp(u2n$unidiff$layer$qvframe$estimate),
+                     iac(fitted(u2n), weighting="none"), check.attributes=FALSE))
> 
> stopifnot(all.equal(iac(yaish + 0.5),
+                     apply(yaish + 0.5, 3, iac,
+                           row.weights=margin.table(yaish + 0.5, 1),
+                           col.weights=margin.table(yaish + 0.5, 2))))
> 
> 
> # Test with all layers having proportional frequencies,
> # where standard and shrunk estimators must give the same result
> data(yaish)
> tab <- yaish[,,c(1, 1, 1)] + 0.5
> tab[,,2] <- tab[,,2] * 3
> tab[,,3] <- tab[,,3] * 8
> 
> phim <- iac(tab, weighting="marginal")
> phims <- iac(tab, weighting="marginal", shrink=TRUE)
> stopifnot(all.equal(phim, phims))
> 
> phiu <- iac(tab, weighting="uniform")
> phius <- iac(tab, weighting="uniform", shrink=TRUE)
> stopifnot(all.equal(phiu, phius))
> 
> phin <- iac(tab, weighting="none")
> phins <- iac(tab, weighting="none", shrink=TRUE)
> stopifnot(all.equal(phin, phins))
> 
> 
> ###
> ## Comparison with mean/sum of all spanning odds ratios
> ###
> # Can be set to arbitrary values
> nr <- 4
> nc <- 5
> 
> or <- function(tab) {
+     or1 <- function(i, j, i2, j2) (tab[i, j] * tab[i2, j2]) / (tab[i, j2] * tab[i2, j])
+     or <- array(NA, c(nrow(tab), ncol(tab), nrow(tab), ncol(tab)))
+ 
+     for(i in 1:nrow(tab))
+       for(j in 1:ncol(tab))
+         for(i2 in 1:nrow(tab))
+           for(j2 in 1:ncol(tab))
+               or[i, j, i2, j2] <- or1(i, j, i2, j2)
+     or
+ }
> 
> wlor2 <- function(tab) {
+     rp <- prop.table(margin.table(tab, 1)) * nrow(tab)
+     cp <- prop.table(margin.table(tab, 2)) * ncol(tab)
+     wlor2 <- w <- array(NA, c(nrow(tab), ncol(tab), nrow(tab), ncol(tab)))
+ 
+     for(i in 1:nrow(tab))
+       for(j in 1:ncol(tab))
+         for(i2 in 1:nrow(tab))
+           for(j2 in 1:ncol(tab)){
+               wlor2[i, j, i2, j2] <- log((tab[i, j] * tab[i2, j2]) / (tab[i, j2] * tab[i2, j]))^2
+               w[i, j, i2, j2] <- rp[i] * cp[j] * rp[i2] * cp[j2]
+           }
+ 
+     wlor2 * w / sum(w)
+ }
> 
> ## Uniform weighting
> # General case
> res <- replicate(10, {
+     tab <- matrix(rpois(nr*nc, 1000), nr, nc) + .5
+     rp <- rep(1/nr, nr)
+     cp <- rep(1/nc, nc)
+     c(iac(tab, weighting="uniform"),
+       iac(tab, row.weights=rp, col.weights=cp),
+       sqrt(sum(log(or(tab))^2)/((nr*nc)^2))/2)
+ })
> 
> stopifnot(all.equal(res[1,], res[2,]))
> stopifnot(all.equal(res[1,], res[3,]))
> 
> # 2x2 table (relation with single odds ratio)
> res <- replicate(10, {
+     tab <- matrix(rpois(2*2, 1000), 2, 2) + .5
+     rp <- cp <- rep(1/2, 2)
+     c(iac(tab, weighting="uniform"),
+       iac(tab, row.weights=rp, col.weights=cp),
+       abs(log(tab[1,1] * tab[2,2] / (tab[1,2] * tab[2,1])))/4)
+ })
> 
> stopifnot(all.equal(res[1,], res[2,]))
> stopifnot(all.equal(res[1,], res[3,]))
> 
> 
> ## No weighting
> # General case
> res <- replicate(10, {
+     tab <- matrix(rpois(nr*nc, 1000), nr, nc) + .5
+     rp <- rep(1, nr)
+     cp <- rep(1, nc)
+     c(iac(tab, weighting="none"),
+       sqrt(sum(log(or(tab))^2)/(nr*nc))/2)
+ })
> 
> stopifnot(all.equal(res[1,], res[2,]))
> 
> # 2x2 table (relation with single log-odds ratio)
> res <- replicate(10, {
+     tab <- matrix(rpois(2*2, 1000), 2, 2) + .5
+     rp <- cp <- rep(1, 2)
+     c(iac(tab, weighting="none"),
+       abs(log(tab[1,1] * tab[2,2] / (tab[1,2] * tab[2,1])))/2)
+ })
> 
> stopifnot(all.equal(res[1,], res[2,]))
> 
> 
> ## Marginal weighting
> # General case
> res <- replicate(10, {
+     tab <- matrix(rpois(nr*nc, 100), nr, nc) +.5
+     rp <- prop.table(margin.table(tab, 1))
+     cp <- prop.table(margin.table(tab, 2))
+     c(iac(tab, weighting="marginal"),
+       iac(tab, row.weights=rp, col.weights=cp),
+       sqrt(sum(wlor2(tab)))/2)
+ })
> 
> stopifnot(all.equal(res[1,], res[2,]))
> stopifnot(all.equal(res[1,], res[3,]))
> 
> # 2x2 table (equality with single odds ratio)
> res <- replicate(10, {
+     tab <- matrix(rpois(2*2, 100), 2, 2) + .5
+     rp <- prop.table(margin.table(tab, 1))
+     cp <- prop.table(margin.table(tab, 2))
+     c(iac(tab, weighting="marginal"),
+       iac(tab, row.weights=rp, col.weights=cp),
+       sqrt(log(tab[1,1] * tab[2,2] / (tab[1,2] * tab[2,1]))^2 *
+            rp[1] * cp[1] * rp[2] * cp[2]))
+ })
> 
> stopifnot(all.equal(res[1,], res[2,]))
> stopifnot(all.equal(res[1,], res[3,]))
> 
> 
> ## Normalized IAC
> data(color)
> 
> stopifnot(abs(iac(color[,,1], normalize=TRUE) - 0.4203209) < 1e-7)
> stopifnot(abs(iac(color[,,2], normalize=TRUE) - 0.3646107) < 1e-7)
> stopifnot(abs(iac(color[,,1], normalize=TRUE, weighting="uniform") - 0.5105931) < 1e-7)
> stopifnot(abs(iac(color[,,2], normalize=TRUE, weighting="uniform") - 0.5275346) < 1e-7)
> stopifnot(abs(iac(color[,,1], normalize=TRUE, weighting="none") - 0.8511365) < 1e-7)
> stopifnot(abs(iac(color[,,2], normalize=TRUE, weighting="none") - 0.8586745) < 1e-7)
> 
> stopifnot(all.equal(iac(color, normalize=TRUE),
+                     c(0.4032049, 0.3698703), check.names=FALSE,
+                     tolerance=1e-7))
> stopifnot(all.equal(iac(color, normalize=TRUE, weighting="uniform"),
+                     c(0.5105931, 0.5275346), check.names=FALSE,
+                     tolerance=1e-7))
> stopifnot(all.equal(iac(color, normalize=TRUE, weighting="none"),
+                     c(0.8511365, 0.8586745), check.names=FALSE,
+                     tolerance=1e-7))
> 
> stopifnot(all.equal(iac(color, normalize=TRUE, shrink=TRUE),
+                     c(0.4114786, 0.3717251), check.names=FALSE,
+                     tolerance=1e-7))
> stopifnot(all.equal(iac(color, normalize=TRUE, shrink=TRUE, weighting="uniform"),
+                     c(0.5188945, 0.5596786), check.names=FALSE,
+                     tolerance=1e-7))
> stopifnot(all.equal(iac(color, normalize=TRUE, shrink=TRUE, weighting="none"),
+                     c(0.854869, 0.872177), check.names=FALSE,
+                     tolerance=1e-7))
> 
> stopifnot(iac(matrix(1, 2, 2), normalize=TRUE) == 0)
> stopifnot(iac(array(1, c(2, 2, 2)), normalize=TRUE) == c(0, 0))
> stopifnot(iac(array(1, c(2, 2, 2)), normalize=TRUE, shrink=TRUE) == c(0, 0))
> 
> stopifnot(tryCatch(iac(color, cell=TRUE, normalize=TRUE),
+                    error=function(e) TRUE))
> 
> proc.time()
   user  system elapsed 
   4.21    0.39    4.59