# Rank AFT: Setup
require(survival)
test <- function(r, dd, xx) {
  n <- length(r)
  ind <- order(r)
  x <- xx[ind]
  d <- dd[ind]
  w <- v <- 0
  for (i in 1:length(x)) {
    w <- w + d[i] * (x[i] - mean(x[i:n]))
    v <- v + d[i] * (x[i] - mean(x[i:n]))^2
  }
  w/sqrt(v)
}

# Pike Rat
Data <- read.delim("http://myweb.uiowa.edu/pbreheny/data/Pike1966.txt")
ind <- order(Data$Time)
Data <- Data[ind,]
y <- log(Data$Time)
d <- 1*(Data$Death)
x <- as.numeric(Data$Group)-1
bb <- seq(-0.2,0.6,len=199)
z <- sapply(bb, function(b) {test(y-x*b, d, x)})
plot(bb, z^2, type='s', bty='n', lwd=3, las=1,
     xlab=expression(beta), ylab=expression(chi^2))
abline(h=qchisq(.95,1), col='gray')

# Comparison table
B <- matrix(NA, 4, 3, dimnames=list(c('Rank-based', 'Weibull', 'Exponential', 'Lognormal'), c('Est', 'Lower', 'Upper')))
B[1,] <- c(bb[which.min(abs(z))],
           bb[which.min(abs(z+1.96))],
           bb[which.min(abs(z-1.96))])
fit <- survreg(Surv(Data$Time, d)~x)
B[2,] <- c(coef(fit)[2], confint(fit, 2))
fit <- survreg(Surv(Data$Time, d)~x, dist='exponential')
B[3,] <- c(coef(fit)[2], confint(fit, 2))
fit <- survreg(Surv(Data$Time, d)~x, dist='lognormal')
B[4,] <- c(coef(fit)[2], confint(fit, 2))
B

# Simulated data
set.seed(3)
n <- 1000
x <- rbinom(n, 1, 0.5)
lam <- exp(-x)
yy <- rweibull(n, 1.5, 1/lam)
cc <- rexp(n)
y <- log(pmin(yy, cc))
d <- yy <= cc
bb <- seq(0,2,len=199)
z <- sapply(bb, function(b) {test(y-x*b, d, x)})
plot(bb, z^2, type='s', bty='n', lwd=3, las=1,
     xlab=expression(beta), ylab=expression(chi^2))
abline(h=qchisq(.95,1), col='gray')

# Comparison table
B <- matrix(NA, 4, 3, dimnames=list(c('Rank-based', 'Weibull', 'Exponential', 'Lognormal'), c('Est', 'Lower', 'Upper')))
B[1,] <- c(bb[which.min(abs(z))],
           bb[which.min(abs(z+1.96))],
           bb[which.min(abs(z-1.96))])
fit <- survreg(Surv(exp(y), d)~x)
B[2,] <- c(coef(fit)[2], confint(fit, 2))
fit <- survreg(Surv(exp(y), d)~x, dist='exponential')
B[3,] <- c(coef(fit)[2], confint(fit, 2))
fit <- survreg(Surv(exp(y), d)~x, dist='lognormal')
B[4,] <- c(coef(fit)[2], confint(fit, 2))
B