Efficient estimation of partially linear additive Cox models and variance estimation under shape restrictions

Shape restrictions (such as monotonicity and convexity) arise naturally in numerous practical scenarios. For instance, the growth curves of animals and plants in ecology and the dose-response in medicine must inherently exhibit non-decreasing characteristics (Chang et al., 2007; Wang and Ghosh, 2012). In the realm of economics, utility and production functions are often concave in income and prices (Matzkin, 1991; Varian, 1984), cost functions are monotone increasing, concave in input prices, and may exhibit non-increasing or non-decreasing returns to scale (Horowitz and Lee, 2017). In genetic epidemiology studies, the cumulative risk of a disease for individuals possess monotonicity (Qin et al., 2014). While in reliability analysis, the bathtub curve describing the failure rate typically displays convexity.

Incorporating shape restrictions into statistical analysis, apart from its exceptional interpretability and ability to enforce domain-specific constraints, often results in an estimation procedure that is devoid of tuning parameters, enhancing its efficiency and robustness. Therefore shape-restricted techniques has become an increasing popular tool for statistical inference or learning in various settings over the past decades. A comprehensive review on shape-restricted nonparametric inferences can be found in Groeneboom and Jongbloed (2014) and references therein. Recently, Chen and Samworth (2016) developed an algorithm for the estimation of the generalized additive model in which each of the additive components is linear or subject to a shape restriction. Balabdaoui et al. (2019) considered the estimation of the index parameter in a single-index model with a monotonically increasing link function. Deng and Zhang (2020) studied minimax and adaptation rates in general multiple isotonic regression. Feng et al. (2022) systematically investigate the theoretical properties of the least squared estimator of a S-shaped regression function.

This paper focus on the statistical inference for righter-censored survival data. Let $T$ denote the survival time and $(Z,X)\in\mathbb{R}^{p}\times\mathbb{R}^{d}$ denote a $(p+d)\times 1$ vector of covariates. We consider the partially linear Cox model of (Sasieni, 1992, PLCM) for modelling the conditional hazard function, i.e.

where $\lambda(\cdot)$ is the unspecified baseline hazard function, $\beta\in\mathbb{R}^{d}$ is unspecified and $g(\cdot):\mathbb{R}^{p}\mapsto\mathbb{R}$ is an unknown function. This model reduces to the renowned Cox proportional hazards model (Cox, 1972, 1975) when the covariate $Z$ disappears, and it becomes the nonparametric Cox model (Sleeper and Harrington, 1990; O’Sullivan, 1993) in the absence of $X$.

Many nonparametric techniques have been developed for the estimation of the PLCM, in particular for the linear covariate effect. Examples include profile partial likelihood together with a kernel technique (Heller, 2001), maximum likelihood estimation with a deep neural network (Zhong et al., 2022), and a kernel machine representation method (Rong et al., 2024), etc. However, these methods are either hampered by the curse of dimensionality or lack interpretability for $g(Z)$, or suffer from tuning parameters, whose selection is not always straightforward. Alternatively, Huang (1999) proposed to model $g(Z)$ by a generalized additive model (Hastie and Tibshirani, 1986, 1990), which effectively avoids the curse of dimensionality and enforces an additive effect for the covariate $Z$. Specifically,

where for $1\leq j\leq p$, $Z_{(j)}$ is the $j$-th component of $Z$ and $g_{j}$ is an unknown function. Huang (1999) proposed the use of polynomial splines to fit the unknown additive components. This method entails a number of tuning parameters, also yields convergence rates that lack conciseness and elegance. Furthermore, the spline method does not provide good interpretability for the additive covariate $Z_{(j)}$.

Our paper is motivated by the work of Deng et al. (2023) which studied a shape-restricted and additive PLCM. Specifically, under models (1) and 2, they assume that each $g_{j}$ is monotonic increasing/decreasing or convex/concave. An active-set optimization algorithm was provided to calculate the shape-restricted maximum likelihood estimator. The shape-restriction strategy facilitates the utilization of prior knowledge regarding the effect of the log conditional hazard function on each covariate $Z_{(j)}$ and leads to a tuning-parameter-free estimation procedure. However, they proved only a consistency result, and did not provide any asymptotic normality results. Qin et al. (2021) studied a PLCM with a single additive component subject to shape restrictions, but they did not establish any $\sqrt{n}$-consistency result. In addition, in general shape-restriction inferences, even if asymptotic normality results can be established, it is generally challenging to construct reasonable estimators for the asymptotic variances with theoretical guarantees (Groeneboom and Hendrickx, 2017).

