Running 'Becker-Clogg1989.R' [166s]
  Running 'Clogg-Shihadeh1994.R' [162s]
  Running 'Goodman1979.R' [2s]
  Running 'Goodman1991.R' [162s]
  Running 'Wong2010-2.7.R' [160s]
  Running 'Wong2010-4.7.R' [160s]
  Running 'Yamaguchi1990.R' [2s]
  Running 'hmskewL.R' [4s]
  Running 'iac.R' [6s]
  Running 'maor.R' [157s]
  Running 'ras.R' [2s]
  Running 'supplementary.R' [168s]
  Running 'vanderHeijden-Mooijaart1995.R' [2s]
Running the tests in 'tests/Becker-Clogg1989.R' failed.
Complete output:
  > ## Becker & Clogg (1989), Table 5 (p. 145)
  > # See also ?rc
  > 
  > options(boot.ncpus=2) # For CRAN policies
  > timings <- as.numeric(Sys.getenv("_R_CHECK_TIMINGS_"))
  > notcran <- Sys.getenv("NOT_CRAN")
  > if(notcran != "" || (!is.na(timings) && timings > 60)) {
  + 
  + library(logmult)
  + data(color)
  + 
  + # "Uniform weights" in the authors' terms mean "no weighting" for us,
  + # and "average marginals" means "marginal" with rcL
  + # See ?rc for "marginals"
  + caithness.unweighted <- rc(color[,,1], nd=2, weighting="none",
  +                            se="jackknife", start=NA)
  + caithness.marginal <- rc(color[,,1], nd=2, weighting="marginal",
  +                          se="jackknife", start=NA)
  + aberdeen.unweighted <- rc(color[,,2], nd=2, weighting="none",
  +                           se="jackknife", start=NA)
  + aberdeen.marginal <- rc(color[,,2], nd=2, weighting="marginal",
  +                         se="jackknife", start=NA)
  + 
  + caithness.unweighted
  + caithness.marginal
  + aberdeen.unweighted
  + aberdeen.marginal
  + 
  + se(caithness.unweighted)
  + se(caithness.marginal)
  + se(aberdeen.unweighted)
  + se(aberdeen.marginal)
  + 
  + summary(caithness.unweighted)
  + summary(caithness.marginal)
  + summary(aberdeen.unweighted)
  + summary(aberdeen.marginal)
  + 
  + stopifnot(isTRUE(all.equal(c(round(caithness.unweighted$assoc$phi, d=3)),
  +                            c(3.067, 0.375))))
  + stopifnot(isTRUE(all.equal(c(round(caithness.unweighted$assoc$row, d=3)),
  +                            c(0.411,  0.478, -0.122, -0.767, 0.472,
  +                             -0.030, -0.804,  0.362))))
  + stopifnot(isTRUE(all.equal(c(round(caithness.unweighted$assoc$col, d=3)),
  +                            c(0.574, 0.279,  0.121, -0.261, -0.714,
  +                              0.271, 0.148, -0.838,  0.449, -0.030))))
  + 
  + stopifnot(isTRUE(all.equal(c(round(caithness.marginal$assoc$phi, d=3)),
  +                            c(0.494, 0.110))))
  + stopifnot(isTRUE(all.equal(c(round(caithness.marginal$assoc$row, d=3)),
  +                            c(0.894,  1.045, -0.166, -1.52, 1.246,
  +                              0.271, -1.356,  0.824))))
  + stopifnot(isTRUE(all.equal(c(round(caithness.marginal$assoc$col, d=3)),
  +                            c(c(1.313, 0.425, -0.019, -1.213, -2.567,
  +                                0.782, 0.518, -1.219,  0.945,  0.038)))))
  + 
  + # This is one of the very few articles reporting jackknife standard errors
  + # Our results are very close for the unweighted solution;
  + # but for marginal weighting their errors are probably wrong,
  + # as they are larger than the unweighted ones, while phi are much smaller.
  + stopifnot(isTRUE(all.equal(c(round(se(caithness.unweighted)$phi, d=3)),
  +                            c(0.282, 0.053))))
  + stopifnot(isTRUE(all.equal(c(round(se(aberdeen.unweighted)$phi, d=3)),
  +                            c(0.221, 0.040))))
  + 
  + 
  + 
  + ## Same with rc(M)-L model
  + # See also ?rcL
  + data(color)
  + 
  + # "Uniform weights" in the authors' terms mean "no weighting" for us,
  + # and "average marginals" means "marginal" with rcL
  + # See ?rc for "marginals"
  + unweighted <- rcL(color, nd=2, weighting="none",
  +                   layer.effect="heterogeneous", se="jackknife", start=NA)
  + marginal <- rcL(color, nd=2, weighting="marginal",
  +                 layer.effect="heterogeneous", se="jackknife", start=NA)
  + unweighted
  + marginal
  + 
  + # (our standard errors are much smaller for the marginal-weighted case)
  + summary(unweighted)
  + summary(marginal)
  + 
  + opar <- par(mfrow=c(1, 2))
  + plot(marginal, layer="Caithness", conf.int=0.95)
  + plot(marginal, layer="Aberdeen", conf.int=0.95)
  + par(opar)
  + 
  + stopifnot(isTRUE(all.equal(c(round(unweighted$assoc$phi, d=3)),
  +                            c(3.067,  2.854, 0.375, 0.294))))
  + stopifnot(isTRUE(all.equal(c(round(unweighted$assoc$row, d=3)),
  +                            c(0.411, 0.478, -0.122, -0.767, 0.472,
  +                             -0.030, -0.804, 0.362, 0.492, 0.399, -0.128,
  +                             -0.763, 0.452, -0.053, -0.797, 0.397))))
  + stopifnot(isTRUE(all.equal(c(round(unweighted$assoc$col, d=3)),
  +                            c(0.574, 0.279,  0.121, -0.261, -0.714, 0.271, 0.148, 
  +                             -0.838, 0.449, -0.030,  0.595,  0.225, 0.135, -0.231,
  +                             -0.724, 0.420,  0.094, -0.867,  0.206, 0.147))))
  + 
  + stopifnot(isTRUE(all.equal(c(round(marginal$assoc$phi, d=3)),
  +                            c(0.465, 0.415, 0.111, 0.085))))
  + stopifnot(isTRUE(all.equal(c(round(marginal$assoc$row, d=3)),
  +                            c(0.883, 1.035, -0.188, -1.554,  1.244,  0.268,
  +                             -1.337, 0.863,  1.175,  0.936, -0.246, -1.512,
  +                              1.184, 0.184, -1.317, 0.976))))
  + stopifnot(isTRUE(all.equal(c(round(marginal$assoc$col, d=3)),
  +                            c(1.361, 0.424, -0.051, -1.301, -2.731,  0.803,  0.544,
  +                             -1.175, 0.970,  0.075,  1.398,  0.212, -0.067, -1.260,
  +                             -2.848, 0.811,  0.463, -1.177,  0.942, 1.144))))
  + 
  + # This is one of the very few articles reporting jackknife standard errors
  + # Our results are very close for the unweighted solution;
  + # but for marginal weighting their errors are probably wrong,
  + # as they are larger than the unweighted ones, while phi are much smaller.
  + stopifnot(isTRUE(all.equal(c(round(se(unweighted)$phi, d=3)),
  +                            c(c(0.282, 0.221, 0.053, 0.04)))))
  + 
  + }
  > 
  > proc.time()
     user  system elapsed 
     0.14    0.09    0.21 
