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 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 209.59 220.19   47.418 15  3.155e-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.11583182 0.11066557  0.06118469 0.31909827  0.22696540
##   Mild:Moderate     -0.07145896 0.07473766 -0.12538139 0.01922754  0.31203110
##   Moderate:Impaired -0.06830298 0.22006188  0.27953207 0.16823542 -0.09402895
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.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")