Count data models: School absenteeism

library(car)       # for Anova()
library(AER)       # dispersiontest
library(lmtest)    # coeftest()
library(effects)   # effect plots

Load the data

data("quine", package="MASS")

Fixup the data for ease of use

Age should be treated as an ordered factor; use better labels for Lrn

quine$Age <- ordered(quine$Age)
levels(quine$Lrn) <- c("Average", "Slow")

examine sample sizes in table of predictors

quine.tab <- xtabs(~ Lrn + Age + Sex + Eth, data=quine)
ftable(Age + Sex ~ Lrn + Eth, data=quine.tab)

##             Age F0    F1    F2    F3   
##             Sex  F  M  F  M  F  M  F  M
## Lrn     Eth                            
## Average A        4  5  5  2  1  7  9  7
##         N        4  6  6  2  1  7 10  7
## Slow    A        1  3 10  3  8  4  0  0
##         N        1  3 11  7  9  3  0  0

Fit model with all predictors

Age should be treated as an ordered factor

quine$Age <- ordered(quine$Age)
quine.mod1 <- glm(Days ~ ., data=quine, family=poisson)
summary(quine.mod1)

## 
## Call:
## glm(formula = Days ~ ., family = poisson, data = quine)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -6.808  -3.065  -1.119   1.818   9.909  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.80329    0.04343  64.551  < 2e-16 ***
## EthN        -0.53360    0.04188 -12.740  < 2e-16 ***
## SexM         0.16160    0.04253   3.799 0.000145 ***
## Age.L        0.41922    0.04710   8.900  < 2e-16 ***
## Age.Q        0.25188    0.04974   5.064 4.11e-07 ***
## Age.C       -0.30131    0.04103  -7.345 2.07e-13 ***
## LrnSlow      0.34894    0.05204   6.705 2.02e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2073.5  on 145  degrees of freedom
## Residual deviance: 1696.7  on 139  degrees of freedom
## AIC: 2299.2
## 
## Number of Fisher Scoring iterations: 5

Anova(quine.mod1)

## Analysis of Deviance Table (Type II tests)
## 
## Response: Days
##     LR Chisq Df Pr(>Chisq)    
## Eth  166.845  1  < 2.2e-16 ***
## Sex   14.404  1  0.0001475 ***
## Age  168.324  3  < 2.2e-16 ***
## Lrn   45.798  1  1.311e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plot effects

plot(allEffects(quine.mod1), ci.style='bands')

Test for overdispersion

dispersiontest(quine.mod1)

## 
##  Overdispersion test
## 
## data:  quine.mod1
## z = 5.469, p-value = 2.263e-08
## alternative hypothesis: true dispersion is greater than 1
## sample estimates:
## dispersion 
##   12.53013

Re-fit as quasi-poisson

quine.mod1q <- glm(Days ~ ., data=quine, family=quasipoisson)
coeftest(quine.mod1q)

## 
## z test of coefficients:
## 
##             Estimate Std. Error z value  Pr(>|z|)    
## (Intercept)  2.80329    0.15758 17.7893 < 2.2e-16 ***
## EthN        -0.53360    0.15198 -3.5111 0.0004463 ***
## SexM         0.16160    0.15434  1.0470 0.2950968    
## Age.L        0.41922    0.17092  2.4528 0.0141760 *  
## Age.Q        0.25188    0.18050  1.3955 0.1628792    
## Age.C       -0.30131    0.14886 -2.0240 0.0429650 *  
## LrnSlow      0.34894    0.18884  1.8478 0.0646345 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(quine.mod1q)

## Analysis of Deviance Table (Type II tests)
## 
## Response: Days
##     LR Chisq Df Pr(>Chisq)    
## Eth  12.6715  1  0.0003713 ***
## Sex   1.0940  1  0.2955952    
## Age  12.7839  3  0.0051281 ** 
## Lrn   3.4783  1  0.0621799 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Visualize effects

plot(allEffects(quine.mod1q), ci.style="bands")

further analyses: test for interactions

quine.mod2q <- glm(Days ~ .^2, data=quine, family=quasipoisson)

summary(quine.mod2q)

## 
## Call:
## glm(formula = Days ~ .^2, family = quasipoisson, data = quine)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -7.6533  -2.7796  -0.5301   1.5749   8.1955  
## 
## Coefficients: (1 not defined because of singularities)
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    2.861979   0.199381  14.354  < 2e-16 ***
## EthN          -0.767818   0.273599  -2.806  0.00579 ** 
## SexM          -0.163705   0.292751  -0.559  0.57701    
## Age.L         -0.112322   0.277112  -0.405  0.68591    
## Age.Q         -0.011937   0.339923  -0.035  0.97204    
## Age.C         -0.004925   0.369837  -0.013  0.98940    
## LrnSlow        0.654578   0.608579   1.076  0.28414    
## EthN:SexM      0.439024   0.301523   1.456  0.14783    
## EthN:Age.L    -0.165996   0.311575  -0.533  0.59512    
## EthN:Age.Q     1.075223   0.352102   3.054  0.00275 ** 
## EthN:Age.C     0.246600   0.298024   0.827  0.40952    
## EthN:LrnSlow   0.264152   0.372599   0.709  0.47965    
## SexM:Age.L     1.036092   0.316983   3.269  0.00139 ** 
## SexM:Age.Q     0.057859   0.426556   0.136  0.89232    
## SexM:Age.C    -0.515580   0.334144  -1.543  0.12530    
## SexM:LrnSlow   0.041427   0.449197   0.092  0.92666    
## Age.L:LrnSlow  1.126061   1.172251   0.961  0.33857    
## Age.Q:LrnSlow  0.633778   0.921995   0.687  0.49308    
## Age.C:LrnSlow        NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 10.72308)
## 
##     Null deviance: 2073.5  on 145  degrees of freedom
## Residual deviance: 1368.7  on 128  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5

car::Anova(quine.mod2q)

## Analysis of Deviance Table (Type II tests)
## 
## Response: Days
##         LR Chisq Df Pr(>Chisq)    
## Eth      13.6944  1  0.0002151 ***
## Sex       1.7636  1  0.1841801    
## Age      15.5242  3  0.0014193 ** 
## Lrn       6.5018  1  0.0107768 *  
## Eth:Sex   2.1371  1  0.1437742    
## Eth:Age  12.0090  3  0.0073524 ** 
## Eth:Lrn   0.5091  1  0.4755290    
## Sex:Age  13.9482  3  0.0029765 ** 
## Sex:Lrn   0.0085  1  0.9265531    
## Age:Lrn   1.0828  2  0.5819202    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Try adding interactions

quine.mod3 <- update(quine.mod1q, . ~ . + Eth:Age)
plot(allEffects(quine.mod3))

quine.mod4 <- update(quine.mod3, . ~ . + Sex:Age)
plot(allEffects(quine.mod4))

Compare models

anova(quine.mod1q, quine.mod3, quine.mod4, quine.mod2q, test="Chisq")

## Analysis of Deviance Table
## 
## Model 1: Days ~ Eth + Sex + Age + Lrn
## Model 2: Days ~ Eth + Sex + Age + Lrn + Eth:Age
## Model 3: Days ~ Eth + Sex + Age + Lrn + Eth:Age + Sex:Age
## Model 4: Days ~ (Eth + Sex + Age + Lrn)^2
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)   
## 1       139     1696.7                        
## 2       136     1542.8  3  153.870 0.002466 **
## 3       133     1410.0  3  132.782 0.006180 **
## 4       128     1368.7  5   41.386 0.569820   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1