Running the tests in 'tests/Clogg-Shihadeh1994.R' failed.
Complete output:
  > ## Clogg & Shihadeh (1994), Tables 3.5a and b (p. 55-61)
  > # See also ?rc
  > 
  > library(logmult)
  Loading required package: gnm
  
  Attaching package: 'logmult'
  
  The following object is masked from 'package:gnm':
  
      se
  
  > data(gss88)
  > model <- rc(gss88, start=NA)
  Running base model to find starting values...
  Running real model...
  Initialising
  Running main iterations.................................
  Done
  > 
  > # Unweighted scores
  > summary(model, weighting="none")
  Call:
  rc(tab = gss88, start = NA)
  
  Deviance Residuals:
      Min       1Q   Median       3Q      Max  
  -2.7608  -0.4975   0.0749   0.4650   1.9374  
  
  Association coefficients:
                                  Normalized  Adjusted
  Dim1                              8.075599        NA
  Dim1:OccupationProfessional-I     0.565514    1.6071
  Dim1:OccupationProfessional-II    0.564018    1.6028
  Dim1:OccupationManagers           0.146153    0.4153
  Dim1:OccupationClerical          -0.013199   -0.0375
  Dim1:OccupationCraftsmen-I       -0.232452   -0.6606
  Dim1:OccupationCraftsmen-II       0.037196    0.1057
  Dim1:OccupationOperatives (NT)   -0.417400   -1.1861
  Dim1:OccupationOperatives        -0.196173   -0.5575
  Dim1:OccupationFarmers           -0.189421   -0.5383
  Dim1:OccupationService           -0.160650   -0.4565
  Dim1:OccupationMissing           -0.103587   -0.2944
  Dim1:Years.of.Schooling0-7       -0.561801   -1.5965
  Dim1:Years.of.Schooling8-11      -0.409544   -1.1638
  Dim1:Years.of.Schooling12        -0.084106   -0.2390
  Dim1:Years.of.Schooling13-15      0.065232    0.1854
  Dim1:Years.of.Schooling16         0.408340    1.1604
  Dim1:Years.of.Schooling17+        0.581879    1.6536
  
  No standard errors available
  (jackknife and bootstrap disabled, or weighting changed since fitting).
  
  Normalization weights: none
  Deviance:              53.47085
  Pearson chi-squared:   52.12034
  Dissimilarity index:   7.179518%
  Residual df:           36
  BIC:                   -200.8605
  AIC:                   -18.52915
  > # Marginally weighted scores
  > summary(model, weighting="marginal")
  Call:
  rc(tab = gss88, start = NA)
  
  Deviance Residuals:
      Min       1Q   Median       3Q      Max  
  -2.7608  -0.4975   0.0749   0.4650   1.9374  
  
  Association coefficients:
                                  Normalized  Adjusted
  Dim1                             0.7320591        NA
  Dim1:OccupationProfessional-I    1.7834155    1.5259
  Dim1:OccupationProfessional-II   1.7783737    1.5216
  Dim1:OccupationManagers          0.3700858    0.3166
  Dim1:OccupationClerical         -0.1669616   -0.1429
  Dim1:OccupationCraftsmen-I      -0.9058900   -0.7751
  Dim1:OccupationCraftsmen-II      0.0028784    0.0025
  Dim1:OccupationOperatives (NT)  -1.5291990   -1.3084
  Dim1:OccupationOperatives       -0.7836197   -0.6705
  Dim1:OccupationFarmers          -0.7608640   -0.6510
  Dim1:OccupationService          -0.6639001   -0.5680
  Dim1:OccupationMissing          -0.4715891   -0.4035
  Dim1:Years.of.Schooling0-7      -1.8890420   -1.6163
  Dim1:Years.of.Schooling8-11     -1.3906757   -1.1899
  Dim1:Years.of.Schooling12       -0.3254488   -0.2785
  Dim1:Years.of.Schooling13-15     0.1633656    0.1398
  Dim1:Years.of.Schooling16        1.2864265    1.1007
  Dim1:Years.of.Schooling17+       1.8544550    1.5867
  
  No standard errors available
  (jackknife and bootstrap disabled, or weighting changed since fitting).
  
  Normalization weights: marginal
  Deviance:              53.47085
  Pearson chi-squared:   52.12034
  Dissimilarity index:   7.179518%
  Residual df:           36
  BIC:                   -200.8605
  AIC:                   -18.52915
  > # Uniformly weighted scores
  > summary(model, weighting="uniform")
  Call:
  rc(tab = gss88, start = NA)
  
  Deviance Residuals:
      Min       1Q   Median       3Q      Max  
  -2.7608  -0.4975   0.0749   0.4650   1.9374  
  
  Association coefficients:
                                  Normalized  Adjusted
  Dim1                              0.994038        NA
  Dim1:OccupationProfessional-I     1.875598    1.8700
  Dim1:OccupationProfessional-II    1.870636    1.8651
  Dim1:OccupationManagers           0.484735    0.4833
  Dim1:OccupationClerical          -0.043775   -0.0436
  Dim1:OccupationCraftsmen-I       -0.770957   -0.7687
  Dim1:OccupationCraftsmen-II       0.123365    0.1230
  Dim1:OccupationOperatives (NT)   -1.384358   -1.3802
  Dim1:OccupationOperatives        -0.650631   -0.6487
  Dim1:OccupationFarmers           -0.628237   -0.6264
  Dim1:OccupationService           -0.532814   -0.5312
  Dim1:OccupationMissing           -0.343560   -0.3425
  Dim1:Years.of.Schooling0-7       -1.376125   -1.3720
  Dim1:Years.of.Schooling8-11      -1.003175   -1.0002
  Dim1:Years.of.Schooling12        -0.206016   -0.2054
  Dim1:Years.of.Schooling13-15      0.159786    0.1593
  Dim1:Years.of.Schooling16         1.000224    0.9972
  Dim1:Years.of.Schooling17+        1.425306    1.4211
  
  No standard errors available
  (jackknife and bootstrap disabled, or weighting changed since fitting).
  
  Normalization weights: uniform
  Deviance:              53.47085
  Pearson chi-squared:   52.12034
  Dissimilarity index:   7.179518%
  Residual df:           36
  BIC:                   -200.8605
  AIC:                   -18.52915
  > 
  > unweighted <- assoc(model, weighting="none")
  > marginal <- assoc(model, weighting="marginal")
  > uniform <- assoc(model, weighting="uniform")
  > 
  > stopifnot(all.equal(round(c(unweighted$row, unweighted$row * sqrt(unweighted$phi[1]),
  +                             marginal$row, marginal$row * sqrt(marginal$phi[1]),
  +                             uniform$row, uniform$row * sqrt(uniform$phi[1])), 2),
  +                     c(0.57, 0.56, 0.15, -0.01, -0.23, 0.04, -0.42, -0.2, -0.19, -0.16, -0.1,
  +                       1.61, 1.6, 0.42, -0.04, -0.66, 0.11, -1.19, -0.56, -0.54, -0.46, -0.29,
  +                       1.78, 1.78, 0.37, -0.17, -0.91, 0, -1.53, -0.78, -0.76, -0.66, -0.47,
  +                       1.53, 1.52, 0.32, -0.14, -0.78, 0, -1.31, -0.67, -0.65, -0.57, -0.4,
  +                       1.88, 1.87, 0.48, -0.04,  -0.77, 0.12, -1.38, -0.65, -0.63, -0.53, -0.34, 
  +                       1.87, 1.87, 0.48, -0.04, -0.77, 0.12, -1.38, -0.65, -0.63, -0.53, -0.34)))
  > stopifnot(all.equal(round(c(unweighted$col, unweighted$col * sqrt(unweighted$phi[1]),
  +                             marginal$col, marginal$col * sqrt(marginal$phi[1]),
  +                             uniform$col, uniform$col * sqrt(uniform$phi[1])), 2),
  +                     c(-0.56, -0.41, -0.08, 0.07, 0.41, 0.58,
  +                       -1.6, -1.16, -0.24, 0.19, 1.16, 1.65,
  +                       -1.89, -1.39, -0.33, 0.16, 1.29, 1.85,
  +                       -1.62, -1.19, -0.28, 0.14, 1.1, 1.59,
  +                       -1.38, -1, -0.21, 0.16, 1, 1.43,
  +                       -1.37, -1, -0.21, 0.16, 1, 1.42)))
  > 
  > # Test anova
  > indep <- gnm(Freq ~ Occupation + Years.of.Schooling, data=gss88, family=poisson)
  > anova(indep, unweighted, test="LR")
  Analysis of Deviance Table
  
  Model: poisson, link: log
  
  Response: Freq
  
  Terms added sequentially (first to last)
  
  
                     Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
  NULL                                  65    1359.84              
  Occupation         10   464.10        55     895.75 < 2.2e-16 ***
  Years.of.Schooling  5   440.69        50     455.06 < 2.2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  > 
  > proc.time()
     user  system elapsed 
     1.31    0.17    1.46 
