BJnoint.R
``````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 30.  minimum iteration
#                           is internally set to 3.
#                error = when the consecutive iteration changes less then
#                        error, the iteration will stop. default is .00001
#                beta0 = initial estimator.
# 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)
}
``````