Output from arthritis-logistic.R

library(vcd)
library(car)
data(Arthritis)

# define Better
Arthritis$Better <- Arthritis$Improved > "None"

# simple linear regression
arth.mod0 <- glm(Better ~ Age, data = Arthritis, family = "binomial")
anova(arth.mod0)

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Better
## 
## Terms added sequentially (first to last)
## 
## 
##      Df Deviance Resid. Df Resid. Dev
## NULL                    83        116
## Age   1     7.29        82        109

# plot, with +-1 SE
plot(Better ~ Age, data = Arthritis, ylab = "Prob (Better)")
pred <- predict(arth.mod0, type = "response", se.fit = TRUE)
ord <- order(Arthritis$Age)

lines(Arthritis$Age[ord], pred$fit[ord], lwd = 3)
upper <- pred$fit + pred$se.fit
lower <- pred$fit - pred$se.fit
lines(Arthritis$Age[ord], upper[ord], lty = 2, col = "blue")
lines(Arthritis$Age[ord], lower[ord], lty = 2, col = "blue")
# smoothed non-parametric curve
lines(lowess(Arthritis$Age[ord], Arthritis$Better[ord]), lwd = 2, col = "red")

# main effects model
arth.mod1 <- glm(Better ~ Age + Sex + Treatment, data = Arthritis, family = "binomial")
Anova(arth.mod1)

## Analysis of Deviance Table (Type II tests)
## 
## Response: Better
##           LR Chisq Df Pr(>Chisq)    
## Age           6.16  1    0.01308 *  
## Sex           6.98  1    0.00823 ** 
## Treatment    12.20  1    0.00048 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# same, using update()
arth.mod1 <- update(arth.mod0, . ~ . + Sex + Treatment)
Anova(arth.mod1)

## Analysis of Deviance Table (Type II tests)
## 
## Response: Better
##           LR Chisq Df Pr(>Chisq)    
## Age           6.16  1    0.01308 *  
## Sex           6.98  1    0.00823 ** 
## Treatment    12.20  1    0.00048 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# plots, using effects package
library(effects)
arth.eff <- allEffects(arth.mod1, xlevels = list(Age = seq(25, 75, 5)))
plot(arth.eff, ylab = "Prob(Better)")

# forward selection from the main effects model
step(arth.mod1, direction = "forward", scope = . ~ (Age + Sex + Treatment)^2 + Age^2)

## Start:  AIC=100.1
## Better ~ Age + Sex + Treatment
## 
##                 Df Deviance   AIC
## + Age:Sex        1     88.6  98.6
##                  92.1 100.1
## + Age:Treatment  1     91.5 101.5
## + Sex:Treatment  1     91.6 101.6
## 
## Step:  AIC=98.64
## Better ~ Age + Sex + Treatment + Age:Sex
## 
##                 Df Deviance   AIC
##                  88.6  98.6
## + Sex:Treatment  1     88.4 100.4
## + Age:Treatment  1     88.6 100.6

## 
## Call:  glm(formula = Better ~ Age + Sex + Treatment + Age:Sex, family = "binomial", 
##     data = Arthritis)
## 
## Coefficients:
##      (Intercept)               Age           SexMale  TreatmentTreated       Age:SexMale  
##          -4.5548            0.0773            2.7507            1.7970           -0.0794  
## 
## Degrees of Freedom: 83 Total (i.e. Null);  79 Residual
## Null Deviance:	    116 
## Residual Deviance: 88.6 	AIC: 98.6