Running the tests in 'tests/Goodman1991.R' failed.
Complete output:
  > ## Goodman (1991), Table 17.1 (p. 1097)
  > # See also ?rc
  > 
  > library(logmult)
  Loading required package: gnm
  
  Attaching package: 'logmult'
  
  The following object is masked from 'package:gnm':
  
      se
  
  > data(criminal)
  > 
  > model <- rc(criminal, start=NA)
  Running base model to find starting values...
  Running real model...
  Initialising
  Running main iterations.....................
  Done
  > 
  > model$assoc # These are the phi (.07), mu and nu
  $phi
             Dim1
  [1,] 0.07051354
  
  $row
  , , 1
  
             Dim1
  1955  1.2591887
  1956  0.8225979
  1957 -0.5602021
  1958 -1.2074486
  
  
  $col
  , , 1
  
           Dim1
  15 -1.4350272
  16 -1.2965748
  17 -0.3324492
  18  0.9995791
  19  0.8889776
  
  
  $diagonal
  numeric(0)
  
  $weighting
  [1] "marginal"
  
  $row.weights
       [,1]
  1955 1614
  1956 1615
  1957 1889
  1958 1907
  
  $col.weights
     [,1]
  15  693
  16 1381
  17 1485
  18 1786
  19 1680
  
  $vars
  [1] "Year" "Age" 
  
  $covmat
  numeric(0)
  
  $adj.covmats
  numeric(0)
  
  $boot.results
  numeric(0)
  
  $jack.results
  numeric(0)
  
  $covtype
  [1] "none"
  
  attr(,"class")
  [1] "assoc.rc" "assoc"   
  > model$assoc$row[,1,1] * model$assoc$phi[1,1] # These are the mu'
         1955        1956        1957        1958 
   0.08878985  0.05800429 -0.03950183 -0.08514147 
  > model$assoc$col[,1,1] * model$assoc$phi[1,1] # These are the nu'
           15          16          17          18          19 
  -0.10118884 -0.09142608 -0.02344217  0.07048386  0.06268496 
  > 
  > stopifnot(round(model$assoc$phi[1,1], d=2) == 0.07)
  > # Scores have reversed signs compared to the orignal article
  > stopifnot(isTRUE(all.equal(round(model$assoc$row[,1,1], d=2),
  +                            c(1.26, 0.82, -0.56, -1.21),
  +                            check.attributes=FALSE)))
  > stopifnot(isTRUE(all.equal(round(model$assoc$col[,1,1], d=2),
  +                            c(-1.44, -1.30, -0.33, 1.00, 0.89),
  +                            check.attributes=FALSE)))
  > 
  > proc.time()
     user  system elapsed 
     1.17    0.17    1.32 
