Semiparametric Odds Rate Model for Modeling Short-Term and Long-Term Effects with Application to a Breast Cancer Genetic Study

1 Introduction

The proportional odds (PO) model [1–3] is commonly used for analyzing survival data. The PO model specifies that the survival function of the failure time T given the covariates X takes the form

where H(t) is a strictly increasing function with H(0)=0, and β is a vector of regression parameters. Unlike the Cox proportional hazards model [4], the hazard ratio between two sets of covariates is time-varying under the PO model. Instead, the odds ratio of survival is assumed to be constant over time under the PO model. Consequently, the hazard ratio converges to one as time goes to infinity.

Much attention has been paid to the development of inference procedures under the PO model and its generalizations. Murphy et al. [5] provided statistical inference for the maximum likelihood estimation with right-censored data under the PO model and demonstrated that the estimators of the regression parameters provided by Bennett [1] are efficient. The sieve maximum likelihood estimation and the estimation based on weighted empirical odds functions were developed by Shen [6] and Yang and Prentice [7], respectively. Computational issues in the PO models are discussed in Hunter and Lange [8]. The PO model has also been extensively studied for current status data and interval censored data, see Rossini and Tsiatis [9], Huang and Rossini [10], Rabinowitz et al. [11] and Sun et al. [12], among others. For correlated failure time data, the PO model with random effects has been investigated in literature [13, 14].

All the aforementioned methods are developed under the assumption of constant odds ratio over time. However, in some applications, this assumption is violated. Violation of this assumption may lead to biased and misleading results. A motivating example is from a copy number alteration (CNA) study with breast cancer patients [15]. The objective of the study was to determine whether the CNAs influence the time to develop distant metastasis among breast cancer patients. We expect the log odds curves of survival to be parallel for two groups if the PO assumption is true. Figure 1 presents the Kaplan–Meier curves for two groups partitioned by the median of copy number variation data at two SNP loci. It appears that the PO assumption is reasonable at SNP #2124; however, the two curves cross at SNP #39687 indicating a violation of the PO assumption. The PO model failed to capture the phenomenon of the crossing survival at SNP #39687.

Figure 1

The Kaplan–Meier survival curves of time to develop distant metastasis for the CNA study at SNP #2124 and SNP #39687

To test the PO assumption, Dauxois and Kirmani [16] developed a formal analytical test for two-sample problems. Resampling approaches using observed score processes for the regression coefficients have been developed by Scheike [17] to check the goodness-of-fit of the PO model. Liu et al. [18] provided graphical diagnostics to check model misspecification for the PO model. To alleviate the PO assumption, Scheike [17] proposed a flexible transformation model which allows for time-varying coefficients for some covariates in the model. They developed estimation and inference procedures based on martingale estimating equations. Additionally, they proposed to use a robust sandwich-type covariance matrix for the estimates of the regression coefficients. More recently, Chen, Hu, Cheng, Musoke, and Zhao [19] proposed an extended PO model by incorporating “external” time-varying covariates and developed appropriate inference procedures for the regression parameters. Note that both Liu et al. [18] and Chen et al. [19] do not incorporate time-varying coefficients.

In this article, we propose an alternative approach for relating covariates of interest to the odds of survival for right-censored data. Specifically, we propose a novel semiparametric odds rate model to incorporate the short-term and long-term covariate effects. This model includes the PO model as a special case, leading to a new diagnostic tool for checking the PO assumption. We establish efficient likelihood-based estimation and inference procedures and derive the asymptotic properties of the proposed estimators. Results from numerical studies demonstrate that the proposed methods perform well in practical settings.

The remaining of the article is organized as follows. In Section 2, we propose the new model and derive the nonparametric maximum likelihood estimation. The model assumptions and the large sample properties are presented in Section 3. In Section 4, we conduct extensive simulation studies to investigate the finite sample properties of the proposed estimators. An application to the motivating CNA study is provided in Section 5 and we conclude with a brief discussion in Section 6. Theoretical proofs are provided in the Appendix.

2 Method

Suppose there are n independent subjects in a random sample. For the ith subject, let Ti be the failure time, Ci be the censoring time, and Xi be a p×1 vector of time-independent covariates. The observed data contain {(Yi=min(Ti,Ci),Δi=I(Ti≤Ci),Xi),i=1,…,n}, where I(⋅) is the indicator function.

