source("http://192.38.117.59/~linearpredictors/datafiles/readSurgery.R")

## remove surgery type "orthopedic"
surgery2 <- subset(surgery, surgtype %in% 2:3)
surgery2$surgtype <- surgery2$surgtype[drop = TRUE]

logistic <- glm(complication ~ surgtype,
                family = binomial,
                data = surgery2)

## Confidence limits for discrepancies between complication frequencies

## OR, abdominal vs gynecolocical

exp(logistic$coef[2] + 1.96*summary(logistic)$coef[2,2])
exp(logistic$coef[2] + 1.96*summary(logistic)$coef[2,2])

## Risk ratio

Log <- glm(complication ~ surgtype,
           family = binomial(link = "log"),
           data = surgery2)

list(
     RR = exp(Log$coef[2]),
     lower = exp(Log$coef[2] - 1.96*summary(Log)$coef[2, "Std. Error"]),
     upper = exp(Log$coef[2] + 1.96*summary(Log)$coef[2, "Std. Error"])
     )

## Risk difference (note that R does not willingly use the identity link for
## the binomial family)

Id <- glm(complication ~ surgtype,
          family = binomial(link = make.link("identity")),
          data = surgery2)

list(
     diff = Id$coef[2],
     lower = Id$coef[2] - 1.96*summary(Id)$coef[2, "Std. Error"],
     upper = Id$coef[2] + 1.96*summary(Id)$coef[2, "Std. Error"]
     )
## roughly the same as
confint.default(Id)

## other methods for obtaining the risk difference:
p <- predict(logistic,
             ## predicting the probability for both surgery types
             data.frame(surgtype = factor(2:3)),
             type = "response",
             se.fit = T)

## Difference of probabilities 
p1 <- p$fit[2]
p0 <- p$fit[1]
p1-p0

p1-p0 - 1.96*sqrt(p1*(1-p1)/245 + p0*(1-p0)/240)
p1-p0 + 1.96*sqrt(p1*(1-p1)/245 + p0*(1-p0)/240)

## same as
p1-p0 - 1.96*sqrt(sum((p$se.fit)^2))
p1-p0 + 1.96*sqrt(sum((p$se.fit)^2))