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

## we wish to use "parity = 1+" as reference level
fever$parity <- relevel(fever$parity, ref = "1")

## weights at 0 is really missing values
fever$weight <- ifelse(fever$weight == 0, NA, fever$weight)
fever$ga <- cut(fever$gwbirth, c(0,31:40, Inf))

## cutting weight at the 5% percentile for each ga-group
pct5 <- tapply(fever$weight, fever$ga, quantile, probs = .05, na.rm=T)

## Small gestational age?
fever$sga <- (fever$weight < pct5[fever$ga])

## Each observation without small gestational age is selected with probability 0.15

U <- runif(length(fever$sga[fever$sga == 0]))
casecon <- rbind(subset(fever, sga == 1),
                 subset(fever, sga == 0)[(U < 0.15),])

## NB new control sample (and a bit different number of sga cases)
## gives different results to those of the book

## logistic regression on Case/Control dataset
summary(model.casecon <- glm(sga ~ parity + parity:smoke, family = binomial, data = casecon))

## logistic regression on the full Fever dataset
summary(model.fever <- glm(sga ~ parity + parity:smoke, family = binomial, data = fever))