Let F(t|X) and S(t|X) be the conditional distribution function and survival function of the failure time T given the covariates X, respectively. We propose the following generalized odds model

where H(t) is a non-decreasing function with H(0)=0. In model (1), H(t) is the baseline odds function of failure, S(t) is the baseline survival function which is the survival function for a subject with X=0, and F(t)=1−S(t) is the baseline cumulative distribution function.

Under the proposed model (1), the regression coefficients β and γ have attractive and intuitive interpretations. For two sets of covariates X and X∗, we can show that

Therefore, β and γ can be interpreted as the long-term and short-term log-odds ratios, respectively. When β=γ, model (1) reduces to the PO model. Hence, this model provides a useful diagnostic tool for testing the PO model. To be more specific, we can check the PO model by testing the null hypothesis H0:β=γ. The test can be performed by using Wald, score or likelihood ratio test statistics.

We now derive the likelihood function for the unknown parameters (β,γ,H). It can be shown that

It follows that

Consequently, the hazard function is given by

where H′(t) is the first derivative of H(t). Under the conditional independent censoring assumption, the likelihood function for (β,γ,H) has the form

We would like to make inference based on the likelihood function. Note that the maximum of the likelihood function does not exist since we can always choose H′(Yi)=0 for some Yi with Δi=1. Therefore, following the techniques of nonparametric maximum likelihood estimation in the Cox proportional hazards model [20, 21], the PO model [5, 14], and the semiparametric transformation model [22], we allow H to be a right continuous function. Therefore, we allow the estimator of H to be a step function with positive jumps only at the observed failure time points. Let H{Yi} denote the jump size of H(t) at Yi. We replace H′(Yi) with H{Yi} in eq. (2) to obtain the nonparametric likelihood function, still denoted by Ln(β,γ,H) for simplicity,

Assume the data have m unique observed failure time points Y(1)<Y(2)<⋯<Y(m). Easily we have H(t)=∑k:Y(k)≤tH{Y(k)}. Then

In this case, the odds function H(t|X)=F(t|X)/S(t|X) is given by

Hence, to estimate the unknown parameters, we maximize the nonparametric likelihood function with respect to the regression parameters β and γ, and jump sizes H{Y(1)}, H{Y(2)}, …, H{Y(m)}. The resulting nonparametric maximum likelihood estimators (NPMLEs) are denoted by (βˆn, γˆn, Hˆn).

The maximization procedure can be performed by using the quasi-Newton optimization algorithm in Press et al. [23]. Based on Fisher information theory, the covariance matrix for (βˆn, γˆn, Hˆn) can be obtained by inverting the observed information matrix, which is the negative second derivatives of the nonparametric likelihood function with respect to the regression parameters β and γ, and the jump sizes H{Y(1)}, H{Y(2)}, …, H{Y(m)}. The justification follows the arguments in theorem 2 of Parner [24].

3 Asymptotic properties

Before we derive the asymptotic properties of the NPMLEs, we impose the following regularity conditions:

1. There exists some positive constant δ0 such that Pr(Ci≥τ|Xi)≥δ0 almost everywhere, where τ is the end of the study. In addition, Pr(Ti>τ|Xi)>0.

2. Conditional on the vector of covariates Xi, the censoring time Ci is independent of Ti for i=1,…,n.

3. Xi are bounded for all i=1,…,n. Furthermore, aTXi=0 almost everywhere if and only if a=0.

4. The true values β0 and γ0 for parameters β and γ belong to the interior of a known compact set B0 in R2p.

5. The true baseline odds function H0 belongs to

   H0={H:Hisstrictlyincreasingin[0,τ]andiscontinuouslydifferentiablewithH(0)=0andH(τ)<∞}.

We first establish the consistency of the NPMLEs.

The existence of (β,γ) can be proved from the compactness of the parameter space and the boundedness of covariates. Suppose for some Yi with Δi=1, we have Hˆn{Yi}=∞. In that case, the likelihood function must be 0, which implies Hˆn{Yi} must be finite. Therefore, the NPMLEs exist. The detailed proof for Theorem 1 is provided in Appendix B. The key part of the proof is to establish the boundedness of Hˆn(τ). Then we can construct a limit, (β∗,γ∗,H∗), for the sequences (βˆn,γˆn,Hˆn) and prove that the limit is actually (β0,γ0,H0) by using the Kullback–Leibler information.

