Arthritis data: Logistic regression

define Better from ordinal Improved

Arthritis$Better <- Arthritis$Improved > 'None'

simple linear logistic 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
## Better ~ Age + Sex + Treatment
## 
##                 Df Deviance   AIC
## + Age:Sex        1     88.6  98.6
## <none>                 92.1 100.1
## + Age:Treatment  1     91.5 101.5
## + Sex:Treatment  1     91.6 101.6
## 
## Step:  AIC=98.6
## Better ~ Age + Sex + Treatment + Age:Sex
## 
##                 Df Deviance   AIC
## <none>                 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