R : Copyright 2003, The R Development Core Team
Version 1.6.2  (2003-01-10)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type `license()' or `licence()' for distribution details.

R is a collaborative project with many contributors.
Type `contributors()' for more information.

Type `demo()' for some demos, `help()' for on-line help, or
`help.start()' for a HTML browser interface to help.
Type `q()' to quit R.

> invisible(options(echo = TRUE))
> X <- read.table(url("http://www.stat.umn.edu/geyer/8701/parm/kyphosis.txt"),
+     header = TRUE)
> names(X)
[1] "present" "Age"     "Start"   "Number" 
> 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)

Call:
glm(formula = present ~ Number + Start, family = binomial)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9035  -0.5938  -0.3886  -0.2490   2.2141  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept) -1.02890    1.22533  -0.840  0.40108   
Number       0.35745    0.20016   1.786  0.07412 . 
Start       -0.18495    0.06412  -2.884  0.00392 **
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 83.234  on 80  degrees of freedom
Residual deviance: 64.536  on 78  degrees of freedom
AIC: 70.536

Number of Fisher Scoring iterations: 4

> out.big <- glm(present ~ Age + Number + Start, family=binomial)
> summary(out.big)

Call:
glm(formula = present ~ Age + Number + Start, family = binomial)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.3123  -0.5484  -0.3632  -0.1659   2.1613  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)   
(Intercept) -2.036708   1.443836  -1.411  0.15836   
Age          0.010929   0.006416   1.703  0.08849 . 
Number       0.410564   0.223772   1.835  0.06654 . 
Start       -0.206501   0.067498  -3.059  0.00222 **
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 83.234  on 80  degrees of freedom
Residual deviance: 61.380  on 77  degrees of freedom
AIC: 69.38

Number of Fisher Scoring iterations: 4

> 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"]
+     }
+ }
Error in (if (is.empty.model(mt)) glm.fit.null else glm.fit)(x = X, y = Y,  : 
	(converted from warning) Algorithm did not converge
Error in (if (is.empty.model(mt)) glm.fit.null else glm.fit)(x = X, y = Y,  : 
	(converted from warning) fitted probabilities numerically 0 or 1 occurred
> 
> 2 * pnorm(- abs(z.hat))
[1] 0.08848842
> p.val <- (sum(abs(z.star) >= z.hat, na.rm = TRUE) +
+     sum(is.na(z.star)) + 1) / (nboot + 1)
> p.val
[1] 0.08708
> p.val + c(-1, 1) * qnorm(0.975) * sqrt(p.val * (1 - p.val) / nboot)
[1] 0.08533247 0.08882753
> 
> postscript(file = "parm.ps")
> qqnorm(z.star)
> abline(0, 1)
> dev.off()
null device 
          1 
> system("pstopnm -nocrop -yborder 0 -xborder 0 -portrait parm.ps")
pstopnm: Writing ppmraw file
> system("pnmcrop parm001.ppm | pnmflip -rotate270 | ppmtogif > parm.gif")
ppmtogif: computing colormap...
ppmtogif: 2 colors found
> 
> cat("Computer name is", system("uname -n", intern = TRUE), "\n")
Computer name is prior 
> cat("Calculation took", proc.time()[1], "seconds\n")
Calculation took 2608.68 seconds
> proc.time()
[1] 2608.68   27.13 2650.21    2.94    0.13
>