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.45
## Age   1   7.2852        82     109.16

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.1587  1  0.0130764 *  
## Sex         6.9829  1  0.0082293 ** 
## Treatment  12.1965  1  0.0004788 ***
## ---
## 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.1587  1  0.0130764 *  
## Sex         6.9829  1  0.0082293 ** 
## Treatment  12.1965  1  0.0004788 ***
## ---
## 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.06
## Better ~ Age + Sex + Treatment
## 
##                 Df Deviance    AIC
## + Age:Sex        1   88.640  98.64
## <none>               92.063 100.06
## + Age:Treatment  1   91.452 101.45
## + Sex:Treatment  1   91.636 101.64
## 
## Step:  AIC=98.64
## Better ~ Age + Sex + Treatment + Age:Sex
## 
##                 Df Deviance    AIC
## <none>               88.640  98.64
## + Sex:Treatment  1   88.369 100.37
## + Age:Treatment  1   88.631 100.63

## 
## Call:  glm(formula = Better ~ Age + Sex + Treatment + Age:Sex, family = "binomial", 
##     data = Arthritis)
## 
## Coefficients:
##      (Intercept)               Age           SexMale  TreatmentTreated       Age:SexMale  
##         -4.55477           0.07734           2.75066           1.79705          -0.07945  
## 
## Degrees of Freedom: 83 Total (i.e. Null);  79 Residual
## Null Deviance:       116.4 
## Residual Deviance: 88.64     AIC: 98.64