This paper makes two main contributions to the literature of additive and shape-restricted PLCMs for survival data. The first contribution is to provide powerful statistical guarantees for the shape-restricted maximum partial likelihood estimator (SMPLE) and the induced Breslow-type estimator for the baseline cumulative hazard function under the model assumption of Deng et al. (2023). This includes a thorough convergence rate analysis for the estimators of the infinitely dimensional parameters, as well as establishing asymptotic normality and semiparametric efficiency for the estimator of the linear covariate effect. Our second contribution is to offer an easy-to-use estimator for the asymptotic variance of the linear covariate effect estimator. We show that this variance estimation method always provide consistent estimators once the corresponding asymptotic normality result holds. This method is very flexible and is applicable for general purpose especially in shape-restricted inferences, where theoretical guarantee of a bootstrap variance estimator is generally rather challenging (Groeneboom and Hendrickx, 2017). Our simulation results and an analysis of the Rotterdam Breast Cancer dataset demonstrate that the SMPLE has comparable performance with the maximum likelihood estimator under the Cox model when the Cox model is correct, and outperforms the latter and Huang (1999)’s method when the Cox model is violated and the hazard is nonsmooth. Meanwhile, the proposed variance estimation method usually leads to reliable interval estimates for the SMPLE and its competitors.

The rest of this paper is organized as follows. Section 2 introduces notations, data, and the shape-restricted maximum partial likelihood estimators (SMPLE). Section 3 investigates the convergence rates of the SMPLEs for all the unknown parameters, including $\beta$, the unknown additive components, and the baseline cumulative hazard function. Section 4 establishes the asymptotic normality and semiparametric efficiency of the SMPLE for $\beta$. A novel estimation method is also provided to estimate the asymptotic variance of the SMPLE of $\beta$. A simulation study and real data analysis are presented in Section 5 and 6, respectively. Section 7 contains concluding remarks. For clarity, all technical proofs are postponed to the supplementary material.

The consistency property of the SMPLEs $\hat{\beta}$, $\hat{g}$, and $\hat{g}_{j}$ was established by Deng et al. (2023). In this section, we establish their convergence rates. We make the following assumptions.

Let $\|\cdot\|$ denote the usual Euclidean norm and $\|f(\cdot)\|_{\infty}$ the supreme norm of a real-valued function $f$. For any constant $M>0$, define

Assumptions 1–2 are standard in the theoretical analysis of traditional Cox model and its variants (Huang, 1999; Zhong et al., 2022). Assumption 3 ensures that the probability of being uncensored is positive regardless of the covariate values, and it is used to establish the convergence rate results in Theorem 1. Assumption 4 is used in the calculation of the semiparametric efficiency lower bound (Huang, 1999). In Assumption 5, the upper bound requirement is used to calculate some entropy results needed in the proof of Proposition 1, and the lower bound requirement guarantees that the approximation errors of piecewise linear approximations of the convex/concave additive components to themselves are small enough in the proof of Theorem 3.

The Proposition below establishes the consistency of $R_{\hat{\eta}}(\cdot)$, as an estimator of $R_{\eta_{0}}(\cdot)$, which roughly implies the consistency of $\hat{\eta}=(\hat{\beta},\hat{g})$. We would have proved the consistencies of $\hat{\beta}$ and each $\hat{g}_{j}$ separately. However, the latter results are not needed in the proofs of the subsequent convergence rate results given the consistency of $R_{\hat{\eta}}(\cdot)$.

Define $d^{2}(\eta,\eta_{0})=\mathbb{E}_{U}\left\{(R_{\eta}(U)-R_{\eta_{0}}(U))^{2}\right\},$ where $\mathbb{E}_{U}$ denotes the expectation with respect to $U$. Let $\|\cdot\|_{L_{2}}$ denote the $L_{2}(P)$ norm and $\rho(\boldsymbol{k}_{0})=0.5+0.5\cdot\mathbf{1}(\cup_{i=1}^{p}({\rm sha}(g_{0,i})\in\{1,2\}))$. One of our main results is to establish the convergence rate of the SMPLE $\hat{\eta}$.

According to Theorem 1, the rates of convergence of all $\hat{g}_{j}$ are $O_{p}(n^{-2/5})$ if none of the additive components of $g$ is monotonic. Conversely, if one additive component of $g$ is monotonic, then their convergence rates all slow down to $O_{p}(n^{-1/3})$. An explanation for this finding is that the complexity of the class of bounded and monotonic functions is much larger than that of the class of bounded and convex (or concave) functions. These convergence rate results are free from the covariate dimensionality and exhibit a much more elegant form than those in Huang (1999) and Zhong et al. (2022). Theorem 1 also establishes the convergence rate of $\hat{\beta}$, although it is sub-optimal.