Running the tests in 'tests/Wong2010-2.7.R' failed.
Complete output:
  > ## Wong (2010), Table 2.7 (p. 48-49)
  > # See also ?rc
  > 
  > options(boot.ncpus=2) # For CRAN policies
  > timings <- as.numeric(Sys.getenv("_R_CHECK_TIMINGS_"))
  > notcran <- Sys.getenv("NOT_CRAN")
  > if(notcran != "" || (!is.na(timings) && timings > 60)) {
  + 
  + library(logmult)
  + data(gss8590)
  + 
  + # The table used in Wong (2010) is not perfectly consistent
  + # with that of Wong (2001)
  + tab <- margin.table(gss8590[,,c(2,4)], 1:2)
  + tab[2,4] <- 49
  + 
  + model <- rc(tab, nd=2, weighting="none", se="jackknife", start=NA)
  + 
  + model
  + summary(model) # Jackknife standard errors are slightly different
  +                # from their asymptotic counterparts
  + 
  + # Compare with bootstrap standard errors
  + model2 <- rc(tab, nd=2, weighting="none", se="bootstrap", start=NA)
  + plot(model, conf.int=0.95)
  + summary(model2)
  + 
  + # A few scores differ from reported results by .001
  + stopifnot(isTRUE(all.equal(round(c(model$assoc$phi), d=3),
  +                            c(2.601, 1.522))))
  + stopifnot(isTRUE(all.equal(round(c(model$assoc$row), d=3),
  +                            c(0.743,  0.088, -0.200, -0.632,
  +                             0.276, -0.061, -0.777,  0.562))))
  + stopifnot(isTRUE(all.equal(round(c(model$assoc$col), d=3),
  +                            c(0.765,  0.020, -0.322, -0.55, 0.088,
  +                             -0.137, -0.549, -0.120, -0.01, 0.816))))
  + 
  + }
  > 
  > proc.time()
     user  system elapsed 
     0.17    0.03    0.18 
