Example: Loglinear model for contigency tables

This example studies the relationship of taking mariguana, alcohol and cigarette among high school seniors (Agresti book chapter 7.2.6).

Drugs <- read.table("Drugs.dat", header = T)
Drugs

A: Alcohol use C: Cigarette Use M: Marijuana Use

mutual.indep <- glm(count ~ A + C + M, family = poisson(link = log), data = Drugs)
homo.assoc <- update(mutual.indep, .~. + A:C + A:M + C:M)
summary(homo.assoc)

## 
## Call:
## glm(formula = count ~ A + C + M + A:C + A:M + C:M, family = poisson(link = log), 
##     data = Drugs)
## 
## Deviance Residuals: 
##        1         2         3         4         5         6         7         8  
##  0.02044  -0.02658  -0.09256   0.02890  -0.33428   0.09452   0.49134  -0.03690  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  5.63342    0.05970  94.361  < 2e-16 ***
## Ayes         0.48772    0.07577   6.437 1.22e-10 ***
## Cyes        -1.88667    0.16270 -11.596  < 2e-16 ***
## Myes        -5.30904    0.47520 -11.172  < 2e-16 ***
## Ayes:Cyes    2.05453    0.17406  11.803  < 2e-16 ***
## Ayes:Myes    2.98601    0.46468   6.426 1.31e-10 ***
## Cyes:Myes    2.84789    0.16384  17.382  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2851.46098  on 7  degrees of freedom
## Residual deviance:    0.37399  on 1  degrees of freedom
## AIC: 63.417
## 
## Number of Fisher Scoring iterations: 4

The residual devience is \(0.37\) with \(df = 1\). The Chi-squire p-value is \(0.54\), so no three-term interactions are needed.
Interaction terms can be interpreted as odds ratios.

We can also check the residuals of different models.

pres.mutual <- residuals(mutual.indep, type = "pearson")
pres.homo <- residuals(homo.assoc, type = "pearson")
## standardized residuals
stres.homo <- rstandard(homo.assoc, type = "pearson")
fit.mutual <- fitted(mutual.indep)
fit.homo <- fitted(homo.assoc)
signif(cbind(Drugs$count, fit.mutual, pres.mutual, fit.homo, pres.homo, stres.homo), 3)

##       fit.mutual pres.mutual fit.homo pres.homo stres.homo
## 1 911      540.0       16.00   910.00    0.0204      0.633
## 2 538      740.0       -7.43   539.00   -0.0266     -0.633
## 3  44      282.0      -14.20    44.60   -0.0923     -0.633
## 4 456      387.0        3.52   455.00    0.0289      0.633
## 5   3       90.6       -9.20     3.62   -0.3240     -0.633
## 6  43      124.0       -7.29    42.40    0.0947      0.633
## 7   2       47.3       -6.59     1.38    0.5240      0.633
## 8 279       64.9       26.60   280.00   -0.0369     -0.633

We can use a logistic regression to study how taking marijuana associate with taking alcohols and cigarettes. We need to reorganize the data frame to a different form.

Drugs2 <- read.table("Drugs2.dat")
Drugs2

fit.logistic <- glm(M_yes/n ~ A + C, weights = n, family = binomial(link = logit), data= Drugs2)
#out <- cbind(Drugs2$M_yes, Drugs2$M_no)
#fit.logistic <- glm(out ~ A + C, family = binomial(link = logit), data= Drugs2)
summary(fit.logistic)

## 
## Call:
## glm(formula = M_yes/n ~ A + C, family = binomial(link = logit), 
##     data = Drugs2, weights = n)
## 
## Deviance Residuals: 
##        1         2         3         4  
##  0.03353  -0.09697  -0.34739   0.49273  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -5.3090     0.4752 -11.172  < 2e-16 ***
## Ayes          2.9860     0.4647   6.426 1.31e-10 ***
## Cyes          2.8479     0.1638  17.382  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 843.82664  on 3  degrees of freedom
## Residual deviance:   0.37399  on 1  degrees of freedom
## AIC: 25.101
## 
## Number of Fisher Scoring iterations: 4

The intercept is the main effect of \(M\) in the log-linear model.
The coefficients of \(A\) and \(C\) are the same as the interaction terms \(A:M\) and \(C:M\) in the log-linear model.
The residual deviance is the same as in the log-linear model.
The null model here is the model that \(M\) is independent from \((A, C)\) (see below).

joint.indep <- glm(count ~ A + C + M + A:C, family = poisson(link = log), data = Drugs)
anova(joint.indep, homo.assoc)