With Theorem 1, we are able to establish the uniformly rate of convergence for the SMPLE $\hat{\Lambda}(y;\hat{\eta})$ in (4) of the baseline cumulative hazard function $\Lambda_{0}(y)$. It turns out that $\hat{\Lambda}(y;\hat{\eta})$ has the same convergence rate as $\hat{\eta}$, although their convergence rates are quantified by different distances.

Based on the convergence rate results in the previous section, in this section, we further show that our SMPLE in (3) for the linear covariate effect $\hat{\beta}$ is asymptotically normal and semiparametric efficient, in the sense that its asymptotic variance achieves the semiparametric efficiency lower bound (Bickel et al., 1993) or the information bound of estimating $\beta_{0}$ under models (1) and (2).

We begin with presenting the information bound for $\beta_{0}$. Recall that $U=(X^{\top},Z^{\top})^{\top}$, and define

which is a counting process martingale associated with the Cox model. The log-likelihood of model (1) based on one observation $(X,Z,Y,\Delta)$ is (up to constant)

Conisder a parametric smooth sub-model $\{\lambda_{(\nu)}:\nu\in\mathbb{R}\}$ and $\{g_{j,(\nu)}:\nu\in\mathbb{R}\}$, $1\leq j\leq p$, with $\lambda_{(0)}=\lambda_{0}$ and $g_{j,(0)}=g_{0,j}$. Define $L_{2}(P_{Y})$ to be the set of $a(\cdot)$ satisfying ${\mathbb{E}}\{\Delta a^{2}(Y)\}<\infty$ and $a(y)=\partial\log\lambda_{(\nu)}(y)/\partial\nu|_{\nu=0},$ for some submodel. Similarly, for $1\leq j\leq p$, define $L_{2}^{0}(P_{Z_{(j)}})$ to be the set of $h_{j}$ satisfying ${\mathbb{E}}\{\Delta h_{j}(Z_{(j)})\}=0$, ${\mathbb{E}}\{\Delta h_{j}^{2}(Z_{(j)})\}<\infty$ and $h_{j}(z_{(j)})=\partial g_{j,(\nu)}(z_{(j)})/\partial\nu|_{\nu=0}$ for some submodel. The following lemma, which is Theorem 3.1 of Huang (1999), gives the information bound of $\beta_{0}$.

Additional assumptions are needed to obtain the asymptotic normality and efficiency of $\hat{\beta}$. Denote $\boldsymbol{h}_{j}^{\star}=(\boldsymbol{h}_{j,1}^{\star},\cdots,\boldsymbol{h}_{j,d}^{\star})^{\top}$ for $1\leq j\leq p$.

Assumption 6 is used to control the fluctuation of the score function corresponding to the monotone components in the direction of the projection defined in Lemma 1. (Huang, 2002; Cheng, 2009) adopted a similar assumption. Assumptions 7–8, which are analogues of Assumption B1 of Kuchibhotla et al. (2023), are used for approximations of the convex (concave) components.

By Theorem 3, $\hat{\beta}$ has an asymptotically normal distribution with asymptotic variance $I^{-1}(\beta_{0})$. When making inference about $\beta_{0}$ based on this theorem, we need to construct a consistent estimator for $I^{-1}(\beta_{0})$. However, $I^{-1}(\beta_{0})$ or equivalently $I(\beta_{0})$ has a rather complicated form, making its plug-in estimator not easy to use. To crack this nut, we propose a novel data-splitting estimation method to estimate $I^{-1}(\beta_{0})$.

In this section, we conduct simulations to assess the finite-sample performance of the proposed SMPLE $\hat{\beta}$ and the proposed confidence region (9) for the linear covariate effect $\beta$. To generate data, we take $X$ and $Z$ to be two scalar random variables, which are iid from the standard normal distribution, and take the conditional distribution of $T$ given $(X,Z)$ to be an exponential distribution with mean $1/\exp(\beta_{0}X+g_{0}(Z))$. Therefore the conditional hazard function of $T$ given $(X,Z)$ is $\exp(\beta_{0}X+g_{0}(Z))$. We set the censoring time $C$ to follow a uniform distribution on $(0,c)$. We consider three scenarios for $g_{0}$: (I) $g_{0}(z)=-2z$, (II) $g_{0}(z)=-|z|^{3}/2$ and (III) $g_{0}(z)=2|z|$. We set $\beta_{0}=-2$ and consider two choices for $c$: $5$ and $10$, and two sample sizes: $600$ and $800$. A larger $c$ results in a smaller censoring proportion.