We next establish the asymptotic normality of the NPMLEs. Let l∞[0,τ] be a space consisting of all uniformly bounded functions on [0,τ] and the norm on it is defined as the supremum absolute value on [0,τ].

4 Simulation studies

We conducted simulation studies to examine the finite-sample performance of our proposed methods. We generate the failure time from the following model

where X was set to be uniform(−0.5,0.5). Under the above model, H(t)=t.

Results were obtained under four scenarios of regression parameters (β,γ): (0.5, 0.5), (0, 0.5), (–0.5, 0.5), and (–0.5, 0). The long-term and short-term effects are equal under the first scenario where the model reduces to the PO model; there is no long-term effect under scenario 2; the long-term and short-term effects are in opposite directions under scenario 3; and there is no short-term effect under scenario 4. We generated the censoring time from an exponential distribution with mean 6. The censoring time was truncated at 4. The censoring rates were about 30% under all scenarios. For each setting, we considered the sample sizes of 100 and 200 and generated 10,000 replicates. The variance–covariance matrix of the NPMLEs was approximated by the inverse of the observed information matrix.

Table 1 presents the summary statistics of the NPMLEs of the unknown regression parameters and the baseline odds functions at the time points 0.3, 0.8, 2, and 4, the first three of which corresponds to the approximate first quartile, median and third quartile of the failure time, and the last one is the end of the study. In the table, Bias is the sampling average of the biases of the estimates; SEE denotes the average of the standard error estimates; SE denotes the empirical standard deviations of the parameter estimates; CP is the coverage probability of 95% confidence interval estimates; and PO denotes the average of the parameter estimates under the PO model.

From Table 1, we observed that the biases of the NPMLEs were small under all simulation settings. The standard deviation estimate exhibited the true variation. The coverage probabilities of the 95% confidence interval estimates were close to the nominal level. As expected, as sample size increases, the coverage probabilities improve. The PO model yielded biased results when the PO assumption is violated. Particularly, when the long-term and short-term effects are in opposite directions, the PO model had little power to detect covariate effects.

The next set of simulation studies were concerned with the performance of our proposed procedures for testing the covariate effects and the PO assumption. We used the same simulation settings as the previous simulation study except that the censoring time was from a uniform distribution in (0,20), yielding approximate censoring rates of 20% under all scenarios. We considered five hypothesis testings here: (H1) H0:β=0; (H2) H0:γ=0; (H3) H0:β=γ=0; (H4) H0:β=γ; (H5) H0:βPO=0. The first four tests were conducted based on our proposed model. The last one was conducted based on the PO model. We used Wald test statistics for testing (H1), (H2), (H4), and (H5) and likelihood ratio test statistics for testing (H3). The data were generated under the same model as above with n=400.

Table 2 presents the sizes/powers at the nominal significance level of 0.05 based on 10,000 replicates. Our method provides accurate control of type I error rates. We also observed reasonable powers of the proposed tests for testing (H3) compared to the test using the PO model when the PO assumption is true; the powers for testing the covariate effect under both scenarios β=γ=−1 and β=γ=1 are approximately 84.2% and 90.7% using the proposed model and the PO model, respectively. However, when the PO assumption is violated, the proposed tests tend to be more powerful than the test based on the PO model, especially when the long-term and short-term effects are in opposite directions. For example, when (β,γ)=(−1,1), the powers for testing (H1), (H2), and (H3) were 0.4553, 0.5213, and 0.5804, respectively, compared to the power of 0.0569 from the PO model. The proposed test for the PO assumption also yielded reasonable powers; for example, we obtained powers of 0.6755 and 0.6713 for (β,γ)=(−1,1) and (β,γ)=(1,−1), respectively.

In the last set of simulation studies, we compared the performance of our method with the method by Scheike [17]. We considered the simulation settings similar to those in Scheike [17]. Specifically, we generated data from the following model

