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) # 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) # same, using update() arth.mod1 <- update(arth.mod0, . ~ . + Sex + Treatment) Anova(arth.mod1) # 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 arth.mod2 <- step(arth.mod1, direction="forward", scope=.~ (Age+Sex+Treatment)^2 + Age^2) Anova(arth.mod2) plot(allEffects(arth.mod2))