Mental health data

library(gnm)
library(vcdExtra)
library(corrplot)
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 the independence model

Residual deviance: 47.418 on 15 degrees of freedom

indep <- glm(Freq ~ mental+ses,
                family = poisson, data = Mental)
vcdExtra::LRstats(indep)
## Likelihood summary table:
##       AIC BIC LR Chisq Df Pr(>Chisq)    
## indep 210 220     47.4 15    3.2e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Local odds ratios

(LMT <- loddsratio(t(mental.tab)))
## log odds ratios for mental and ses 
## 
##                    ses
## mental                  1:2    2:3     3:4    4:5    5:6
##   Well:Mild          0.1158 0.1107  0.0612 0.3191  0.227
##   Mild:Moderate     -0.0715 0.0747 -0.1254 0.0192  0.312
##   Moderate:Impaired -0.0683 0.2201  0.2795 0.1682 -0.094
corrplot(as.matrix(LMT), method="square", is.corr = FALSE,
         tl.col = "black", tl.srt = 0, tl.offset=1)

Mosaic plot

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.4            
## 2        12        6.3  3     41.1
## 3        10        6.8  2     -0.5
## 4        14        9.9 -4     -3.1
AIC(indep, roweff, coleff, linlin)
##        df AIC
## indep   9 210
## roweff 12 174
## coleff 14 179
## linlin 10 174
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")