Running the tests in 'tests/Wong2010-4.7.R' failed.
Complete output:
  > ## Wong (2010), Table 4.7 (p. 103), model 9
  > 
  > options(boot.ncpus=2) # For CRAN policies
  > timings <- as.numeric(Sys.getenv("_R_CHECK_TIMINGS_"))
  > notcran <- Sys.getenv("NOT_CRAN")
  > if(notcran != "" || (!is.na(timings) && timings > 60)) {
  + 
  + library(logmult)
  + data(gss7590)
  + 
  + set.seed(1)
  + model <- rcL(gss7590, nd=2, weighting="none", se="jackknife")
  + 
  + model
  + summary(model) # Jackknife standard errors are slightly different
  +                # from their asymptotic counterparts
  + 
  + stopifnot(all.equal(round(c(model$assoc$phi), d=3),
  +                     c(3.075, 3.474,  2.949, 2.460,
  +                       0.539, 1.686, -0.770, 1.891)))
  + stopifnot(all.equal(round(c(model$assoc$row[,,1]), d=3),
  +                     c(0.640,  0.239, -0.168, -0.711,
  +                       0.731, -0.217, -0.636,  0.121)))
  + stopifnot(all.equal(round(c(model$assoc$col[,,1]), d=3),
  +                     c(0.765,  0.250, -0.398, -0.273, -0.344,
  +                       0.071, -0.480, -0.198, -0.216,  0.824)))
  + 
  + # Check scores of heterogeneous model (scores not given by Wong)
  + model.heterog <- rcL(gss7590, nd=2, layer.effect="heterogeneous",
  +                      weighting="none")
  + sep.models <- lapply(1:4, function(i) rc(gss7590[,,i], nd=2,
  +                                          weighting="none"))
  + 
  + # Level 2 does not converge properly because there are too few farmers:
  + # values are always slightly different and cannot be checked
  + stopifnot(isTRUE(all.equal(model.heterog$assoc$phi[-2,],
  +                            rbind(sep.models[[1]]$assoc$phi,
  +                                  sep.models[[3]]$assoc$phi,
  +                                  sep.models[[4]]$assoc$phi),
  +                            check.attributes=FALSE,
  +                            tolerance=1e-6)))
  + stopifnot(isTRUE(all.equal(model.heterog$assoc$row[,,-2],
  +                            array(c(sep.models[[1]]$assoc$row,
  +                                  sep.models[[3]]$assoc$row,
  +                                  sep.models[[4]]$assoc$row),
  +                                  dim=c(nrow(gss7590), 2, 3)),
  +                            check.attributes=FALSE,
  +                            tolerance=1e-6)))
  + stopifnot(isTRUE(all.equal(model.heterog$assoc$col[,,-2],
  +                            array(c(sep.models[[1]]$assoc$col,
  +                                  sep.models[[3]]$assoc$col,
  +                                  sep.models[[4]]$assoc$col),
  +                                  dim=c(ncol(gss7590), 2, 3)),
  +                            check.attributes=FALSE,
  +                            tolerance=1e-6)))
  + 
  + # Test anova
  + indep <- gnm(Freq ~ Education*Group + Occupation*Group, data=gss7590, family=poisson)
  + a <- anova(indep, model, model.heterog, test="LR")
  + 
  + }
  > 
  > proc.time()
     user  system elapsed 
     0.18    0.06    0.23 
