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)