VisualAcuity: Quasi-independence and quasi-symmetry

Load packages & data

library(vcdExtra)
library(gnm)      # for Diag(), Symm()

data(VisualAcuity, package="vcd")

Work with the data for females

women <- subset(VisualAcuity, gender=="female", select=-gender)

Fit the independence model

indep <- glm(Freq ~ right + left,  data = women, family=poisson)
mosaic(indep, residuals_type="rstandard", gp=shading_Friendly,
       main="Vision data: Independence (women)"  )

Quasi-independence: ignore the diagonal cells by fitting them exactly

quasi.indep <- glm(Freq ~ right + left + Diag(right, left), 
       data = women, family = poisson)
mosaic(quasi.indep, residuals_type="rstandard", gp=shading_Friendly,
       main="Quasi-Independence (women)"  )

Symmetry: test F[i,j] = F[j,i].

Note that this model does not include the ‘main’ effects of right and left, so assumes marginal homogeneity

symmetry <- glm(Freq ~ Symm(right, left), 
       data = women, family = poisson)
mosaic(symmetry, residuals_type="rstandard", gp=shading_Friendly,
       main="Symmetry model (women)"  )

Quasi-symmetry: allow different marginal frequencies for left and right

quasi.symm <- glm(Freq ~ right + left + Symm(right, left), 
       data = women, family = poisson)
mosaic(quasi.symm, residuals_type="rstandard", gp=shading_Friendly,
       main="Quasi-Symmetry model (women)")

model comparisons: for nested models

anova(indep, quasi.indep, quasi.symm, test="Chisq")

## Analysis of Deviance Table
## 
## Model 1: Freq ~ right + left
## Model 2: Freq ~ right + left + Diag(right, left)
## Model 3: Freq ~ right + left + Symm(right, left)
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1         9     6671.5                          
## 2         5      199.1  4   6472.4 < 2.2e-16 ***
## 3         3        7.3  2    191.8 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(symmetry, quasi.symm, test="Chisq")

## Analysis of Deviance Table
## 
## Model 1: Freq ~ Symm(right, left)
## Model 2: Freq ~ right + left + Symm(right, left)
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)   
## 1         6    19.2492                        
## 2         3     7.2708  3   11.978 0.007457 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model summaries, with AIC and BIC

models <- glmlist(indep, quasi.indep, symmetry, quasi.symm)
LRstats(models)

## Likelihood summary table:
##                AIC    BIC LR Chisq Df Pr(>Chisq)    
## indep       6802.9 6808.3   6671.5  9  < 2.2e-16 ***
## quasi.indep  338.5  347.0    199.1  5  < 2.2e-16 ***
## symmetry     156.6  164.4     19.2  6   0.003763 ** 
## quasi.symm   150.6  160.7      7.3  3   0.063751 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1