where Z is a standard normal random variable, X is a log-normal random variable truncated at 30 with mean 0.3Z and standard deviation 0.4, A0(t)=0.1t, and A1(t)=0.05t2. We varied the value of β from –0.2 to 0.2. The failure times were censored at 12 yielding approximate censoring rates of 30%. For each setting, we considered the sample sizes of 200 and 400. We fit the proposed model with the constraint of equal short-term and long-term effects of Z and compared its performance with that of Scheike [17].

Table 3 summarizes the results for estimating the effect of Z. The NPMLE of the effect of Z under the proposed model appears to have little bias suggesting that the proposed method is robust to moderate model mis-specifications. Compared to the method of Scheike [17], the NPMLE has slightly larger biases but the empirical standard deviations are smaller. Additionally, the NPMLE appears to be more efficient in terms of the mean squared error since the method of Scheike [17] involves the estimation of two unknown functions A0(t) and A1(t) instead of one unknown function H(t) in the proposed model.

5 Example

In this Section, we apply the proposed methods to the motivating CNA study [15] described in Section 1. In this study, the copy number data for 313 lymph node-negative breast cancer patients who were selected from the tumor bank at the Erasmus Medical Center are available. Survival analysis on time to develop distant metastasis (DM) was conducted at 56859 SNP loci in the breast cancer genome. The follow-up time ranged from 20 to 169 months with a median of 99 months. Among the 313 breast cancer patients, 114 patients developed distant metastasis yielding a censoring rate of 63.6%.

There are two types of estrogen receptor (ER) for tumors, ER positive or ER negative. Among 313 patients, 199 tumors were ER positive; and 161 patients experienced menopause. The ages (years) of patients at the time of surgery were recorded ranging from 26 to 83 with a median of 54. As suggested by Zhang et al. [15], we conduct the analysis for ER-positive group and ER-negative group separately. For each group, we include surgery age and menopause status as confounders and evaluate the long-term and short-term effects of the CNV data at each of the 56859 SNP loci. As a common practice, we standardize all the covariates in the numerical implementation.

Using the proposed method, the PO assumption holds at about 99.5% of SNP loci for age and about 66.4% for menopause status. Table 4 presents the results for the ER-positive group. The first part of Table 4 presents the top 10 SNPs with the smallest P-values for testing the PO assumption. It evidently shows that the PO model is inappropriate when the long-term and short-term effects are in opposite directions. At these loci, the proposed methods are able to detect significant short-term and/or long-term genetic effects missed by the PO model. The second part of Table 4 presents the top 10 SNPs with the most significant copy number effect under the PO model. The results demonstrate that our proposed model can also detect significant long-term and/or short-term effects at those SNP loci. In contrast to the first part, we present the first 10 SNPs with the largest P-values for testing the PO assumption in the third part of Table 4. The proposed methods yield short-term and long-term effects similar to the estimates obtained from the PO model, suggesting that the proposed methods work well when the PO assumption is true.

As described in Section 1, the PO assumption appears to be true at SNP #2124, but is violated at SNP #39687. In Table 4, at SNP #2124, the estimates for the long-term and short-term effects were both −0.203 with standard error estimates of 0.285 and 0.242, respectively. The estimate for the effect under the PO model was –0.193, which is close to the estimates under the proposed model. At SNP #39687, we obtained significant long-term and short-term effects in opposite directions with estimates 1.482 and −1.130. However, the PO model yielded insignificant results with regression parameter estimate of 0.006 (P-value =0.967). These results suggest that in this case the proposed model works as well as the PO model when the PO assumption is true and outperforms the PO model when the assumption is violated.

We also calculated the adjusted odds ratio comparing patients with copy number value of 1 against 0 while fixing the two confounders at their sample means. Figure 2 displays the estimated curves of the log odds ratios with their pointwise 95% confidence bands at SNP #2124 and SNP #39687. It appears that the log odds ratio curve at SNP #2124 is close to a horizontal line, whereas the log odds ratio curve at SNP #39687 is non-decreasing and crosses 0 around 60 months suggesting violation of the PO assumption.

Figure 2

Log odds ratios (solid curves) and the 95% pointwise confidence bands (dashed curves) at SNP #2124 and SNP #39687

6 Discussion

We have proposed a novel odds rate model incorporating short-term and long-term covariate effects and have established efficient likelihood-based nonparametric estimation and inference procedures. The proposed model allows one to explicitly model the long-term and short-term covariate effects on the odds of survival without the use of smoothing. The simulation studies demonstrate that our model can detect the short-term and long-term effects when the PO assumption is violated, while the PO model is inappropriate in this circumstance. We have implemented our new method in an efficient computer program. This program is freely available upon request.

