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

## model with splines
model <-  lm(log10(bw) ~  log10(ad) + log10(bpd) + I(log10(bpd)^2) +
             I((log10(bpd)-log10(85))^2 * ifelse(bpd >= 85, 1, 0)) +
             I((log10(bpd)-log10(90))^2 * ifelse(bpd >= 90, 1, 0)) +
             I((log10(bpd)-log10(95))^2 * ifelse(bpd >= 95, 1, 0)),
             data = secher)

## predicting on n values of bpd for 3 values of ad.
n <- 1000 
r2 <- range(secher$bpd, na.rm=TRUE)
newbpd <- seq(min(r2), max(r2), length.out = n)

linpred90 <- predict(model, data.frame(ad = rep(90, n), bpd = newbpd))
linpred100 <- predict(model, data.frame(ad = rep(100, n), bpd = newbpd))
linpred110 <- predict(model, data.frame(ad = rep(110, n), bpd = newbpd))

## the 3 curves
plot(linpred90~log10(newbpd), type = "l",
     ylim = c(3,3.6),
     xlab = expression(paste(log[10],"(biparietal diameter)", sep="")),
     ylab = "Linear predictor"
     )
lines(linpred100~log10(newbpd), lty = "32")
lines(linpred110~log10(newbpd), lty = "11")