A penalized likelihood approach for efficiently estimating a partially linear additive transformation model with current status data

TL;DR

This paper introduces a penalized likelihood method for efficiently estimating a flexible partially linear additive transformation model tailored for current status data, addressing challenges in medical and epidemiological research.

Contribution

It proposes a novel penalized likelihood approach with constrained B-splines for modeling monotone transformations and nonlinear effects, along with an efficient hybrid algorithm for estimation.

Findings

Estimators of regression coefficients are root-n consistent and asymptotically normal.

Nonparametric estimators achieve optimal convergence rates.

Method performs well in simulations and real data analysis.

Abstract

Current status data are commonly encountered in medical and epidemiological studies in which the failure time for study units is the outcome variable of interest. Data of this form are characterized by the fact that the failure time is not directly observed but rather is known relative to an observation time; i.e., the failure times are either left- or right-censored. Due to its structure, the analysis of such data can be challenging. To circumvent these challenges and to provide for a flexible modeling construct which can be used to analyze current status data, herein, a partially linear additive transformation model is proposed. In the formulation of this model, constrained $B$-splines are employed to model the monotone transformation function and nonlinear covariate effects. To provide for more efficient estimates, a penalization technique is used to regularize the estimation of all unknown functions. An easy to implement hybrid algorithm is developed for model fitting and a simple estimator of the large-sample variance-covariance matrix is proposed. It is shown theoretically that the proposed estimators of the finite-dimensional regression coefficients are root-$n$ consistent, asymptotically normal, and achieve the semi-parametric information bound while the estimators of the nonparametric components attain the optimal rate of convergence. The finite-sample performance of the proposed methodology is evaluated through extensive numerical studies and is further demonstrated through the analysis of uterine leiomyomata data.

Figures

Tables

Table 1. Table 1: Summary of the results that were obtained using the spline estimator and the penalized spline estimator when α = 0 𝛼 0 \alpha=0 (PH model), α = 0.5 𝛼 0.5 \alpha=0.5 , and α = 1 𝛼 1 \alpha=1 (PO model). Presented results include the empirical bias (Bias), sample standard deviation (SD) and mean squared error (MSE) of the 500 point estimates, as well as the average of the 500 estimated standard errors (SE) and empirical coverage probabilities associated with 95% confidence intervals (CP)

		Spline approach					Penalized approach
$\alpha =0.0$		Bias	SD	MSE	SE	CP	Bias	SD	MSE	SE	CP
S1	${\beta }_1$	0.065	0.295	0.092	0.252	91.4	-0.012	0.234	0.055	0.225	93.6
	${\beta }_2$	-0.078	0.174	0.036	0.139	88.0	-0.014	0.129	0.017	0.122	94.2
S2	${\beta }_1$	0.070	0.288	0.088	0.252	91.6	0.022	0.250	0.063	0.232	94.8
	${\beta }_2$	0.071	0.162	0.031	0.137	90.0	0.028	0.140	0.020	0.125	93.4
S3	${\beta }_1$	-0.064	0.327	0.111	0.284	91.8	-0.030	0.292	0.086	0.259	93.4
	${\beta }_2$	-0.071	0.177	0.036	0.149	90.6	-0.034	0.158	0.026	0.136	92.0
$\alpha =0.5$		Bias	SD	MSE	SE	CP	Bias	SD	MSE	SE	CP
S1	${\beta }_1$	0.084	0.303	0.099	0.278	93.4	0.023	0.268	0.072	0.257	93.8
	${\beta }_2$	-0.060	0.166	0.031	0.147	91.8	-0.013	0.145	0.021	0.135	93.6
S2	${\beta }_1$	0.049	0.304	0.095	0.283	92.4	0.006	0.275	0.075	0.265	93.4
	${\beta }_2$	0.047	0.171	0.031	0.150	92.8	0.010	0.155	0.024	0.139	93.4
S3	${\beta }_1$	-0.073	0.356	0.132	0.324	93.0	-0.041	0.325	0.107	0.298	93.2
	${\beta }_2$	-0.061	0.177	0.035	0.169	94.4	-0.025	0.159	0.026	0.155	94.4
$\alpha =1.0$		Bias	SD	MSE	SE	CP	Bias	SD	MSE	SE	CP
S1	${\beta }_1$	0.036	0.326	0.107	0.306	93.4	-0.013	0.294	0.087	0.288	94.2
	${\beta }_2$	-0.050	0.158	0.027	0.159	95.4	-0.015	0.144	0.021	0.149	96.6
S2	${\beta }_1$	0.073	0.328	0.113	0.315	93.6	0.032	0.298	0.090	0.297	95.2
	${\beta }_2$	0.041	0.176	0.033	0.164	93.8	0.007	0.164	0.027	0.154	93.6
S3	${\beta }_1$	-0.035	0.393	0.156	0.363	94.2	-0.007	0.351	0.123	0.333	93.4
	${\beta }_2$	-0.064	0.194	0.042	0.188	93.8	-0.021	0.173	0.030	0.172	95.8

Table 2. Table 2: Uterine fibroid data: Presented results include the parameter estimates, estimated standard errors (SE) and 95% Wald type confidence intervals (95% C.I.). The results were obtained using the penalized spline estimator under both the PH and PO model

Model	Parameter	Estimate	SE	95% C.I.
PH	Race	1.479	0.165	(1.154, 1.803)
	Parity	-0.305	0.144	(-0.588, -0.022)
	Obesity	0.114	0.174	(-0.227, 0.455)
