Nothing
# set seed to fix sign of coef
suppressWarnings(RNGversion("3.0.0"))
set.seed(1)
# Agresti A (2002).Categorical Data Analysis. 2nd edition
mentalHealth$MHS <- C(mentalHealth$MHS, treatment)
mentalHealth$SES <- C(mentalHealth$SES, treatment)
RC1model <- gnm(count ~ SES + MHS +
Mult(-1 + SES, -1 + MHS),
family = poisson, data = mentalHealth, verbose = FALSE)
test_that("RC model as expected for mentalHealth data", {
# compare vs results in sec 9.6.2
pearson_chi_sq <- sum(na.omit(c(residuals(RC1model, type = "pearson")))^2)
expect_equal(round(pearson_chi_sq, 1), 3.6)
expect_equal(df.residual(RC1model), 8)
# normalize as in Agresti's eqn 9.15
rowProbs <- with(mentalHealth, tapply(count, SES, sum) / sum(count))
colProbs <- with(mentalHealth, tapply(count, MHS, sum) / sum(count))
mu <- getContrasts(RC1model, pickCoef(RC1model, "[.]SES"),
ref = rowProbs, scaleRef = rowProbs,
scaleWeights = rowProbs)
nu <- getContrasts(RC1model, pickCoef(RC1model, "[.]MHS"),
ref = colProbs, scaleRef = colProbs,
scaleWeights = colProbs)
# change of scale
expect_equal(round(-mu$qvframe$Estimate, 2),
c(-1.11, -1.12, -0.37, 0.03, 1.01, 1.82))
expect_equal(round(-nu$qvframe$Estimate, 2),
c(-1.68, -0.14, 0.14, 1.41))
# association parameter
rowScores <- coef(RC1model)[10:15]
colScores <- coef(RC1model)[16:19]
rowScores <- rowScores - sum(rowScores * rowProbs)
colScores <- colScores - sum(colScores * colProbs)
beta1 <- sqrt(sum(rowScores^2 * rowProbs))
beta2 <- sqrt(sum(colScores^2 * colProbs))
expect_equal(round(beta1 * beta2, 2), 0.17)
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.