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. > ## 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 1.29 0.25 1.53