Running the tests in 'tests/maor.R' failed.
Complete output:
  > # 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 
     3.46    0.48    3.93 
Running the tests in 'tests/supplementary.R' failed.
Complete output:
  > ## Variations based on Becker & Clogg (1989), cf. Becker-Clogg1989.R
  > 
  > library(logmult)
  Loading required package: gnm
  
  Attaching package: 'logmult'
  
  The following object is masked from 'package:gnm':
  
      se
  
  > data(color)
  > 
  > caithness <- color[,,"Caithness"]
  > 
  > 
  > ## Adding the table again as supplementary rows
  > rsup <- caithness
  > csup <- caithness
  > rownames(rsup) <- paste(rownames(caithness), "2")
  > colnames(csup) <- paste(colnames(caithness), "2")
  > 
  > names(dimnames(rsup)) <- names(dimnames(csup)) <- names(dimnames(caithness))
  > 
  > mA1r <- rc(caithness, 1, rowsup=rsup, weighting="u", start=NA)
  Running base model to find starting values...
  Running real model...
  Initialising
  Running main iterations.............................
  Done
  > mA1c <- rc(caithness, 1, colsup=csup, weighting="u", start=NA)
  Running base model to find starting values...
  Running real model...
  Initialising
  Running main iterations.............................
  Done
  > mA1 <- rc(caithness, 1, rowsup=rsup, colsup=csup, weighting="u", start=NA)
  Running base model to find starting values...
  Running real model...
  Initialising
  Running main iterations.............................
  Done
  > 
  > rscoresA1 <- mA1$assoc$row[1:4, 1, 1] * sqrt(mA1$assoc$phi[1, 1])
  > cscoresA1 <- mA1$assoc$col[1:5, 1, 1] * sqrt(mA1$assoc$phi[1, 1])
  > 
  > rmA1 <- gnm(Freq ~ Eye + Hair + Eye:cscoresA1[col(rsup)],
  +                data=as.table(rsup), family=poisson)
  > cmA1 <- gnm(Freq ~ Eye + Hair + Hair:rscoresA1[row(csup)],
  +                data=as.table(csup), family=poisson)
  > 
  > rscA1 <- getContrasts(rmA1, pickCoef(rmA1, "cscores"), ref="mean")$qvframe[,1]
  > cscA1 <- getContrasts(cmA1, pickCoef(cmA1, "rscores"), ref="mean")$qvframe[,1]
  > 
  > stopifnot(all.equal(rscA1, rscoresA1, tolerance=1e-6, check.attributes=FALSE))
  > stopifnot(all.equal(cscA1, cscoresA1, tolerance=1e-6, check.attributes=FALSE))
  > stopifnot(all.equal(rscA1, mA1$assoc$row[5:8,,] * sqrt(mA1$assoc$phi[1, 1]),
  +                     tolerance=1e-6, check.attributes=FALSE))
  > stopifnot(all.equal(cscA1, mA1$assoc$col[6:10,,] * sqrt(mA1$assoc$phi[1, 1]),
  +                     tolerance=1e-6, check.attributes=FALSE))
  > 
  > 
  > ## Merging some categories and re-adding existing ones
  > rsup <- rbind(Fairer=colSums(caithness[1:2,]),
  +               Darker=colSums(caithness[3:4,]),
  +               Dark2=caithness[4,])
  > 
  > csup <- cbind(FairRed=rowSums(caithness[,1:2]),
  +               Darker=rowSums(caithness[,3:5]),
  +               Black2=caithness[,5])
  > 
  > names(dimnames(rsup)) <- names(dimnames(csup)) <- names(dimnames(caithness))
  > 
  > # With one dimension
  > mB1r <- rc(caithness, 1, rowsup=rsup, weighting="u", start=NA)
  Running base model to find starting values...
  Running real model...
  Initialising
  Running main iterations.............................
  Done
  > mB1c <- rc(caithness, 1, colsup=csup, weighting="u", start=NA)
  Running base model to find starting values...
  Running real model...
  Initialising
  Running main iterations.............................
  Done
  > mB1 <- rc(caithness, 1, rowsup=rsup, colsup=csup, weighting="u", start=NA)
  Running base model to find starting values...
  Running real model...
  Initialising
  Running main iterations.............................
  Done
  > 
  > rscoresB1 <- mB1$assoc$row[1:4, 1, 1] * sqrt(mB1$assoc$phi[1, 1])
  > cscoresB1 <- mB1$assoc$col[1:5, 1, 1] * sqrt(mB1$assoc$phi[1, 1])
  > 
  > rmB1 <- gnm(Freq ~ Eye + Hair + Eye:cscoresB1[col(rsup)],
  +             data=as.table(rsup), family=poisson)
  > cmB1 <- gnm(Freq ~ Eye + Hair + Hair:rscoresB1[row(csup)],
  +             data=as.table(csup), family=poisson)
  > 
  > rscB1 <- getContrasts(rmB1, pickCoef(rmB1, "cscores"), ref="mean")$qvframe[,1]
  > cscB1 <- getContrasts(cmB1, pickCoef(cmB1, "rscores"), ref="mean")$qvframe[,1]
  > 
  > stopifnot(all.equal(rscB1, mB1$assoc$row[5:7,,] * sqrt(mB1$assoc$phi[1, 1]),
  +                     tolerance=1e-6, check.attributes=FALSE))
  > stopifnot(all.equal(cscB1, mB1$assoc$col[6:8,,] * sqrt(mB1$assoc$phi[1, 1]),
  +                     tolerance=1e-6, check.attributes=FALSE))
  > 
  > # With two dimensions
  > mB2r <- rc(caithness, 2, rowsup=rsup, weighting="u", start=NA)
  Running base model to find starting values...
  Running real model...
  Initialising
  Running main iterations......................
  Done
  > mB2c <- rc(caithness, 2, colsup=csup, weighting="u", start=NA)
  Running base model to find starting values...
  Running real model...
  Initialising
  Running main iterations......................
  Done
  > mB2 <- rc(caithness, 2, rowsup=rsup, colsup=csup, weighting="u", start=NA)
  Running base model to find starting values...
  Running real model...
  Initialising
  Running main iterations......................
  Done
  > 
  > rscoresB2 <- sweep(mB2$assoc$row[1:4,, 1], 2, sqrt(mB2$assoc$phi[1,]), "*")
  > cscoresB2 <- sweep(mB2$assoc$col[1:5,, 1], 2, sqrt(mB2$assoc$phi[1,]), "*")
  > 
  > rmB2 <- gnm(Freq ~ Eye + Hair + Eye:cscoresB2[col(rsup), 1] + Eye:cscoresB2[col(rsup), 2],
  +             data=as.table(rsup), family=poisson)
  > cmB2 <- gnm(Freq ~ Eye + Hair + Hair:rscoresB2[row(csup), 1] + Hair:rscoresB2[row(csup), 2],
  +             data=as.table(csup), family=poisson)
  > 
  > rscB2 <- cbind(getContrasts(rmB2, pickCoef(rmB2, "cscores.*1\\]"), ref="mean")$qvframe[,1],
  +                getContrasts(rmB2, pickCoef(rmB2, "cscores.*2\\]"), ref="mean")$qvframe[,1])
  > cscB2 <- cbind(getContrasts(cmB2, pickCoef(cmB2, "rscores.*1\\]"), ref="mean")$qvframe[,1],
  +                getContrasts(cmB2, pickCoef(cmB2, "rscores.*2\\]"), ref="mean")$qvframe[,1])
  > 
  > stopifnot(all.equal(rscB2, sweep(mB2$assoc$row[5:7,,], 2, sqrt(mB2$assoc$phi[1,]), "*"),
  +                     tolerance=1e-6, check.attributes=FALSE))
  > stopifnot(all.equal(cscB2, sweep(mB2$assoc$col[6:8,,], 2, sqrt(mB2$assoc$phi[1,]), "*"),
  +                     tolerance=1e-6, check.attributes=FALSE))
  > 
  > # With two dimensions and marginal weighting
  > mB3r <- rc(caithness, 2, rowsup=rsup, weighting="m", start=NA)
  Running base model to find starting values...
  Running real model...
  Initialising
  Running main iterations......................
  Done
  > mB3c <- rc(caithness, 2, colsup=csup, weighting="m", start=NA)
  Running base model to find starting values...
  Running real model...
  Initialising
  Running main iterations......................
  Done
  > mB3 <- rc(caithness, 2, rowsup=rsup, colsup=csup, weighting="m", start=NA)
  Running base model to find starting values...
  Running real model...
  Initialising
  Running main iterations......................
  Done
  > 
  > rscoresB3 <- sweep(mB3$assoc$row[1:4,, 1], 2, sqrt(mB3$assoc$phi[1,]), "*")
  > cscoresB3 <- sweep(mB3$assoc$col[1:5,, 1], 2, sqrt(mB3$assoc$phi[1,]), "*")
  > 
  > rmB3 <- gnm(Freq ~ Eye + Hair + Eye:cscoresB3[col(rsup), 1] + Eye:cscoresB3[col(rsup), 2],
  +             data=as.table(rsup), family=poisson)
  > cmB3 <- gnm(Freq ~ Eye + Hair + Hair:rscoresB3[row(csup), 1] + Hair:rscoresB3[row(csup), 2],
  +             data=as.table(csup), family=poisson)
  > 
  > rscB3 <- cbind(getContrasts(rmB3, pickCoef(rmB3, "cscores.*1\\]"),
  +                             ref=rowSums(rsup)/sum(rsup))$qvframe[,1],
  +                getContrasts(rmB3, pickCoef(rmB3, "cscores.*2\\]"),
  +                             ref=rowSums(rsup)/sum(rsup))$qvframe[,1])
  > cscB3 <- cbind(getContrasts(cmB3, pickCoef(cmB3, "rscores.*1\\]"),
  +                             ref=colSums(csup)/sum(csup))$qvframe[,1],
  +                getContrasts(cmB3, pickCoef(cmB3, "rscores.*2\\]"),
  +                             ref=colSums(csup)/sum(csup))$qvframe[,1])
  > 
  > stopifnot(all.equal(rscB3, sweep(mB3$assoc$row[5:7,,], 2, sqrt(mB3$assoc$phi[1,]), "*"),
  +                     tolerance=1e-6, check.attributes=FALSE))
  > stopifnot(all.equal(cscB3, sweep(mB3$assoc$col[6:8,,], 2, sqrt(mB3$assoc$phi[1,]), "*"),
  +                     tolerance=1e-6, check.attributes=FALSE))
  > 
  > 
  > proc.time()
     user  system elapsed 
     2.15    0.25    2.37