On a simple estimation of the proportional odds model under right truncation

Retrospective sampling can be useful in epidemiological research for its convenience to explore an etiological association. One particular retrospective sampling is that disease outcomes of the time-to-event type are collected subject to right truncation, along with other covariates of interest. For regression analysis of the right-truncated time-to-event data, the so-called proportional reverse-time hazards model has been proposed, but the interpretation of its regression parameters tends to be cumbersome, which has greatly hampered its application in practice. In this paper, we instead consider the proportional odds model, an appealing alternative to the popular proportional hazards model. Under the proportional odds model, there is an embedded relationship between the reverse-time hazard function and the usual hazard function. Building on this relationship, we provide a simple procedure to estimate the regression parameters in the proportional odds model for the right truncated data. Weighted estimations are also studied.


Introduction
Truncation is common in survival analysis where the incomplete nature of the observations is due to a systematic biased selection process originated in the study design. Right truncated data arise naturally when an incubation period (i.e., the time between disease incidence and the onset of clinical symptoms) cannot be observed completely B Peng Liu p.liu@kent.ac.uk in a retrospective study. In survival analysis, right truncation will lead to biased sampling in which shorter observations will be oversampled (Gürler 1996). For example, to study AIDS caused by blood transfusion (Lagakos et al. 1988), the incubation period is the time from a contaminated blood transfusion to the time when symptoms and signs of AIDS are first apparent. However, in those studies, the following-up period are usually limited. Therefore, only those developed AIDS before the end of study can be identified.
Many authors have studied right truncated data: Woodroofe (1985) and Wang et al. (1986) focused on the asymptotic properties the product limit estimator under random truncation. Keiding and Gill (1990) studied asymptotic properties of random left truncation estimator by a reparametrization of the left truncation model as a threestate Markov process. Lagakos et al. (1988) considered nonparametric estimation and inference of right truncated data by treating the process in reverse time, they showed that λ B (t) = λ(τ − t), where τ is the study duration, λ B (t) and λ(t) are reverse-time hazard and forward-time hazard, respectively. The authors also discussed the implications and limitations of introducing reverse time hazard to analyze right truncated data. Gross (1992) further explained the necessity of reverse time hazard in the Cox model setting.
However, in most of the current literature, researchers study right truncated data in nonparametric setting, fairly few studied semiparametric models, among them, Kalbfleisch and Lawless (1989) formulating the Cox model on the reverse time hazard (or retro hazard, Lagakos et al. (1988); Keiding and Gill (1990)). For other related work on reverse time hazard, please refer to Gross (1992); Chen et al. (2004), among others.
In this paper, we study right truncated data under a semiparametric proportional odds model. Different from a proportional hazards model, the reverse-time hazard in proportional odds model has a simple log-linear relationship with the forwardtime hazard, which leads to an intuitive estimator. While Sundaram (2009)'s method can also be adapted to proportional odds model for right truncated data, she focused on applying a reversed-time argument to an estimator for left truncated data. Our estimator, on the other hand, utilize a direct relationship between the reverse-time hazard, the forward-time hazard and the baseline odds function, so that we obtain a simpler estimator. Weighted functions are also being inserted into the estimating equation to obtain more efficient estimates.
The rest of the paper is organized as follows. Section 2 describes the inference procedure as well as asymptotic results, Sect. 3 shows simulation and real data results, Sect. 5 provides some discussion. Proof of theorems are left into the Appendix part.

Inference procedure
Assume that the failure time of interest T follows the semiparametric proportional odds model: and the observed failure time is subject to a right truncation time variable R. The observed data is (T i , R i ), i = 1, . . . , n, where T i ≤ R i . Let τ be the study duration, which is greater than max{T 1 , T 2 , . . . , T n }. An (observed) reverse-time sample, (T * i , R * i ), i = 1, . . . , n can be constructed, where T * = τ − T , R * = τ − R, so that T * is left truncated by the variable R * . Denote (T * ,R * ) as the reverse-time sample (potentially truncated). Then the hazard function ofT * is a quantity originated in τ and counts backward in time. The reverse hazard and cumulative reverse hazard function of backward recurrence time is defined as We would like to mention that a similar definition of the reverse hazard can also be found in Kalbfleisch and Lawless (1989) and Jiang (2011). Denote v(t) = exp(α(t)), and λ(t) = f (t)/S(t) as the forward-time hazard, then dt .

