# R script for lecture 7 in STK4900/9900

# Slide 4: Logistic regression -------------------------------------------------

wcgs <- read.table("http://www.uio.no/studier/emner/matnat/math/STK4900/data/wcgs.txt",sep="\t",header=T,na.strings=".")
fit <- glm(chd69~age, data=wcgs, family=binomial)
summary(fit)

a <- predict(fit, type="response", data.frame(age=30))
b <- predict(fit, type="response", data.frame(age=30+1))
(b/(1-b))/(a/(1-a))
exp(fit$coefficients["age"])

a <- predict(fit, type="response", data.frame(age=30))
b <- predict(fit, type="response", data.frame(age=30+10))
(b/(1-b))/(a/(1-a))
exp(fit$coefficients["age"]*10)


# Slide 10: OR confidence interval ----------------------------------------------

or.ci <- function(fit){
  tmp <- summary(fit)$coef
  est <- exp(tmp[,1])
  lower <- exp(tmp[,1]-1.96*tmp[,2])
  upper <- exp(tmp[,1]+1.96*tmp[,2])
  return(cbind(est,lower,upper)) 
}

or.ci(fit)


# Slide 15: Multiple logistic regression ---------------------------------------

fit.mult <- glm(chd69~age+chol+sbp+bmi+smoke, data=wcgs, family=binomial, subset=(chol<600))
summary(fit.mult)

or.ci(fit.mult)


# Slide 17: Multiple logistic regression (cont) --------------------------------

fit.mult.2 <- glm(chd69~age_10+chol_50+sbp_50+bmi_10+smoke, data=wcgs, family=binomial, subset=(chol<600))
summary(fit.mult.2)

or.ci(fit.mult.2)


# Slide 18: Multiple logistic regression (cont) --------------------------------

wcgs$behcat <- factor(wcgs$behpat)
fit.mult.3 <- glm(chd69~age_10+chol_50+sbp_50+bmi_10+smoke+behcat, data=wcgs, family=binomial, subset=(chol<600))
summary(fit.mult.3)


# Slide 28: Model comparison test - Deviance -----------------------------------

anova(fit.mult.2,fit.mult.3,test="Chisq")


# Slide 32: Checking linearity - Grouping variables ----------------------------

fit.catage <- glm(chd69~factor(agec), data=wcgs, family=binomial)
summary(fit.catage)

wcgs$agem <- 39.5*(wcgs$agec==0)+42.9*(wcgs$agec==1)+47.9*(wcgs$agec==2)+52.8*(wcgs$agec==3)+57.3*(wcgs$agec==4)
fit.linage <- glm(chd69~agem, data=wcgs,family=binomial)
summary(fit.linage)

anova(fit.linage, fit.catage,test="Chisq")


# Slide 34: Checking linearity - Square/log transform --------------------------

fit <- glm(chd69~age, data=wcgs,family=binomial)

fit.a2 <- glm(chd69~age+I(age^2), data=wcgs, family=binomial)
anova(fit,fit.a2,test="Chisq")

fit.log <- glm(chd69~age+log(age), data=wcgs,family=binomial)
anova(fit, fit.log, test="Chisq")


# Slide 36: Checking linearity - GAM -------------------------------------------

library(gam)

fit.gam <- gam(chd69~s(age), data=wcgs, family=binomial)
plot(fit.gam, se=TRUE)

anova(fit, fit.gam, test="Chisq")