One limitation of the proposed method is that we assume inexplicitly that the odds ratio function is monotonic over time. Consider two univariate covariates X>X∗. It can be shown that

When β<γ, the odds ratio in eq. (3) is decreasing, whereas when β>γ, the odds ratio is increasing. On the other hand, the transformation model proposed by Scheike [17] postulates that

where Z is another set of covariates which satisfy the PO assumption, A(⋅) is a vector of functions, and H(t|X)≡XTA(t) is assumed to be increasing for all X. This model is more general than the proposed model in the sense that it does not impose any assumption on the odds ratio functions A(t). However, it is not clear whether the estimators of A(t) are constrained to be monotonic (see equation (7) in Scheike [17]). Alternatively, one may consider another model described by Scheike [17]

This model is not identifiable as t→0. Furthermore, it is not clear how one can ensure that the estimator of S(t|X,Z) is non-increasing.

We assume in this article that the failure time and the censoring time are conditionally independent given the covariates. Violation of this assumption may lead to biased results. For example, in the analysis of time to develop distant metastasis in the CNA study, censoring is likely to be informative because some local recurrences might censor the observations while being driven by the same tumor. To account for such dependent censoring, one may assume certain dependence structures between the failure time and the censoring time (e.g., see Tanaka and Rao [26] and Li et al. [27]). Alternatively, we can adopt the technique proposed by Othus et al. [28] through the use of the conditional crude hazard function of the censoring time given the covariates. Future research along this direction is warranted.

In this article, we assume independent observations. In several applications, we may encounter clustered or family data. One interesting future work is to extend our proposed method by incorporating random effects to account for the correlations within a cluster/family. Suppose Xij is the set of covariates for the jth subject in the ith family. We can incorporate random effects in the following model

where ωi denotes the shared frailty for the ith family. We are currently investigating the theoretical and computational aspects of this model.

Yang and Prentice [29] proposed a two-sample model that accommodates the short-term and long-term covariate effects on the hazard function. More recently, Diao et al. [30] extended the above two-sample model to a regression setting accommodating possibly external time-dependent covariates. While the models of Yang and Prentice [29] and Diao et al. [30] are extensions of the Cox proportional hazards model, the proposed model renders direct interpretations of short-term and long-term covariate effects on the survival function and includes the PO model as a special case. It is well known that both the proportional hazards model and the PO model belong to a general class of semiparamtric transformation models [22]. One may be interested in investigating the short-term and long-term covariate effects under this class of transformation models. Future research is warranted.

The flat lines after about 80 months in Figure 1 suggest that a proportion of patients in the CNA study may never develop distant metastasis, i.e., there exists a cured sub-population. To account for the presence of a cured sub-population, we can combine the proposed model and the standard mixture cure model. Specifically, given the covariates Zi and Xi, the mixture cure model postulates that the population survival function takes the form

where

and S(t|Xi) is the conditional survival function in the un-cured sub-population. Note that 1−θ(Zi) is the cure fraction. We assume the proposed model for S(t|Xi). The likelihood function is then given by

where f(t|Xi) is the conditional density function of the failure time in the un-cured sub-population. We can perform nonparametric likelihood estimation and inference using techniques similar to those in Zeng et al. [31] and Diao and Yin [32].

In this study, we performed hypothesis testing at each SNP locus to search for the significant location at which the CNV impacts the time to develop distant metastasis, leading to the so-called large p small n problem. Appropriate procedures are needed to control for the family-wise error rate introduced by the multiple comparisons. Many statistical methods have been proposed in literature to account for multiple testing. Particularly, Lin [33] developed a Monte Carlo approach to approximate the joint distribution of the test statistics along the genome. It would be interesting to adopt this approach in our model setting to adjust for multiple testing.

An interesting area is focused on the selection of important variables that impact the failure time. Recently, Lu and Zhang [34] proposed variable selection procedures using adaptive LASSO for the PO model. It would be desirable to extend the adaptive LASSO procedure to select variables with significant long-term and/or short-term effects.

7 Software

Program in the form of C code is available on request from the corresponding author.