When implementing our SMPLE, we set $\boldsymbol{k}_{0}$ to be $3,4,3$ in Scenarios I–III, respectively. For comparison, we also consider the traditional Cox regression estimator (TCR) of $\beta$ and the partial likelihood estimator of $\beta$ with the $r$-order polynomial splines under the partially linear additive model (Huang, 1999, SPLA-$r$), where $r$ may be $2,3$ or $4$. We generate 1000 random samples to evaluate the performance of the above five estimation methods.

We apply the proposed SMPLE method to analyze the Rotterdam Breast Cancer (RBC) dataset, which is publicly available from the R package survival. This dataset was used in Royston and Altman (2013) to perform external validation of a Cox prognostic model. The RBC dataset comprises 2982 primary breast cancers patients whose records were included in the Rotterdam tumor bank. We focus on studying potential factors that affect the survival time ($T$) of breast cancers patients, defined as the days from primary surgery to the death. We consider four covariates: progesterone receptors (unit: fmol/l, $X$), age at surgery ($Z_{1}$), number of positive lymph nodes ($Z_{2}$), estrogen receptors (unit: fmol/l, $Z_{3}$). Let $Z=(Z_{1},Z_{2},Z_{3})^{\mathrm{\scriptscriptstyle\top}}$. All covariates have been scaled to between $0$ and $1$. The impact of the content of progesterone receptors on the treatment of breast cancer is not clear. We consider modelling the conditional hazard of $T$ given $(X=x,Z=z)$ by the following partially linear Cox additive model

where $\lambda(\cdot)$, $f_{1}(\cdot)$ $f_{2}(\cdot)$ and $f_{3}(\cdot)$ are all left un-specified. Since elderly patients are typically predisposed to a higher risk of cancer deterioration or recurrence, we postulate that $f_{1}(\cdot)$ is monotonically increasing. Similarly, $f_{2}(\cdot)$ and $f_{3}(\cdot)$ are also assumed to be monotonically increasing functions. Inspired by our extensive simulation results and discussion in Remark 1, we further impose a convexity or concavity restriction on $f_{1}(\cdot)$–$f_{3}(\cdot)$.

Because menopause generally affects the therapeutic effect and the mortality of breast cancer, we categorize the patients into two groups based on whether they have reached menopause or not. Group I (premenopausal) comprises $1312$ samples, exhibiting a censor rate of $64.33\%$, while Group II (postmenopausal) consists of $1670$ samples with a censor rate of $51.86\%$. We apply the point and interval estimation methods considered in the simulation section to estimate the linear covariate effect under model (11). We set the confidence level to 95% when constructing confidence intervals. In the proposed variance estimation method, we take $\tilde{\alpha}=0.3$ and construct variance estimates by repeating the splitting and estimating procedure twenty times and taking average. Table 3 summarizes the analysis results.

The point estimates of $\beta$ are about $-5$ or less in Group I, which are much less than those in Group II, namely $-0.7$ or larger. This implies that the effect of progesterone receptors ($X$) on the hazard rate of survival time among breast cancer patients has a significant increase in Group II than in Group I, which potentially stems from the substantial reduction in progesterone levels in women’s bodies post-menopause. For patients who have not reached menopause, the $95\%$ confidence intervals derived from all five methods exclude $0$, suggesting a significant effect of the progesterone receptors on the survival time of breast cancer patients. Among these techniques, our proposed SMPLE method produces the shortest confidence interval. Conversely, for postmenopausal patients, only our SMPLE-based confidence interval excludes $0$, which together with the point estimate $-0.62$ suggests that the progesterone receptors still has a significant effect in increasing the survival time of breast cancer patients. Once again, our confidence interval exhibits the shortest confidence interval. According to our reasonable shape restriction assumptions and our simulation experience especially those in Scenario III, we believe that the analysis results based on our SMPLE are more reliable, which also demonstrate its priority over the competing methods.

Figure 3 depicts the estimates of $f_{1}(\cdot)$–$f_{3}(\cdot)$ in model (3) under our shape constraints. All the estimated curves are continuous and piecewise linear functions in both groups. We see that both number of positive lymph nodes ($Z_{2}$) and estrogen receptors ($Z_{3}$) have very similar effects in the two groups on the survival time $T$. Exceptionally, the covariate age at the time of surgery ($Z_{1}$) has a significantly more pronounced effect on $T$ in Group II than in Group I. This discrepancy may stem from the fact that premenopausal women tend to be younger, and thus, age has not yet exerted a substantial influence on their physical condition.