Consider the counting process
and denote Then M i (t, β) is a martingale with respect to the self-exciting (canonical) filtration (Keiding and Gill 1990;Stralkowska-Kominiak and Stute 2009) and Multiply both sides of (2) by {exp(Z i β)v(t)+1} and summing over n observations, Divide both left-hand side and right-hand side by n i=1 Y i (t), we obtain: .
which is equivalent to: Denote the left-hand side of (4) as: .
From standard counting process arguments (Anderson and Gill, 1982;Aalen10), we know that the stochastic integral with respect to the counting process martingale M i (dt, β) is also a martingale, motivate by the following equation We construct the following estimating equation Only v(t) is unknown in (5), let the estimate of v(t) bev n (t, β). Denote Multiply (2) by Z i {exp(Z i β)v(t) + 1}/n and summing over n observations, we obtain By virtue of the same idea of (5), take integration on both sides of (7), we can also construct another equation: Substituting (6) into (8), we can obtain the estimate of β by solving the following equation: .
Finally, let and denote the solution of S n (β) = 0 beβ n , we have the following theorem: Remark For proportional odds model with the normal logit link: DefineM as a result, we can obtain .
Following the definition of reverse hazard in Sect. 2, we can write the reverse hazard asλ From the general definition of martingale in Fleming and Harrington (1991) (pp. 25), we can easily show thatM i (t, β) is a martingale. While for model (1), , The corresponding estimating equation under model (10) has the following form .
Equation (11) also can be used to estimate β, however, comparing with (9), (11) is more complicated and more computational intensive, while the derivative of (9) with respect to β can be easily obtained. As a result, (9) can be easily solved by the newton raphson algorithm. In the following simulations, we will use estimating equation (9). In addition to the unweighted object function (9), weighted object function can also being included to obtain a class of weighted estimators of β 0 . This procedure is often used to minimize the sandwich estimate as well as improve the efficiency. The weighted version of object function is here W n (t) is a predictable weight function with respect to the canonical filtration which converges to a non-random function w(t). One of the common used weight function is the Prentice-Wilcoxon type function W n1 (t) =Ŝ L B (t), whereŜ L B (·) is the Lynden Bell estimate of the baseline survival function for right truncated failure time data. Denote the corresponding estimate of β asβ n,w . Then we have the following theorem: Theorem 2 Under the same assumptions as Theorem 1, when n → ∞, for a prespecified weight function W n (·) → w(·), √ n(β n,w − β 0 ) converges weakly to a mean-zero normal distribution, with covariance matrix Recently, many people considered problem of finding the optimal weight in a weighted estimating equation, including Chen and Cheng (2005); Chen and Wang (2000); Chen et al. (2012), among others. To achieve this goal, we only need to find the w(t) such that U w (β 0 ) −1 V w (β 0 )U w (β 0 ) −1 achieves the minimum. Since both the empirical weight function W n (t) and its limit w(t) do not rely on unknown parameter β 0 , it is reasonable to set β 0 = 0. Another explanation for letting β 0 = 0 is that it represents the baseline distribution. Therefore, let β 0 = 0, then we have: Apply the Cauchy-Schwarz inequality to U w (β 0 ) −1 V w (β 0 )U w (β 0 ) −1 and let β 0 = 0, then it follows that the optimal weight is proportional to which minimize the variance ofβ n . Since when (15) holds, we have In simulation, let W n2 (t) =Ŝ L B (t) 1 −Ŝ L B (t) , the results are shown in Table 1 , it can be seen that the weight W n2 (t) achieve the minimal variance among the three estimators.

