Cox model for interval censored data
We give two methods for computing the baseline cumulative hazard function and the regression parameter by maximum likelihood for the Cox proportional hazards model: a profile likelihood method and a method proposed by Wei Pan (Journal of Computational and Graphical Statistics Vol. 8, No. 1 (1999), pp. 109-120). Both methods use the iterative convex minorant algorithm, as proposed in Groeneboom (Technical Report 378, Department of Statistics, Stanford University, 1991), and also discussed in Groeneboom and Wellner, Information bounds and nonparametric maximum likelihood estimation (Volume 19 of DMV Seminar. Birkhäuser Verlag, Basel, 1992). The profile likelihood method uses the derivative-free search algorithm of Hooke and Jeeves (https://en.wikipedia.org/wiki/Pattern_search_(optimization)) for computing the estimate of the finite dimensional regression parameter. There has been an R package intcox, implementing the method of Wei Pan, but this package is no longer supported on CRAN.
We demonstrate the method for the following model. The baseline cumulative hazard function is the function H_0(z)=z^3, so we use the Weibull(3,1) distribution function for the baseline hazard. The finite dimensional parameter beta is the vector (-1,1,1/2) and the covariates are of the form (U,V,W), where U and V are Uniform(0,1) and W is a Bernoulli random variable with parameter p=1/2. The R file Cox_model.R generates 100 samples of size 500 from this Cox model, but the only information we have about the generated values comes via interval censoring. If X is a generated value of the Cox model, our only information are the values of the indicators 1_{X<=T}, 1_{T<X<=U} and 1_{X>U}, where in our example T and U are the order statistics of a uniform sample of size 2 on the square [0,2]^2.
In our example the profile likelihood method is behaving somewhat better than the method of Wei Pan. In general one can expect the profile likelihood method to be more robust, since there are no problems with diverging Hessians, which quasi-Newton methods like Pan's method are prone to. Convergence properties of the Hooke-Jeeves method are proved in several interesting papers of Virginia Torczon and co-authors. One can run the example via Cox_model.R. This produces Boxplots of the rescaled Euclidean norm of the difference between the estimate of beta and beta itself, the time needed to compute the estimates and also the estimate of the baseline cumulative hazard and the baseline distribution function of the last sample of the simulation. The R file Cox_model.R uses the R package Rcpp, so one must have Rcpp installed to run the example (no other packages are needed). The author welcomes comments on the program, which is just the start of a more comprehensive package.