PO	Race	1.648	0.200	(1.256, 2.040)
	Parity	-0.341	0.163	(-0.661 -0.022)
	Obesity	0.155	0.200	(-0.237 0.547)

Equations

$$\eta (T_i)=-Z_i^{⊺}\beta -\sum _{j=1}^J{\phi }_j(W_{ij})+{\epsilon }_i,$$

$$F(t\mid X_i)=\text{textsl}{g}^{-1}\left\{Z_i^{⊺}\beta +\eta (t)+\sum _{j=1}^J{\phi }_j(W_{ij})\right\},$$

$$l_{\lambda }(\tau )=\sum _{i=1}^n[{\Delta }_i\log F(Y_i\mid X_i)+$$

$${\mathcal{J}}^2(\eta )=\int \{{\eta }^{(r)}(t){\}}^{2}dt,{\mathcal{J}}^2({\phi }_j)=\int \{{\phi }_j^{(r)}(W){\}}^{2}dW.$$

$${\lambda }_j^2=o_p(n^{-1/2})\text{mbox}and{\lambda }_j^{-1}=O_p(n^{r/(2r+1)}),\text{mbox}for1\le j\le J.$$

$$n^{1/2}(\hat{\beta }-{\beta }_0)\stackrel{d}{⟶}\mathcal{N}\{0,{\mathcal{I}}^{-1}({\beta }_0)\},\text{mbox}asn\to \infty .$$

$$\eta (\cdot )\approx {\eta }_n(\cdot \mid \gamma )=\sum _{k=1}^{p_0}{\gamma }_k{b}_k(\cdot )\text{ and }{\phi }_j(\cdot )\approx {\phi }_{nj}(\cdot \mid {\alpha }_j^{\ast })=\sum _{k=1}^{p_j}{\alpha }_{jk}^{\ast }{b}_{jk}(\cdot ),$$

$$F_n(Y_i\mid X_i)=\text{textsl}{g}^{-1}\left\{Z_i^{⊺}\beta +M_i^{⊺}\gamma +\sum _{j=1}^J{B}_{ji}^{⊺}{\alpha }_j\right\},$$

$$l_{n,\lambda }(\theta )=\sum _{i=1}^n[{\Delta }_i\log F_n(Y_i\mid X_i)$$

$$-\frac{{\lambda }_0^2}{2}\Vert D_0^{(r)}\gamma {\Vert }^2-\sum _{j=1}^J\frac{{\lambda }_j^2}{2}\Vert D_j^{(r)}{\alpha }_j{\Vert }^2,$$

$$\displaystyle{D}_{j}^{(1,r)}=\left[\begin{tabular}[]{rrrrrr}$-1$&$1$&$0$&$\ldots$&$0$&$0$\\ $0$&$-1$&$1$&$\ldots$&$0$&$0$\\ $\vdots$&$\vdots$&$\vdots$&&$\vdots$&$\vdots$\\ $0$&$0$&$0$&$\ldots$&$-1$&$1$\end{tabular}\right]\in\mathbb{R}^{(p_{j}-r-2)\times(p_{j}-r-1)},$$

$$\Vert D_0^{(1)}\gamma {\Vert }^2$$

$$\Vert D_0^{(2)}\gamma {\Vert }^2$$

$${\theta }^{\ast }={\theta }^{(m)}+{\mathcal{I}}_{n,\lambda }^{-1}({\theta }^{(m)})\nabla l_{n,\lambda }({\theta }^{(m)}),$$

$$\nabla l_{n,\lambda }(\theta )$$

$${\mathcal{I}}_{n,\lambda }(\theta )$$

$$\Delta (\theta )$$

$$\Omega (\theta )$$

$${\gamma }^{(m+1)}=\underset{\gamma \in \mathcal{G}}{\mathrm{argmin}}\left\{\sum _{k=1}^{p_0}{\sigma }_k^2({\gamma }_k-{\gamma }_k^{\ast }{)}^2\right\},$$

$${\lambda }_j^{\ast 2}=\frac{\text{tr}\{S^{-}{S}_j\}-\text{tr}\left\{{\mathcal{I}}_{n,\lambda }({\tilde{\theta }}_{\lambda }{)}^{-1}{S}_j\right\}}{{{\tilde{\theta }}_{\lambda }}^{⊺}{S}_j{\tilde{\theta }}_{\lambda }}{\lambda }_j^2,\text{ for }0\le j\le J,$$

$$\Sigma (\tilde{\theta })=({\mathcal{X}}^{⊺}\Omega (\tilde{\theta })\mathcal{X}+S{)}^{-1}{\mathcal{X}}^{⊺}\Omega (\tilde{\theta })\mathcal{X}({\mathcal{X}}^{⊺}\Omega (\tilde{\theta })\mathcal{X}+S{)}^{-1}.$$

$$\text{textsl}{g}_{\alpha }\{F(t\mid X)\}=Z_1{\beta }_1+Z_2{\beta }_2+\eta (t)+{\phi }_1(W_1)+{\phi }_2(W_2),$$

$$\text{textsl}{g}_{\alpha }(u)=\{\begin{array}{ll}\log \frac{(1-u{)}^{-\alpha }-1}{\alpha }& \text{if }\alpha >0,\\ \log \{-\log (1-u)\}& \text{if }\alpha =0.\end{array}$$

$$\text{textsl}g\{F(t\mid X)\}=Z_1{\beta }_1+Z_2{\beta }_2+Z_3{\beta }_3+\eta (t)+\phi (W),$$

Full text

Keywords: $B$-splines; Current status data; Isotonic regression; Partially linear additive transformation model; Penalized estimation.

1 Introduction

Current status data is characterized by the fact that the failure/event time cannot be directly measured, but rather is known relative to a single observation time. That is, the failure time is known to be smaller or larger than an observation time leading to left- or right-censored data, respectively. Data of this structure commonly arise as a part of medical and epidemiological studies, among many others, in which study participants are observed only once due to destructive/invasive testing, resource limitations, etc. For example, as part of the Right From The Start (RFTS) study, participants completed an ultrasound examination as a means to determine their fibroid status (Laughlin et al.,, 2009). In this study, fibroid onset time could not be directly observed but was known relative to the participant’s age at the ultrasound examination time.

For the regression analysis of current status data, many parametric and semi-parametric methods have been proposed. In particular, many of these efforts have focused on the proportional hazards (Cox,, 1972) and proportional odds models (Dabrowska and Doksum,, 1988); for a review see McMahan et al., (2013), Lu and Li, (2017), and the references therein. To provide a more general modeling framework, Gu et al., (2005), Sun and Sun, (2005), Zhang et al., (2005), Zhang and Zhao, (2013), and Zeng et al., (2016) developed estimation techniques under linear transformation models. A potentially prohibitive assumption made in all of the aforementioned works is that the covariate effects on the failure time are “linear.” To relax this assumption, Cheng et al., (2011) developed a semi-parametric additive transformation model and Lu and McMahan, (2018) presented a partially linear PH model for current status data. For a comprehensive review of methods to analyze interval-censored data see Sun, (2007) and Banerjee, (2012).

In this article, we seek to develop a novel methodology which can be used to conduct the regression analysis of current status data. In particular, a very general partially linear additive transformation model is proposed. The proposed model makes use of constrained $B$-splines to approximate all unknown functions; i.e., the transformation and regression functions. In order to provide for more efficient estimators, a penalization strategy is utilized to regularize the estimation of the unknown functions. To complete model fitting a computationally efficient and easy to implement hybrid algorithm is developed. Theoretically, it is shown that the penalized estimators of the regressions coefficients are asymptotically normal and efficient; i.e., they attain the semi-parametric information bound. Moreover, we establish that the estimators of the nonparametric components achieve the optimal rate of convergence for an appropriately chosen order of the smoothing parameters. Finite sample techniques of determining the smoothing parameter and estimating the asymptotic covariance matrix of the estimated regression coefficients are provided and evaluated.

2 Model

Consider a study which examines $n$ subjects for a failure time of interest. Let $T_{i}$ denote the time that the $i$th subject experiences the failure, for $i=1,....,n$. To provide for modeling flexibility, we assume that $T_{i}$ obeys a partial linear additive transformation model which is given by

[TABLE]

where $\eta(\cdot)$ is an unknown monotone increasing transformation function on $(0,\infty)$ with $\eta(0)=-\infty$ and $\varepsilon_{i}$ is a random variable having cumulative distribution function (CDF) $H(\cdot)$. In the above specification, $\beta$ is a $q$-dimensional vector of regression coefficients corresponding to the covariates ${Z}_{i}=(Z_{i1},...,Z_{iq})^{\intercal}$ and $\varphi_{j}(\cdot)$ is an unknown smooth regression function taking the covariate $W_{ij}$ as its argument, for $j=1,...,J$. Throughout, we assume that $\varepsilon_{i}\perp\varepsilon_{j}$, for $i\neq j$, and $\varepsilon_{i}\perp{X}_{i}$, where $\perp$ denotes statistical independence and ${X}_{i}=({Z}_{i}^{\intercal},W_{i1},...,W_{iJ})^{\intercal}$. This formulation provides for a very general modeling construct that holds several popular models as special cases. For example, specifying $H(\cdot)$ to be the CDF of the extreme value distribution causes (1) to reduce to the standard proportional hazards (PH) model, while taking $H(\cdot)$ to be the CDF of the standard logistic distribution results in the proportional odds (PO) model.

Following the work of Murphy et al., (1997), Sun and Sun, (2005), and Cheng et al., (2011), we note that (1) can be rewritten in terms of the conditional CDF of $T_{i}$, given the covariate vector ${X}_{i}$, and can be expressed as

[TABLE]

where $\textsl{g}(\cdot)$ is a known, smooth, and strictly increasing link function on $(0,1)$. Though equivalent to (1), the representation provided in (2) is more convenient for reasons that will shortly become apparent. Here $\textsl{g}(\cdot)$ assumes the role of $H(\cdot)$ with respect to determining the final form of the proposed model; e.g., specifying $\textsl{g}(x)=\log\{-\log(1-x)\}$ provides for the PH model and $\textsl{g}(x)=\log\{x/(1-x)\}$ results in the PO model.

A key feature of current status data is that the failure times (i.e., $T_{i}$, for $i=1,\dots,n$) are not directly observed but rather are known relative to observation times; i.e., the $T_{i}$ is either left- or right-censored indicating that $T_{i}$ occurs before or after the observation time, respectively. To mathematically frame the structure of this data, let $Y_{i}$ denote the observation time for the $i$th individual and let $\Delta_{i}=I(T_{i}<Y_{i})$ denote the corresponding censoring indicator, where $I(\cdot)$ is the usual indicator function. Note, $\Delta_{i}=1(0)$ represents the event that the $i$th individual’s failure time is left(right)-censored. Thus, the observed data consist of $\{(Y_{i},\Delta_{i},{X}_{i})\}_{i=1}^{n}$. Following the works of Mammen and van de Geer, (1997), Ma et al., (2005), and Ma, (2009), we propose to estimate the unknown parameters via the following penalized log-likelihood

[TABLE]

where $\tau=(\beta,\eta,\varphi_{1},...,\varphi_{J})$ represents the collection of unknown model parameters, $\lambda=(\lambda_{0},\ldots,\lambda_{J})^{\intercal}$ is a vector of smoothing parameters, such that $\lambda_{j}\geq 0$ for all $j$, and

[TABLE]

In the penalties above, $\eta^{(r)}$ and $\varphi_{j}^{(r)}$ denote the $r$th derivative of $\eta(\cdot)$ and $\varphi_{j}(\cdot)$, respectively, $r\geq 1$. This penalized log-likelihood is derived under the assumption that the covariates are time-invariant and that the failure and censoring times are conditionally independent, given the covariates. These assumptions are common among the literature; e.g., see Ma et al., (2005), Cheng et al., (2011), Banerjee, (2012), and the references therein.

3 Large sample properties

Let $\tau_{0}=(\beta_{0},\eta_{0},\varphi_{10},...,\varphi_{J0})$ and $\widehat{\tau}=(\widehat{\beta},\widehat{\eta},\widehat{\varphi}_{1},...,\widehat{\varphi}_{J})$ denote the true value and the penalized maximum likelihood estimator of $\tau$, respectively. We established the asymptotic properties of the penalized estimators $\widehat{\tau}$ under the following regularity conditions:

Condition 1

The true parameter $\beta_{0}$ is an interior point of a compact subset of $\mathbb{R}^{q}$.

Condition 2

The $r$th derivatives of $\eta_{0}$ and $\varphi_{0j}$, for $1\leq j\leq J$, satisfy the Lipchitz condition on $[c,d]$ for $c>0$ and $[a_{j},b_{j}]$, respectively, and $\eta_{0}$ is strictly increasing on $[c,d]$.

Condition 3

The support of the observation times is an interval within $[l_{y},u_{y}]$ with $0<l_{y}<u_{y}<\inf\{t:F_{0}(t)=1\}$, where $F_{0}$ is the CDF of the observation times.

Condition 4

The covariate vector ${Z}$ has bounded support.

Condition 5

The first derivative of the link function g is bounded away from 0 and the second derivative of g is bounded on (0,1).

Condition 6

The joint density of $(Y,{Z},{W},\Delta)$ is bounded away from 0 and $\infty$.

Condition 7

For any $\beta\neq\beta_{0}$, $\Pr({Z}^{\intercal}\beta\neq{Z}^{\intercal}\beta_{0})>0$ and $E\{\varphi_{j}(W_{j})\}=0$,   for $1\leq j\leq J$.

Condition 8

The smoothing parameters are specified such that

[TABLE]

Condition 9

The efficient information $\mathscr{I}(\beta_{0})$ defined in Appendix 1 of the Supplementary Material is non-singular.

Remark 1

Condition 1 is a standard assumption in semi-parametric estimation. The smoothness assumption of $\eta_{0}$ and $\varphi_{0j}$ in condition 2 guarantees that the unknown nonparametric components can be well approximated by $B$-splines. Conditions 3–6 are necessary in entropy calculation used to derive the rate of convergence of $\widehat{\tau}$ and the asymptotic normality of $\widehat{\beta}$. Condition 7 is required to establish the identifiability of the proposed semi-parametric model. Condition 8 is a typical assumption used to derive the asymptotic property of penalized estimators. Finally, condition 9 is common for the asymptotic normality and the efficiency of $\widehat{\beta}$.

Using modern empirical process theory, we establish the following asymptotic properties:

Theorem 1

(Consistency and rate of convergence)  Under conditions 1–8, the penalized estimator $\widehat{\beta}$ is asymptotically consistent for $\beta_{0}$, $\|\widehat{\eta}\|_{\infty}=\|\widehat{\varphi}_{j}\|_{\infty}=O_{p}(1)$, and $\|\widehat{\eta}-\eta_{0}\|_{2}=\|\widehat{\varphi}_{j}-\varphi_{j0}\|_{2}=O_{p}(n^{-r/(1+2r)})$, for $1\leq j\leq J$.

Theorem 2

(Asymptotic normality and efficiency) Under conditions 1–9,

[TABLE]

The proofs of these results are relegated to Appendix 1 of the Supplementary Material. It is worthwhile to point out that if the smoothing parameters are chosen to be of the right order (see condition 8), our functional estimates attain the optimal rate of convergence; i.e., $O_{p}(n^{r/(1+2r)})$. Moreover, $\widehat{\beta}$ is efficient in the semi-parametric sense since it achieves the information bound.

4 Estimation

4.1 Penalized spline log-likelihood

In general, estimating $\beta$, $\eta$ and $\varphi_{j}$ directly from the penalized log-likelihood can be a tumultuous task. Thus, motivated by Xiang and Wahba, (1997), Ma et al., (2005), and Lu and Li, (2017), we propose to strike a balance between modeling flexibility and complexity by approximating all unknown functions in (2) through the use of constrained $B$-splines functions. In particular,

[TABLE]

where $b_{k}(\cdot)$ and $b_{jk}(\cdot)$ are $B$-spline basis functions with unknown coefficients $\gamma=(\gamma_{1},\ldots,\gamma_{p_{0}})^{\intercal}$ and $\alpha_{j}^{*}=(\alpha_{j1}^{*},...,\alpha_{jp_{j}}^{*})^{\intercal}$, respectively, $1\leq j\leq J$. Once a knot set is chosen and the degree of the spline function is specified these basis functions are uniquely determined; for further discussion see Schumaker, (2007). Constraints are placed on the spline coefficients to insure the monotonicity of the transformation function and the identifiability of the regression functions. The monotonicity of $\eta_{n}(\cdot~{}|\gamma)$ can be insured by requiring that $\gamma_{1}\leq\ldots\leq\gamma_{p_{0}}$; see Theorem 5.9 of Schumaker, (2007). As with all additive models, $\varphi_{j}(\cdot)$ is only identifiable up to an additive constant. Thus, sum-to-zero constraints are imposed such that $\sum_{i=1}^{n}\varphi_{nj}(W_{ji}\mid\alpha_{j}^{*})=0$, for $1\leq j\leq J$.

For implementation purposes, the sum-to-zero constraints can be expressed as ${B_{j}^{*}}^{\intercal}\alpha_{j}^{*}=0$, where $B_{j}^{*}=\sum_{i=1}^{n}B_{ji}^{*}$, $B_{ji}^{*}=\{b_{j1}(W_{ji}),\ldots,b_{jp_{j}}(W_{ji})\}^{\intercal}$, and $\alpha_{j}^{*}=(\alpha_{j1}^{*},...,\alpha_{jp_{j}}^{*})^{\intercal}$. Motivated by Wood, (2017), the model is reparameterized such that the constrained set of parameters $\alpha_{j}^{*}$ can be uniquely determined by an unconstrained reduced set $\alpha_{j}=(\alpha_{j1},...,\alpha_{jp_{j}-1})^{\intercal}$. In particular, let $Q_{j}$ denote a $p_{j}\times(p_{j}-1)$ semi-orthogonal matrix whose columns are orthogonal to $B_{j}^{*}$. Obviously, ${B_{j}^{*}}^{\intercal}Q_{j}\alpha_{j}=0$, for all $\alpha_{j}\in\mathbb{R}^{p_{j}-1}$, and for any $\alpha_{j}^{*}$ satisfying the sum-to-zero constraint there exists $\alpha_{j}$ such that $\alpha_{j}^{*}=Q_{j}\alpha_{j}$. Note, $Q_{j}$ can easily be found via a QR-decomposition of $B_{j}^{*}$. Thus, the log-likelihood under the proposed spline model is obtained by replacing $F(Y_{i}\mid X_{i})$ in the penalized log-likelihood function by

[TABLE]

where ${M}_{i}=\{b_{1}(Y_{i}),\ldots,b_{p_{0}}(Y_{i})\}^{\intercal}$ and ${B}_{ji}=Q_{j}^{\intercal}B_{ji}^{*}$, for $1\leq i\leq n$. Thus, the set of unknown parameters to be estimated are given by $\theta=(\beta^{\intercal},\gamma^{\intercal},\alpha^{\intercal})^{\intercal}$, where $\alpha=(\alpha_{1}^{\intercal},\ldots,\alpha_{J}^{\intercal})^{\intercal}$.

To further alleviate computational burdens, we follow the proposal of Eilers and Marx, (1996) and replace the penalty terms $\mathcal{J}^{2}(\eta)$ and $\mathcal{J}^{2}(\varphi_{j})$ by discrete approximations. Under this strategy, the penalized spline log-likelihood function is given by

[TABLE]

where ${D}_{j}^{(r)}$ is the discrete difference operator of order $r$ (Eilers and Marx,, 1996; Tibshirani et al.,, 2014) and $\|\cdot\|$ is the Euclidean norm. This operator is best defined recursively as ${D}_{j}^{(r+1)}={D}_{j}^{(1,r)}{D}_{j}^{(r)}$, where

[TABLE]

where ${D}_{j}^{(1)}={D}_{j}^{(1,0)}$. Proceeding in this fashion for $r=1$ and $2$ yields penalties of the form

[TABLE]

respectively. Through extensive numerical studies, see Section 5, we have found that setting $r=2$ is sufficient when implementing the proposed approach.

Due to these approximations, the proposed penalized spline estimator is different from the penalized estimator, but as long as the number of basis functions for each of the unknown functions grows at least at the rate $O(n^{1/5})$, this difference will vanish and the statistical performance of the two estimators should be indistinguishable; for further discussion see Eilers and Marx, (1996), Ma et al., (2005), Wood, (2017), and the references therein. That is, the penalized spline estimator should inherit the theoretical properties of the penalized maximum likelihood estimator. This assertion is supported by the results of the numerical experiments reported in Section 5. To ensure the accuracy of the approximation, it is generally advised that the number of knots for each function be chosen in a neighborhood of $n^{1/3}$; for further discussion see Section 5.

4.2 Parameter estimation

To complete model fitting a computationally efficient hybrid algorithm is developed which can be used to identify both the penalized spline estimators and the smoothing parameters via a nested iterative process. The inner process is aimed at identifying the value of $\theta$ that maximizes (4.1) for a fixed value of $\lambda$, while adhering to the monotonicity constraints; i.e., the inner step identifies $\widetilde{\theta}_{\lambda}=\operatorname*{\arg\!\max}_{\theta}l_{n,\lambda}(\theta)$ s.t. $\gamma_{1}\leq\ldots\leq\gamma_{p_{0}}$. To solve this constrained optimization problem, an iterative algorithm is developed. In each iteration of this algorithm, the unconstrained updates of the parameters, say $\theta^{*}$, are first determined based on the current parameter values, say $\theta^{(m)}$, via the following Fisher-scoring step

[TABLE]

where $\nabla l_{n,\lambda}(\theta^{(m)})$ and $\mathcal{I}_{n,\lambda}(\theta^{(m)})$ are the gradient and the expected value of the negative hessian matrix of (4.1) evaluated at $\theta^{(m)}$, respectively. Closed form expressions for these quantities are given by

[TABLE]

where $\mathcal{X}=(Z,M,B_{1},...,B_{J})$, ${Z}^{\intercal}=({Z}_{1},\ldots,{Z}_{n})$, ${M}^{\intercal}=({M}_{1},\ldots,{M}_{n})$, and ${B}_{j}^{\intercal}=({B}_{j1},\ldots,{B}_{jn})$, and ${S}=\sum_{j=0}^{J}{\lambda_{j}^{2}}{S}_{j}$, with $S_{j}=\textrm{diag}\{0_{r_{1j}\times r_{1j}},{D}_{j}^{(r)^{\intercal}}{D}_{j}^{(r)},0_{r_{2j}\times r_{2j}}\}$, $r_{1j}=q+p_{0}+...+p_{j-1}-(j-1)$, and $r_{2j}=p_{j+1}+...+p_{J}-(J-j)$. In the expressions above, the transformed response vector and weight matrix for a specific value of $\theta$ are given by

[TABLE]

where $\mathcal{X}_{i}$ is the $i$th row of $\mathcal{X}$, $G^{\prime}(\cdot)$ is the derivative of $G(\cdot)=\textsl{g}^{-1}(\cdot)$, and $\pi_{i}=\textsl{g}^{-1}(\mathcal{X}_{i}^{\intercal}\theta)$.

The updated values of $\gamma$ available in $\theta^{*}$ do not necessarily satisfy the monotonicity constraints; see Section 2. To enforce these conditions, isotonic regression is implemented to project the unconstrained estimates into the constrained space; i.e., into $\mathscr{G}=\{\gamma:\gamma_{1}\leq\ldots\leq\gamma_{p_{0}}\}$. This is accomplished by solving the following minimization problem

[TABLE]

where $\sigma_{k}^{2}$ is the diagonal element of $\mathcal{I}^{-1}_{n,\lambda}(\theta^{*})$ corresponding to $\gamma_{k}$, for $k=1,...,p_{0}$. This problem can be solved using the pava function in the R package Iso, which takes $\gamma^{*}$ and $(\sigma_{1}^{2},\ldots,\sigma_{p_{0}}^{2})^{\intercal}$ as inputs. Substituting $\gamma^{(m+1)}$ for $\gamma^{*}$ in $\theta^{*}$ one obtains the constrained update $\theta^{(m+1)}$. These two steps are completed in turn until convergence is attained; i.e., the difference between $\theta^{(m)}$ and $\theta^{(m+1)}$ is less than some specified stopping criterion. At the point of convergence of this inner process, one obtains $\widetilde{\theta}_{\lambda}$; i.e., the maximizer of the penalized spline log-likelihood function for a fixed value of $\lambda$. Once the inner process is complete, the outer process updates the smoothing parameters and the inner process is restarted. To obtain the updated smoothing parameters, say $\lambda^{*}=(\lambda_{0}^{*},...,\lambda_{J}^{*})^{\intercal}$, a generalized Fellner-Schall approach (Wood and Fasiolo,, 2017) is adopted and makes the update as

[TABLE]

where ${S}^{-}$ is the generalized inverse of ${S}$. Using the updated smoothing parameters (i.e., $\lambda^{*}$), the inner process is then utilized to identify $\widetilde{\theta}_{\lambda^{*}}$. This nested iterative process continues until the difference between $\widetilde{\theta}_{\lambda}$ and $\widetilde{\theta}_{\lambda^{*}}$ is less than some specified stopping criterion. It is worth while to point out that the update of $\lambda_{j}$ is guaranteed to be positive, which is required, as long as $\mathcal{I}_{n,\lambda}(\widetilde{\theta}_{\lambda})$ is positive definite (Wood and Fasiolo,, 2017), which is the case for the proposed model. At convergence of the proposed hybrid algorithm, let $\widetilde{\theta}=(\widetilde{\beta}^{\intercal},\widetilde{\gamma}^{\intercal},\widetilde{\alpha}^{\intercal})^{\intercal}$ denote the final update of $\theta$ with $\widetilde{\eta}(\cdot)=\eta_{n}(\cdot~{}\mid\widetilde{\gamma})$ and $\widetilde{\varphi}_{j}(\cdot)=\varphi_{nj}(\cdot~{}\mid\widetilde{\alpha}_{j})$ being the proposed penalized spline estimators of $\eta(\cdot)$ and $\varphi_{j}(\cdot)$, respectively.

4.3 Variance estimation

For the purposes of conducting large-sample inference an estimator of the variance-covariance matrix of $\widetilde{\beta}$ was developed. This estimator is given by the upper $q\times q$ submatrix of

[TABLE]

By applying the least-square based variance estimation approach (Zhang et al.,, 2010), it can be shown that the proposed estimator consistently estimates $\mathscr{I}^{-1}(\beta_{0})$. Using this estimator one can conduct standard Wald type inference about the regression coefficients.

5 Simulation Study

To assess the finite sample performance of the proposed methodology the following simulation study was conducted. In this study, the true failure time $T$ was generated according to the following model

[TABLE]

where $W_{j}\sim\mathcal{U}[-1,1]$, $Z_{1}\sim\text{Bernoulli}(0.5)$ and $Z_{2}\sim\mathcal{N}(0,1)$. To provide for a broad range of examples, the link function $\textsl{g}_{\alpha}(\cdot)$ was assumed to belong to a family of links indexed by a parameter $\alpha$ as:

[TABLE]

Proceeding in this fashion allows us to capture several common survival models; i.e., $\alpha=0$ leads to the PH model while $\alpha=1$ corresponds to a PO model. To examine different forms of the unknown functions and regression coefficients three scenarios were considered. In Scenario 1 (S1) the three functions were specified as $\eta(t)=\log(2t)$, $\varphi_{1}(w)=\exp(w+0.5)-\{\exp(1.5)-\exp(-0.5)\}/2$, and $\varphi_{2}(w)=2\sin(-\pi w)$ and the regression coefficients were chosen to be $\beta_{1}=0.5$ and $\beta_{2}=-0.5$; in Scenario 2 (S2) the three functions were specified as $\eta(t)=\log\{1.5t-\log(1+1.5t)\}$, $\varphi_{1}(w)=2\sin(-\pi w)$, and $\varphi_{2}(w)=4w^{2}-4/3$ and the regression coefficients were set to be $\beta_{1}=0.5$ and $\beta_{2}=0.5$; in Scenario 3 (S3) the three functions were specified as $\eta(t)=\log\{\log(1+t/10)+\sqrt{t}/10\}$, $\varphi_{1}(w)=4w^{2}-4/3$, and $\varphi_{2}(w)=\exp(w+0.5)-\{\exp(1.5)-\exp(-0.5)\}/2$ and the regression coefficients were set to be $\beta_{1}=-0.5$ and $\beta_{2}=-0.5$. In each of S1-S3, we consider values of $\alpha\in\{0,0.5,1\}$. Note, each of the specified functions adheres to the identifiability constraint; i.e., they integrate to zero over the support of the corresponding covariates. To create current status data, observation times ($Y$) were sampled from an exponential distribution with a mean of 2, 2, and 1 under S1, S2, and S3, respectively, and the censoring indicator was determined as $\Delta=I(T<Y)$. These specifications were made to yield a variety of right-censoring rates; for example when $\alpha=0(1)$ we have right censoring rates of 27%(36%), 44%(51%), and 77%(81%) under S1, S2, and S3, respectively. For each simulation configuration, this data generating process was repeated 500 times to create 500 independent data sets consisting of $n=400$ observations, which is approximately one fourth the sample size available in our motivating data application.

The proposed methodology was used to analyze each of the simulated data sets. For purposes of comparison, we compare the results of the proposed approach to those of the spline estimator that arises from maximizing (4.1) after dropping the penalty term. Note, this competing approach is a generalization of the technique outlined by Lu and McMahan, (2018) which was shown to outperform Cheng et al., (2011). To compute both the penalized and unpenalized spline estimators, cubic $B$-spines were used to approximate all unknown functions with interior knots being placed at equally-spaced quantiles of the covariate values. For each data set, the proposed approach made use of a single knot set consisting of $\lceil n^{1/3}\rceil$ interior knots, while the unpenalized estimator was computed using multiple knot set configurations with the number of interior knots ranging from $\lceil n^{1/3}\rceil-3$ to $\lceil n^{1/3}\rceil+3$, with the final unpenailzed estimator being selected based on the BIC criterion; for further discussion see Lu and McMahan, (2018).

Table 1 summarizes the results of estimating the regression coefficients across all considered simulation configurations for both the penalized and unpenalized spline estimator. Provided in this summary are the empirical bias, sample standard deviation, and mean squared error of the 500 point estimates, as well as the average of the 500 estimated standard errors and empirical coverage probability associated with 95% confidence interval. From these results it is apparent that the proposed approach works well; i.e., the penalized spline estimators of the regression coefficients exhibit little bias, the standard deviations of the point estimates are in agreement with the average estimated standard errors, and the empirical coverage probabilities are at their nominal level. The latter two findings tend to suggest that the proposed method of estimating the variance-covariance matrix works well and that large-sample inference is possible even for relatively small sample sizes. Further, the proposed penalized spline estimator outperformed in virtually every aspect its unpenalized counterpart even though this competing estimator was chosen as the best from amongst multiple candidates. Moreover, the performance of this competing technique varied dramatically across the different knot set configurations; i.e., the unpenalized spline estimator was found to be sensitive to the specification of the number and placement of the basis functions. However, penalized estimator was robust to the specification of the basis functions.

One of the primary goals of this study was to assess the performance of the proposed method with respect to estimating the unknown functions. To this end, Figure 1 summarizes the estimates of the unknown functions for S2 for all considered values of $\alpha$ that were obtained from both approaches. This summary includes the 0.025, 0.5, and 0.975 point-wise quantiles of the estimates. The analogous figures under S1 and S3 can be found in Appendix 2 of the Supplementary Material. These results suggest that the proposed methodology can be used to accurately estimate all of the unknown functions; i.e. very little bias can be observed when one compares the the median of the functional estimates to the true underlying function. Additionally, these findings demonstrate that the penalized spline estimator results in functional estimates that have smaller variability when compared to the unpenalized spline estimator. These finding reinforce the assertion that the proposed approach performs better with respect to both estimation and inference than its unpenalized counterpart, and moreover, that it is robust to the specification of the spline functions.

6 Application to Uterine leiomyomata data

For further illustration, the proposed methodology is now used to analyze uterine leiomyomata (fibroid) data collected as a part of a prospective cohort study concerned with early pregnancy. Uterine fibroids are noncancerous growths that appear during childbearing years and have been found to be associated with complications in pregnancy and delivery (Coronado et al.,, 2000). As a part of this study, a patient’s fibroid status was determined during an ultrasound examination; i.e., during this examination fibroids were either found to be present or absent indicating that the fibroid onset time was smaller or larger than the examination time, respectively. Thus, the fibroid onset time is either left- or right-censored with the patient’s age at examination being the censoring time. Of the $n=1604$ participants considered in this analysis 216 were found to be fibroid positive at the time of their examination, leading to a left-censoring rate of approximately 13%.

In this analysis, we focus on four risk factors purported to be related to fibroid development; i.e., race ($Z_{1}$ = 1 if African American and $Z_{1}$ = 0 Caucasian), parity status ($Z_{2}$ = 1 if they had given birth before and $Z_{2}$ = 0 otherwise), obesity status ($Z_{3}$ = 1 indicating obesity and $Z_{3}$ = 0 otherwise), and standardized age of menarche ($W$). To relate these risk factors to the fibroid onset time, the following model was specified

[TABLE]

where $\textsl{g}(u)=\log\{-\log(1-u)\}$. This link function leads to the usual PH model. Cubic $B$-splines with a knot set consisting of $\lceil n^{1/3}\rceil=12$ interior knots were used to approximate all unknown functions. Model fitting was completed using the hybrid algorithm outlined in Section 4. The estimator of the variance-covariance matrix of $\widetilde{\beta}$ developed in Section 4 was used to provide for standard errors and $95\%$ Wald confidence intervals for the estimated regression coefficients.

Table 2 summarizes the regression parameter estimates, along with their estimated standard errors and corresponding 95% confidence intervals. The results suggest that race and parity are significant factors, while obesity is not. In particular, these results suggest that African American women have a higher hazard of developing fibroids when compared to Caucasian women, likewise it appears that women who have not given birth before appeared to have a higher hazard of developing fibroids when compared to those that had. These results and conclusions reinforce those presented in Wang and Dunson, (2011) and Lu and McMahan, (2018). Figure 2 provides the estimates of $\eta(\cdot)$ and $\varphi(\cdot)$, along with 95% point-wise confidence intervals. The upper and lower limits of the 95% point-wise confidence intervals were obtained as the 0.025 and 0.975 point-wise quantiles of 1000 bootstrap estimates. The estimate of $\varphi(\cdot)$ clearly exhibits a nonlinear pattern; i.e., standardized age of menarche has a nonlinear impact on the development of fibroids. To examine the effect of the assumed model form, the aforementioned analysis was conducted under the PO model by specifying $\textsl{g}(x)=\log\{x/(1-x)\}$. The estimates of the regression coefficients under the PO model are presented in Table 2. These results provide no contradiction to the results obtained under PH model. A figure summarizing the estimates of $\eta(\cdot)$ and $\varphi(\cdot)$ under the PO model, along with 95% point-wise confidence intervals, is provided in Appendix 2 of the Supplementary Material, and these estimated functions exhibit similar features to those obtained under PH model.

7 Summary and discussion

In this work we developed a computationally efficient penalized approach for fitting flexible partially linear additive transformation models to current status data. The proposed procedure demonstrates appealing properties, both numerical and theoretical. The hybrid algorithm developed for model fitting is easy to implement and robust to initialization. Through the use of penalization, the proposed approach attains superior finite-sample performance when compared to the standard spline alternative. A simple variance-covariance estimation procedure was established and was shown to provide reliable Wald-type inference. In addition to attractive finite sample performance, the estimators of the nonparametric components are shown to attain the optimal rate of convergence and we prove that estimators of regression coefficients are asymptotically normal and efficient in the semi-parametric sense.

Given the successes of this approach, future work will be directed at generalizing this methodology to allow for the analysis of interval-censored data. A key obstacle in the development of this methodology will be the model fitting strategy, which will be far more complicated due to the complex nature of interval-censored data. Another direction for future research could involve further extending this work to allow for multiple end points, both for current status and interval censored data. Lastly, and potentially most challenging, would be the development of goodness-of-fit tests that could be used to guide the specification of the link function. This is a notoriously hard problem in the transformation model setting, and is only made harder when one considers censored data, and/or additive models.

Acknowledgement

This research was supported by the National Institutes of Health grant AI121351 and National Science Foundation grant OIA-1826715.

Supplementary material

Technical details and proofs of Theorems 1 and 2 are provided in the Supplementary Materials, along with the additional results referenced in Sections 5 and 6.