## Mental health data
library(gnm)
library(vcdExtra)
data(Mental)

# display the frequency table
(Mental.tab <- xtabs(Freq ~ mental+ses, data=Mental))
##           ses
## mental       1   2   3   4   5   6
##   Well      64  57  57  72  36  21
##   Mild      94  94 105 141  97  71
##   Moderate  58  54  65  77  54  54
##   Impaired  46  40  60  94  78  71

fit independence model

# Residual deviance: 47.418 on 15 degrees of freedom
indep <- glm(Freq ~ mental+ses,
                family = poisson, data = Mental)
deviance(indep)
## [1] 47.41785
long.labels <- list(set_varnames = c(mental="Mental Health Status", ses="Parent SES"))
mosaic(indep,residuals_type="rstandard", labeling_args = long.labels, labeling=labeling_residuals,
       main="Mental health data: Independence")

# as a sieve diagram
mosaic(indep, labeling_args = long.labels, panel=sieve, gp=shading_Friendly,
       main="Mental health data: Independence")

fit linear x linear (uniform) association.

Use integer scores for rows/cols

Cscore <- as.numeric(Mental$ses)
Rscore <- as.numeric(Mental$mental)

column effects model (ses)

coleff <- glm(Freq ~ mental + ses + Rscore:ses,
                family = poisson, data = Mental)
mosaic(coleff,residuals_type="rstandard", 
 labeling_args = long.labels, labeling=labeling_residuals, suppress=1, gp=shading_Friendly,
 main="Mental health data: Col effects (ses)")

row effects model (mental)

roweff <- glm(Freq ~ mental + ses + mental:Cscore,
                family = poisson, data = Mental)
mosaic(roweff,residuals_type="rstandard", 
 labeling_args = long.labels, labeling=labeling_residuals, suppress=1, gp=shading_Friendly,
 main="Mental health data: Row effects (mental)")

linear x linear model

linlin <- glm(Freq ~ mental + ses + Rscore:Cscore,
                family = poisson, data = Mental)

# compare models
anova(indep, roweff, coleff, linlin)
## Analysis of Deviance Table
## 
## Model 1: Freq ~ mental + ses
## Model 2: Freq ~ mental + ses + mental:Cscore
## Model 3: Freq ~ mental + ses + Rscore:ses
## Model 4: Freq ~ mental + ses + Rscore:Cscore
##   Resid. Df Resid. Dev Df Deviance
## 1        15     47.418            
## 2        12      6.281  3   41.137
## 3        10      6.829  2   -0.549
## 4        14      9.895 -4   -3.066
AIC(indep, roweff, coleff, linlin)
##        df      AIC
## indep   9 209.5908
## roweff 12 174.4537
## coleff 14 179.0023
## linlin 10 174.0681
mosaic(linlin,residuals_type="rstandard", 
 labeling_args = long.labels, labeling=labeling_residuals, suppress=1, gp=shading_Friendly,
 main="Mental health data: Linear x Linear")

Goodman Row-Column association model fits well (deviance 3.57, df 8)

Mental$mental <- C(Mental$mental, treatment)
Mental$ses <- C(Mental$ses, treatment)
RC1model <- gnm(Freq ~ mental + ses + Mult(mental, ses),
                family = poisson, data = Mental)
## Initialising
## Running start-up iterations..
## Running main iterations........
## Done
mosaic(RC1model, residuals_type="rstandard",
 labeling_args = long.labels, labeling=labeling_residuals, suppress=1, gp=shading_Friendly,
 main="Mental health data: RC(1) model")