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 <- maor(fitted(rcm), TRUE, weighting="marginal", norm=2)
> phiu <- maor(fitted(rcu), TRUE, weighting="uniform", norm=2)
> phin <- maor(fitted(rcn), TRUE, weighting="none", norm=2)
> 
> cphim <- maor(fitted(rcm), TRUE, TRUE, weighting="marginal", norm=2)
> cphiu <- maor(fitted(rcm), TRUE, TRUE, weighting="uniform", norm=2)
> cphin <- maor(fitted(rcm), TRUE, TRUE, weighting="none", norm=2)
> 
> maorm <- maor(fitted(rcm), weighting="marginal", norm=2)
> maoru <- maor(fitted(rcu), weighting="uniform", norm=2)
> maorn <- maor(fitted(rcn), weighting="none", norm=2)
> 
> cmaorm <- maor(fitted(rcm), cell=TRUE, weighting="marginal", norm=2)
> cmaoru <- maor(fitted(rcu), cell=TRUE, weighting="uniform", norm=2)
> cmaorn <- maor(fitted(rcn), cell=TRUE, weighting="none", norm=2)
> 
> stopifnot(all.equal(phim, sqrt(sum((rcm$assoc$phi)^2))))
> stopifnot(all.equal(phiu, sqrt(sum(abs(rcu$assoc$phi)^2))))
> stopifnot(all.equal(phin, sqrt(sum(abs(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))))
> 
> stopifnot(all.equal(maorm, exp(sqrt(sum(cmaorm)))))
> stopifnot(all.equal(maoru, exp(sqrt(sum(cmaoru)))))
> stopifnot(all.equal(maorn, exp(sqrt(sum(cmaorn)))))
> 
> 
> # 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 <- maor(fitted(rcm), TRUE, weighting="marginal", norm=2, row.weights=w, col.weights=w)
> phiu <- maor(fitted(rcu), TRUE, weighting="uniform", norm=2)
> phin <- maor(fitted(rcn), TRUE, weighting="none", norm=2)
> 
> sphim <- maor(fitted(rcm), TRUE, component="symmetric", weighting="marginal", norm=2,
+               row.weights=w, col.weights=w)
> sphiu <- maor(fitted(rcu), TRUE, component="symmetric", weighting="uniform", norm=2)
> sphin <- maor(fitted(rcn), TRUE, component="symmetric", weighting="none", norm=2)
> 
> aphim <- maor(fitted(rcm), TRUE, component="antisymmetric", weighting="marginal", norm=2,
+               row.weights=w, col.weights=w)
> aphiu <- maor(fitted(rcu), TRUE, component="antisymmetric", weighting="uniform", norm=2)
> aphin <- maor(fitted(rcn), TRUE, component="antisymmetric", weighting="none", norm=2)
> 
> cphim <- maor(fitted(rcm), TRUE, TRUE, weighting="marginal", norm=2,
+               row.weights=w, col.weights=w)
> cphiu <- maor(fitted(rcu), TRUE, TRUE, weighting="uniform", norm=2)
> cphin <- maor(fitted(rcn), TRUE, TRUE, weighting="none", norm=2)
> 
> maorm <- maor(fitted(rcm), weighting="marginal", norm=2,
+               row.weights=w, col.weights=w)
> maoru <- maor(fitted(rcu), weighting="uniform", norm=2)
> maorn <- maor(fitted(rcn), weighting="none", norm=2)
> 
> smaorm <- maor(fitted(rcm), component="symmetric", weighting="marginal", norm=2,
+                row.weights=w, col.weights=w)
> smaoru <- maor(fitted(rcu), component="symmetric", weighting="uniform", norm=2)
> smaorn <- maor(fitted(rcn), component="symmetric", weighting="none", norm=2)
> 
> amaorm <- maor(fitted(rcm), component="antisymmetric", weighting="marginal", norm=2,
+                row.weights=w, col.weights=w)
> amaoru <- maor(fitted(rcu), component="antisymmetric", weighting="uniform", norm=2)
> amaorn <- maor(fitted(rcn), component="antisymmetric", weighting="none", norm=2)
> 
> cmaorm <- maor(fitted(rcm), cell=TRUE, weighting="marginal", norm=2,
+                row.weights=w, col.weights=w)
> cmaoru <- maor(fitted(rcu), cell=TRUE, weighting="uniform", norm=2)
> cmaorn <- maor(fitted(rcn), cell=TRUE, weighting="none", norm=2)
> 
> 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(abs(rcu$assoc$phi)^2))))
> stopifnot(all.equal(phin, sqrt(sum(abs(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))))
> 
> stopifnot(all(c(amaorm, amaoru, amaorn) - 1 < 1e-10))
> 
> stopifnot(all.equal(maorm, exp(sqrt(sum(cmaorm)))))
> stopifnot(all.equal(maoru, exp(sqrt(sum(cmaoru)))))
> stopifnot(all.equal(maorn, exp(sqrt(sum(cmaorn)))))
> 
> stopifnot(all.equal(maorm, smaorm))
> stopifnot(all.equal(maoru, smaoru))
> stopifnot(all.equal(maorn, smaorn))
> 
> # 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 <- maor(fitted(hmm), TRUE, weighting="marginal", norm=2, row.weights=w, col.weights=w)
> phiu <- maor(fitted(hmu), TRUE, weighting="uniform", norm=2)
> phin <- maor(fitted(hmn), TRUE, weighting="none", norm=2)
> 
> sphim <- maor(fitted(hmm), TRUE, component="symmetric", weighting="marginal", norm=2,
+               row.weights=w, col.weights=w)
> sphiu <- maor(fitted(hmu), TRUE, component="symmetric", weighting="uniform", norm=2)
> sphin <- maor(fitted(hmn), TRUE, component="symmetric", weighting="none", norm=2)
> 
> aphim <- maor(fitted(hmm), TRUE, component="antisymmetric", weighting="marginal", norm=2,
+               row.weights=w, col.weights=w)
> aphiu <- maor(fitted(hmu), TRUE, component="antisymmetric", weighting="uniform", norm=2)
> aphin <- maor(fitted(hmn), TRUE, component="antisymmetric", weighting="none", norm=2)
> 
> cphim <- maor(fitted(hmm), TRUE, TRUE, weighting="marginal", norm=2,
+               row.weights=w, col.weights=w)
> cphiu <- maor(fitted(hmu), TRUE, TRUE, weighting="uniform", norm=2)
> cphin <- maor(fitted(hmn), TRUE, TRUE, weighting="none", norm=2)
> 
> maorm <- maor(fitted(hmm), weighting="marginal", norm=2,
+               row.weights=w, col.weights=w)
> maoru <- maor(fitted(hmu), weighting="uniform", norm=2)
> maorn <- maor(fitted(hmn), weighting="none", norm=2)
> 
> smaorm <- maor(fitted(hmm), component="symmetric", weighting="marginal", norm=2,
+                row.weights=w, col.weights=w)
> smaoru <- maor(fitted(hmu), component="symmetric", weighting="uniform", norm=2)
> smaorn <- maor(fitted(hmn), component="symmetric", weighting="none", norm=2)
> 
> amaorm <- maor(fitted(hmm), component="antisymmetric", weighting="marginal", norm=2,
+                row.weights=w, col.weights=w)
> amaoru <- maor(fitted(hmu), component="antisymmetric", weighting="uniform", norm=2)
> amaorn <- maor(fitted(hmn), component="antisymmetric", weighting="none", norm=2)
> 
> cmaorm <- maor(fitted(hmm), cell=TRUE, weighting="marginal", norm=2,
+                row.weights=w, col.weights=w)
> cmaoru <- maor(fitted(hmu), cell=TRUE, weighting="uniform", norm=2)
> cmaorn <- maor(fitted(hmn), cell=TRUE, weighting="none", norm=2)
> 
> 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(abs(hmu$assoc$phi)^2))))
> stopifnot(all.equal(phin, sqrt(sum(abs(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))))
> 
> stopifnot(all(c(smaorm, smaoru, smaorn) - 1 < 1e-10))
> 
> stopifnot(all.equal(maorm, exp(sqrt(sum(cmaorm)))))
> stopifnot(all.equal(maoru, exp(sqrt(sum(cmaoru)))))
> stopifnot(all.equal(maorn, exp(sqrt(sum(cmaorn)))))
> 
> stopifnot(all.equal(maorm, amaorm))
> stopifnot(all.equal(maoru, amaoru))
> stopifnot(all.equal(maorn, amaorn))
> 
> # 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 <- maor(fitted(hmm), TRUE, weighting="marginal", norm=2, row.weights=w, col.weights=w)
> phiu <- maor(fitted(hmu), TRUE, weighting="uniform", norm=2)
> phin <- maor(fitted(hmn), TRUE, weighting="none", norm=2)
> 
> sphim <- maor(fitted(hmm), TRUE, component="symmetric", weighting="marginal", norm=2,
+               row.weights=w, col.weights=w)
> sphiu <- maor(fitted(hmu), TRUE, component="symmetric", weighting="uniform", norm=2)
> sphin <- maor(fitted(hmn), TRUE, component="symmetric", weighting="none", norm=2)
> 
> aphim <- maor(fitted(hmm), TRUE, component="antisymmetric", weighting="marginal", norm=2,
+               row.weights=w, col.weights=w)
> aphiu <- maor(fitted(hmu), TRUE, component="antisymmetric", weighting="uniform", norm=2)
> aphin <- maor(fitted(hmn), TRUE, component="antisymmetric", weighting="none", norm=2)
> 
> cphim <- maor(fitted(hmm), TRUE, TRUE, weighting="marginal", norm=2,
+               row.weights=w, col.weights=w)
> cphiu <- maor(fitted(hmu), TRUE, TRUE, weighting="uniform", norm=2)
> cphin <- maor(fitted(hmn), TRUE, TRUE, weighting="none", norm=2)
> 
> maorm <- maor(fitted(hmm), weighting="marginal", norm=2,
+               row.weights=w, col.weights=w)
> maoru <- maor(fitted(hmu), weighting="uniform", norm=2)
> maorn <- maor(fitted(hmn), weighting="none", norm=2)
> 
> smaorm <- maor(fitted(hmm), component="symmetric", weighting="marginal", norm=2,
+                row.weights=w, col.weights=w)
> smaoru <- maor(fitted(hmu), component="symmetric", weighting="uniform", norm=2)
> smaorn <- maor(fitted(hmn), component="symmetric", weighting="none", norm=2)
> 
> amaorm <- maor(fitted(hmm), component="antisymmetric", weighting="marginal", norm=2,
+                row.weights=w, col.weights=w)
> amaoru <- maor(fitted(hmu), component="antisymmetric", weighting="uniform", norm=2)
> amaorn <- maor(fitted(hmn), component="antisymmetric", weighting="none", norm=2)
> 
> cmaorm <- maor(fitted(hmm), cell=TRUE, weighting="marginal", norm=2,
+                row.weights=w, col.weights=w)
> cmaoru <- maor(fitted(hmu), cell=TRUE, weighting="uniform", norm=2)
> cmaorn <- maor(fitted(hmn), cell=TRUE, weighting="none", norm=2)
> 
> 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(abs(hmu$assoc$phi)^2))))
> stopifnot(all.equal(sphin, sqrt(sum(abs(hmn$assoc$phi)^2))))
> stopifnot(all.equal(aphim, sqrt(sum((hmm$assoc.hmskew$phi)^2))))
> stopifnot(all.equal(aphiu, sqrt(sum(abs(hmu$assoc.hmskew$phi)^2))))
> stopifnot(all.equal(aphin, sqrt(sum(abs(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))))
> 
> stopifnot(all.equal(maorm, exp(sqrt(sum(cmaorm)))))
> stopifnot(all.equal(maoru, exp(sqrt(sum(cmaoru)))))
> stopifnot(all.equal(maorn, exp(sqrt(sum(cmaorn)))))
> 
> 
> # Test for phi computed from UNIDIFF two-way interaction coefficients
> data(yaish)
> tab <- aperm(yaish[,,-7], 3:1)
> 
> # 1-norm
> # Currently disabled until supported correctly
> # u1m <- unidiff(tab, weighting="marginal", norm=1)
> # rp <- prop.table(margin.table(tab, 1))
> # cp <- prop.table(margin.table(tab, 2))
> # stopifnot(all.equal(u1m$unidiff$phi,
> #                     maor(fitted(u1m)[,,1], TRUE, weighting="marginal", norm=1, rp, cp)))
> # stopifnot(all.equal(u1m$unidiff$phi,
> #                     maor(fitted(u1m)[,,1], TRUE, weighting="marginal", norm=1, rp, cp)))
> # 
> # u1u <- unidiff(tab, weighting="uniform", norm=1)
> # stopifnot(all.equal(u1u$unidiff$phi,
> #                     maor(fitted(u1u)[,,1], TRUE, weighting="uniform", norm=1)))
> # stopifnot(all.equal(u1u$unidiff$phi * exp(u1u$unidiff$layer$qvframe$estimate[2]),
> #                     maor(fitted(u1u)[,,2], TRUE, weighting="uniform", norm=1)))
> # 
> # u1n <- unidiff(tab, weighting="none", norm=1)
> # stopifnot(all.equal(u1n$unidiff$phi,
> #                     maor(fitted(u1n)[,,1], TRUE, weighting="none", norm=1)))
> # stopifnot(all.equal(u1n$unidiff$phi * exp(u1u$unidiff$layer$qvframe$estimate[2]),
> #                     maor(fitted(u1n)[,,2], TRUE, weighting="none", norm=1)))
> 
> 
> # 2-norm
> u2m <- unidiff(tab, weighting="marginal", norm=2)
Initialising
Running start-up iterations..
Running main iterations.........................................................
........
Done
> stopifnot(all.equal(u2m$unidiff$phi,
+                     maor(fitted(u2m), TRUE, weighting="marginal", norm=2), check.attributes=FALSE))
> stopifnot(all.equal(u2m$unidiff$phi[1] * exp(u2m$unidiff$layer$qvframe$estimate),
+                     maor(fitted(u2m), TRUE, weighting="marginal", norm=2), check.attributes=FALSE))
> 
> u2u <- unidiff(tab, weighting="uniform", norm=2)
Initialising
Running start-up iterations..
Running main iterations.........................................................
........
Done
> stopifnot(all.equal(u2u$unidiff$phi,
+                     maor(fitted(u2u), TRUE, weighting="uniform", norm=2), check.attributes=FALSE))
> stopifnot(all.equal(u2u$unidiff$phi[1] * exp(u2u$unidiff$layer$qvframe$estimate),
+                     maor(fitted(u2u), TRUE, weighting="uniform", norm=2), 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,
+                     maor(fitted(u2n), TRUE, weighting="none", norm=2), check.attributes=FALSE))
> stopifnot(all.equal(u2n$unidiff$phi[1] * exp(u2n$unidiff$layer$qvframe$estimate),
+                     maor(fitted(u2n), TRUE, weighting="none", norm=2), check.attributes=FALSE))
> 
> stopifnot(all.equal(maor(yaish),
+                     apply(yaish, 3, maor,
+                           row.weights=margin.table(yaish, 1),
+                           col.weights=margin.table(yaish, 2))))
Warning messages:
1: In maor(yaish) :
  Cells with zero counts found: adding 0.5 to each cell of the table.
2: In FUN(array(newX[, i], d.call, dn.call), ...) :
  Cells with zero counts found: adding 0.5 to each cell of the table.
3: In FUN(array(newX[, i], d.call, dn.call), ...) :
  Cells with zero counts found: adding 0.5 to each cell of the table.
4: In FUN(array(newX[, i], d.call, dn.call), ...) :
  Cells with zero counts found: adding 0.5 to each cell of the table.
5: In FUN(array(newX[, i], d.call, dn.call), ...) :
  More than 25% of cells are empty, the value of the index may not be reliable.
6: In FUN(array(newX[, i], d.call, dn.call), ...) :
  Cells with zero counts found: adding 0.5 to each cell of the table.
7: In FUN(array(newX[, i], d.call, dn.call), ...) :
  Cells with zero counts found: adding 0.5 to each cell of the table.
8: In FUN(array(newX[, i], d.call, dn.call), ...) :
  Cells with zero counts found: adding 0.5 to each cell of the table.
9: In FUN(array(newX[, i], d.call, dn.call), ...) :
  More than 25% of cells are empty, the value of the index may not be reliable.
10: In FUN(array(newX[, i], d.call, dn.call), ...) :
  Cells with zero counts found: adding 0.5 to each cell of the table.
> 
> 
> ###
> ## Comparison with mean/sum of all spanning odds ratios
> ###
> # Can be set to arbitrary values
> nr <- 4
> nc <- 5
> 
> norm <- 2
> 
> 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))
+             if(i2 != i && j2 != j)
+               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))
+             if(i2 != i && j2 != j) {
+               wlor2[i, j, i2, j2] <- log((tab[i, j] * tab[i2, j2]) / (tab[i, j2] * tab[i2, j]))^norm
+               w[i, j, i2, j2] <- rp[i] * cp[j] * rp[i2] * cp[j2]
+             }
+ 
+     wlor2 * w / sum(w, na.rm=TRUE)
+ }
> 
> ## Unweighted mean
> # 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(maor(tab, weighting="uniform", norm=norm),
+       maor(tab, norm=norm, row.weights=rp, col.weights=cp),
+       exp(mean(abs(log(or(tab)))^norm, na.rm=TRUE)^(1/norm)))
+ })
> 
> 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, 1000), 2, 2) + .5
+     rp <- cp <- rep(1/2, 2)
+     c(maor(tab, weighting="uniform", norm=norm),
+       maor(tab, norm=norm, row.weights=rp, col.weights=cp),
+       exp(sqrt(log(tab[1,1] * tab[2,2] / (tab[1,2] * tab[2,1]))^2)))
+ })
> 
> stopifnot(all.equal(res[1,], res[2,]))
> stopifnot(all.equal(res[1,], res[3,]))
> 
> 
> ## Unweighted sum
> # General case
> res <- replicate(10, {
+     tab <- matrix(rpois(nr*nc, 1000), nr, nc) + .5
+     rp <- rep(1, nr)
+     cp <- rep(1, nc)
+     c(maor(tab, weighting="none", norm=norm),
+       exp(sum(abs(log(or(tab)))^norm, na.rm=TRUE)^(1/norm)))
+ })
> 
> stopifnot(all.equal(res[1,], res[2,]))
> 
> # 2x2 table (equality with single odds ratio)
> res <- replicate(10, {
+     tab <- matrix(rpois(2*2, 1000), 2, 2) + .5
+     rp <- cp <- rep(1, 2)
+     c(maor(tab, weighting="none", norm=norm),
+       exp(sqrt(4 * log(tab[1,1] * tab[2,2] / (tab[1,2] * tab[2,1]))^2)))
+ })
> 
> stopifnot(all.equal(res[1,], res[2,]))
> 
> 
> ## Marginal-weighted mean
> # 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(maor(tab, weighting="marginal", norm=norm),
+       maor(tab, norm=norm, row.weights=rp, col.weights=cp),
+       exp(abs(sum(wlor2(tab), na.rm=TRUE))^(1/norm)))
+ })
> 
> 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(maor(tab, weighting="marginal", norm=norm),
+       maor(tab, norm=norm, row.weights=rp, col.weights=cp),
+       exp(abs(log(tab[1,1] * tab[2,2] / (tab[1,2] * tab[2,1])))))
+ })
> 
> stopifnot(all.equal(res[1,], res[2,]))
> stopifnot(all.equal(res[1,], res[3,]))
> 
> res <- replicate(10000, {
+     tab <- matrix(rpois(15, 100), 3, 5)
+     c(log(tab[1,1])+log(tab[2,2])-log(tab[1,2])-log(tab[2,1]),
+       tab[1,1])
+ })
> 
> proc.time()
   user  system elapsed 
   2.76    0.21    2.98