Simulation and real data
We perform simulation studies to evaluate the finite sample properties of the proposed estimator. In simulation, let α(t) = 3 log t, β 0 = (1, 0.5) , Z 1 is a continuous variable follows a uniform distribution from 0 to 2, Z 2 is a discrete variable follows a Bernoulli distribution with probability 0.5. The failure time variable is generated from model (1). The right truncation variable follows a uniform distribution from 0 to 4. This makes the truncation rate equals to 20%. For each simulation, 1000 datasets are generated, in each dataset, there are n observations, n = 300, 400, 500, 600, respectively. W n1 (t) and W n2 (t) are chosen as the weight functions in weighted estimating equations. As is shown in Table 1, three estimation equations yield unbiased estimates and the empirical coverage probability is around nominal level 95%, when weighted function is incorporated into the estimation equation, the efficiency is greatly improved, and the variance achieve minimal for W n2 (t) under three estimates. As pointed out by one of the referees and the associate editor, Shen et al. (2017) also studied right truncated data under linear transformation models, and we know that when the error term in the linear transformation model follows logistic distribution (Fine et al., 1998), the model becomes the proportional odds model. Let then the estimating equations (3) and (4) in Shen et al. (2017) can be written as We recognize that Shen et al. (2017)'s methodology is general and works for all the linear transformation models, including the proportional odds model. However, our approach will be more convenient compared with Shen et al. (2017)'s under the proportional odds model, since our approach has a simpler form, and the estimation of the intercept α(t) can be done beforehand and plugged in the final estimating equation, while Shen et al. (2017) can not achieve this and their estimation produce involves a complicated iteration which increases the risk of non-convergence. Besides, Shen et al. (2017) only deal with the reverse time but not the reverse hazard function, and we utilize the relationship between the reverse hazard function and the forward-time hazard function and produced a more intuitive estimator.
We conduct simulations for Shen et al. (2017)'s method and the results are reported in Table 1. The code was obtained from the authors via personal communication. However, one of the authors, Prof. Pao-Sheng Shen mentioned that they were unable to calculate the asymptotic variance and coverage probabilities, the existing results in their paper contain some errors, and their current code only consists of bias and standard error. As a result, we only report bias and standard error of Shen et al. (2017)'s method. All the simulations were conducted under the same model as ours. We also want to mention that we found the computation speed is very slow for Shen et al. (2017)'s method, though asymptotic variance and coverage probability were not calculated, their method is still more than 3 times slower than ours under the same model setting and the sample size. The SSE of Shen et al. (2017)'s method is smaller than our unweighted estimator, but is bigger than the two weighted estimators. For the second approach in their paper, i.e. the conditional maximum-likelihood approach, since the bias is large, we did not perform further comparisons here. We would like to mention that the large bias of the conditional maximum-likelihood approach is also confirmed in .
As suggested by one of the reviewers, we also perform simulations without accounting for the truncation, and the results are shown in Table 2. We choose the truncation distribution as uniform distributions from 0 to 4, 2 and 1, respectively, which corresponds to 20% truncation rate (mild truncation), 40% truncation rate (moderate truncation) as well as 70% truncation rate (heavy truncation). As we can see from Table 2, all the estimators are biased, and a larger truncation rate will lead to a bigger bias and variance, though for the same truncation, variances will decrease when the sample sizes increase. These results also coincide with Table 2.1 (pp. 20) in Rennert (2018) and Table 1 in Rennert and Xie (2018), though the two articles deal with the doubly truncated data under the Cox model.
To better illustrate how to employ the proposed method in real situation, we analyze the Centers for Disease Control's blood-transfusion data, this data was used by Kalbfleisch and Lawless (1989) and Wang (1989). The data include 494 cases reported to the Center of Disease Control prior January, 1, 1987, and diagnosed before July, 1, 1986. Only 295 of the 494 has consistent data, and they got infection by a single blood transfusion or a short series of transfusions, analyse is restricted to this subset. We obtain the raw observation data via personal communication, Thomas Peterman, Centers for Disease Control and Prevention. The data contains three variables: T is the time from blood transfusion to the diagnosis of AIDS (in months), R is the time from blood transfusion to the end of the study (July, 1986, in months), Age is the age of the person when transfusing blood (in years). Comparing the data with Kalbfleisch and Lawless (1989)'s as well as Wang (1989)'s, the observation (X =16, T =33, Age=34) cannot be found in the raw data, thus is being deleted and the final sample size is 294, and a few fractions of the data are also corrected because these entries are not correct compared to the raw data.
We apply the proposed method to this data and treat Age as the covariate in regression. In Wang (1989)'s paper, the data are categorized into three age groups: 'children' aged 1-4, 'adults' aged 5-59, and 'elderly patients' aged 60 and older because of different patterns of survivorship, the survivor behaviour of groups 'adults' and 'elderly patients' are similar except for the right tail while there is an evident distinction compared with 'children', in current analysis, we delete the data from 'children', and focus on a combined sample of 'adults' and 'elderly patients' with a sample size equal to 260. Finally, the range of T is from 0 to 89, and the range of R is from 0 to 99. For all i ∈ {1, . . . , 260}, we have T i ≤ R i . As a result, our dataset will not have the Table 2 Simulation results when ignoring truncation  identifiability issue as mentioned in Seaman et al. (2022). We also applied Shen et al. (2017)'s method and the result is similar. All the results are shown in Table 3 , where the weights are chosen as W n1 (t) and W n2 (t), the estimated parameter between unweighted and weighted estimation equation does not show much difference, but the variance is reduced when weights are considered. In both situations mentioned above, Age has a very weakening positive effect on the odds ratio, but the effect is not significant.

Discussion
Directly consider the right truncated data in normal time order can be failed because 'at risk' process is not adapt to the history of the process (Gross 1992). Retro hazard solves this problem which transform right truncated data to left truncated in reverse time (Woodroofe 1985). Statistical modelling is even more flexible by incorporating the nature structure of proportional odds model. The usual form of proportional odds model can also be utilized but the theoretical and computational burden for the estimator will be increased, employ (1) can substantially improve the situation.

A4: f (t) is continuous.
Assumption A1 is also used by Chen et al. (2012), A2 is a standard assumption to ensure martingale properties holds (Fleming and Harrington 1991), A3 is also a standard assumption to avoid theoretical discussion, it is also being used in Huang and Qin (2013), A4 is being used in prove the martingale representation ofv n (t, β 0 )−v(t).
Besides that, we also need an condition to ensure that the truncated distribution to be correctly identified, let F(·) and G(t) be the distribution function of T and R, define (a F , b F ) and (a G , b G ) be the support of F(·) and G(·) of T and R under the meaning that a W = inf{x : W (x) > 0}, b W = sup{x : W (x) < 1}, where W is a distribution function. Under right truncation, actually, only conditional distribution P(T ≤ x|T ≤ b G ) and P(R ≤ x|R ≥ a F ) can be estimated, thus we assume a F = a G = 0, b R = ∞, so that the conditional distribution will be the actual distribution of T and R, we also assume P(T ≤ Y ) = α > 0 to ensure that there exist observations satisfy our condition, similar assumption and discussion also appeared in Woodroofe (1985); Wang (1989), and Sundaram (2009), among others.

Proof of Theorem 2
Since Theorem 1 and 2 are quite similar, in this part, we will omit the proof detail and only give the detailed expression of V w . where