BJnoint <- function(x, y, delta, beta0 = NA, maxiter = 30, error = 0.00001)
{
# This is an R function to compute the Buckley-James estimator for
# censored regression WITHOUT intercept term. (but you may still force x
# to have a column of 1's). It calls another function iter(). This function
# use to call a C function. Written by mai zhou, mai@ms.uky.edu
# First (C) version Jan. 1994. Last C revision: May 10, 1995
# R version: 2004 Jan 20. R speed actually is OK, no need of C function.
#
# Input:
# x is a matrix of N rows (covariates).
# y is the observed (censored) responses --- a vector of length N.
# delta is a vector of length N. delta =1 means (y) is not censored.
# delta = 0 means y is right censored, i.e. the true
# response is larger than y.
#
# Optinal input: maxiter = number of maximum iteration, default is 20.
# minimum iteration is internally set to 3.
# error = when the consecutive iteration changes less then
# error, the iteration will stop. default is .00001
# Output:
# the estimate, beta, and an extra integer at the end: number of iterations.
#
# Bug: More careful for ties. Do something better when last obs. is censored.
x <- as.matrix(x) # to insure x is a matrix, not a vector
newtemp <- matrix(NA, ncol = ncol(x), nrow = 3)
newtemp[1,] <- beta0
if(any(is.na(beta0)))
newtemp[1, ] <- lm(y ~ 0 + x)$coef # get initial est.
for(i in 2:3) {
newtemp[i, ] <- iter(x, y, delta, newtemp[i - 1, ])
}
num <- 2 # do at least 2 iterations
while(num <= maxiter && error <= sum(abs(newtemp[2, ] - newtemp[3, ])))
{
newtemp[2, ] <- newtemp[3, ]
newtemp[3, ] <- iter(x, y, delta, newtemp[2, ])
num <- num + 1
}
if(num > maxiter) # always take average? or only when not converge?
{newtemp[3, ] <- (newtemp[2, ] + newtemp[3, ])/2}
list(beta=newtemp[3, ], iteration=num)
}