The hpstat \textit{intcox} command implements Cox proportional hazards regression for interval-censored observations, using an iterative convex minorant algorithm. The general algorithm is proposed by Huang \& Wellner [1]. When computing the baseline hazard, we apply a damped iterative convex minorant algorithm described by Aragón \& Eberly [2] and Pan [3]. This documentation discusses technical details of the implementation.
Let the baseline cumulative hazard be a step function $Λ(t)$ with steps at times $\{t_1, t_2, …, t_m\}$. We seek to optimise for the vector $\blambda=(Λ(t_1), Λ(t_2), …, Λ(t_m))^\mathrm{T}$. The cumulative hazard $Λ(t; \symbf{Z}_i)$ given covariates $\symbf{Z}_i$ is $\exp(\symbf{Z}_i^\mathrm{T}\bbeta) Λ(t)$. The cumulative distribution function for failure times is $F(t; \symbf{Z}_i)=1-\exp\left(-\exp(\symbf{Z}_i^\mathrm{T}\bbeta) Λ(t)\right)$.
Where $A_{ij}=\partial\left( S_{Li}- S_{Ri}\right)/\partial Λ(t_j)= S_{Ri}\exp(\symbf{Z}_i^\mathrm{T}\bbeta)\mathrm{I}(t_j = R_i)- S_{Li}\exp(\symbf{Z}_i^\mathrm{T}\bbeta)\mathrm{I}(t_j = L_i)$. The sum of all $\nablasub{\blambda}\mathcal{L}_i$ yields $\nablasub{\blambda}\mathcal{L}$.
Note that $\partial A_{ij}/\partial Λ(t_j)=-A_{ij}\exp(\symbf{Z}_i^\mathrm{T}\bbeta)$, so applying the quotient rule and simplifying, the Hessian $\nabla^2_{\blambda}\mathcal{L}_i$ has diagonal $(j, j)$-th elements:
Let $\symbf{G}$ be a diagonal matrix of the diagonal elements of $-\nabla^2_{\blambda}\mathcal{L}$. As discussed by Pan [3], we update $\blambda$ by iteratively applying:
Where $\mathrm{Proj}$, as defined by Pan [3], essentially represents monotonic (isotonic) regression. We port an efficient pool adjacent violators algorithm (PAVA) implementation from scikit-learn.
Note $\partial B_{Li}/\partial\bbeta=\exp(\symbf{Z}_i^\mathrm{T}\bbeta) Λ(L_i)(S_{Li}- B_{Li})\symbf{Z}_i$, and $\partial B_{Ri}/\partial\bbeta$ is analogous. Applying the quotient rule and simplifying, the Hessian $\nabla^2_{\blambda}\mathcal{L}_i$ is:
The covariance matrix for $\bbeta$ is computed as suggested by Zeng, Gao \& Lin [4] as the inverse of the empirical Fisher information, which is in turn estimated as the outer product of the gradient of the profile log-likelihood.
\item Huang J, Wellner JA. Interval censored survival data: a review of recent progress. In: Lin DY, Fleming TR, editors. \textit{Proceedings of the First Seattle Symposium in Biostatistics: survival analysis}; 1995 Nov 20–21; University of Washington, Seattle. New York: Springer-Verlag; 1997. p. 123–69. \href{https://doi.org/10.1007/978-1-4684-6316-3_8}{doi: 10.1007\slash 978-1-4684-6316-3\_8}
\item Aragón J, Eberly D. On convergence of convex minorant algorithms for distribution estimation with interval-censored data. \textit{Journal of Computational and Graphical Statistics}. 1992;1(2):129–40. \href{https://doi.org/10.2307/1390837}{doi: 10.2307\slash 1390837}
\item Pan W. Extending the iterative convex minorant algorithm to the Cox model for interval-censored data. \textit{Journal of Computational and Graphical Statistics}. 1999;8(1):109–20. \href{https://doi.org/10.1080/10618600.1999.10474804}{doi: 10.1080\slash 10618600.\allowbreak{}1999.10474804}
\item Zeng D, Gao F, Lin DY. Maximum likelihood estimation for semiparametric regression models with multivariate interval-censored data. \textit{Biometrika}. 2017;104(3):505–25. \href{https://doi.org/10.1093/biomet/asx029}{doi: 10.\allowbreak{}1093\slash biomet\slash asx029}
