X <- read.table(url("http://www.stat.umn.edu/geyer/8701/parm/kyphosis.txt"),
    header = TRUE)
names(X)
attach(X)

.Random.seed <- c(1, 195929957, 784973954)

out.little <- glm(present ~ Number + Start, family=binomial)
p.hat <- predict(out.little, type = "response")
summary(out.little)
out.big <- glm(present ~ Age + Number + Start, family=binomial)
summary(out.big)
coef <- summary(out.big)$coefficients
z.hat <- coef["Age", "z value"]

nboot <- 1e5 - 1
z.star <- rep(NA, nboot)
n <- length(present)

options(warn = 2)
for (i in 1:nboot) {
    present.star <- as.numeric(runif(n) < p.hat)
    out.star <- try(glm(present.star ~ Age + Number + Start, family=binomial))
    if (! inherits(out.star, "try-error")) {
        coef.star <- summary(out.star)$coefficients
        z.star[i] <- coef.star["Age", "z value"]
    }
}

2 * pnorm(- abs(z.hat))
p.val <- (sum(abs(z.star) >= z.hat, na.rm = TRUE) +
    sum(is.na(z.star)) + 1) / (nboot + 1)
p.val
p.val + c(-1, 1) * qnorm(0.975) * sqrt(p.val * (1 - p.val) / nboot)

postscript(file = "parm.ps")
qqnorm(z.star)
abline(0, 1)
dev.off()
system("pstopnm -nocrop -yborder 0 -xborder 0 -portrait parm.ps")
system("pnmcrop parm001.ppm | pnmflip -rotate270 | ppmtogif > parm.gif")

cat("Computer name is", system("uname -n", intern = TRUE), "\n")
cat("Calculation took", proc.time()[1], "seconds\n")