NEJSDS The New England Journal of Statistics in Data Science 2693-71662693-7166 New England Statistical Society NEJSDS108 10.51387/26-NEJSDS108 Methodology Article Statistical Methodology Estimating Nonlinear Exposure and Event Time Association in the Presence of Exposure Measurement Error ThomsonTrevor J.tthomson@fredhutch.org∗ HuangYingyhuang@fredhutch.org Fred Hutchinson Cancer Center, Vaccine and Infectious Disease Division, Seattle, WA, USA. E-mail address: tthomson@fredhutch.org Fred Hutchinson Cancer Center, Vaccine and Infectious Disease Division, Seattle, WA, USA. E-mail address: yhuang@fredhutch.org Corresponding author. 20261510202600114972026 © 2026 New England Statistical Society2026 Open access article under the CC BY license.

Estimating nonlinear exposure effects in Cox regression models is essential for accurately characterizing exposure–disease relationships, yet the presence of measurement error in exposures poses major challenges. While regression calibration has been widely applied under rare disease assumptions with linear exposure effects, recent work has extended the approach to accommodate quadratic and quantile effects. We generalize this framework to handle arbitrary polynomial and spline-based representations of the exposure effect, thereby enabling flexible modelling of complex nonlinear associations. Our method combines regression calibration with Monte Carlo integration to evaluate the conditional expectation of the hazard function when analytical solutions are unavailable. Furthermore, by stacking the estimating functions from both the Cox regression and measurement error models, we construct a sandwich variance estimator that avoids computationally intensive bootstrap procedures. To facilitate practical use, we implement our approach in the R package NonLinearRC.

 Cox regression Measurement error Monte Carlo approximation Nonlinear effect Regression calibration Sandwich variance estimator This research is supported by the National Institute of Health grants R01 CA119171, R01CA277133, and R37AI054165. Introduction

Estimating the effect of a continuous exposure on disease risk is central to understanding disease etiology and guiding prevention strategies. When the outcome of interest is time to disease onset, Cox’s proportional hazards model is a widely used tool for characterizing the relationship between exposure and risk. However, the association between the exposure and the outcome may be complex and nonlinear, making a single regression coefficient insufficient to fully capture the nature of the relationship. To address this, researchers have incorporated nonlinear effects into the Cox model by including higher-order polynomial terms for the exposure, or by representing the exposure using linear combinations of either spline basis functions or indicator variables corresponding to exposure percentiles [15, 5]. Conventional partial likelihood methods can be applied to fit Cox regression models when a random sample of individuals is available with observations on the time-to-event outcome (subject to right censoring) and the exposure of interest. In many practical settings, however, the true exposure is unobserved, and only an error-prone surrogate is available. An example arises in nutritional epidemiology, where traditional dietary exposure measures based on self-reported intake are often subject to both random error and systematic bias [5]. It is well established that naively ignoring measurement error can lead to attenuation of effect estimates and, consequently, bias the assessment of the true exposure–disease relationship [2, 21].

In settings where disease onset is rare, regression calibration is a widely used method for correcting covariate measurement error in Cox regression models [12]. The method replaces the original Cox model with its conditional expectation given the observed covariate, and the parameters of interest are parameterized by the induced hazard function. In nutritional epidemiology, objective biomarkers, often conceptualized as true dietary intake plus random measurement error, have commonly been used to calibrate self-reported intake for subsequent diet–disease association analyses [13, 6, 7]. Under a Cox model with only a linear term for the exposure, the induced hazard function has an analytical form, which makes its implementation compatible with standard statistical software. Recently, [5] developed a regression calibration approach for Cox models that incorporate either exposure quantiles or a quadratic exposure term and applied it to the Women’s Health Initiative Dietary Modification Trial [18]. Using biomarker data to calibrate self-reported sodium and potassium intake, they examined the associations of these exposures with cardiovascular disease (CVD) risk to assess potential departures from linearity. Their work provides a natural epidemiologic setting in which more flexible nonlinear modelling is scientifically relevant. The goal of the present paper is to extend that framework beyond prespecified nonlinear terms to a broader class of exposure–response models, including arbitrary polynomial and spline-based representations. Accordingly, we focus on methodological development and controlled numerical studies, which allow the operating characteristics of the proposed method to be evaluated systematically across a wider range of nonlinear settings. The added flexibility of our proposed modelling enables researchers to capture complex nonlinear associations between the exposure and event time, while maintaining compatibility with standard tools for model assessment and variable selection. Moreover, by stacking the estimating functions for both the Cox regression parameters and the parameters of the measurement error model, we derive a variance estimator for the regression parameters following the approach of [9]. This eliminates the need for computationally intensive bootstrap procedures, such as those used in [5]. The main challenge lies in evaluating the conditional expectation of the hazard function, which defines the induced hazard. While [5] obtained a closed-form expression under a quadratic effect and approximated the expression for quantile effects with a Taylor series approximation, these strategies do not generalize to broader settings. However, because the conditional distribution of the exposure given the observed covariates is specified, the expectation can naturally be approximated using Monte Carlo integration. Since standard survival analysis software cannot accommodate this generalized setting, we have developed an R package, NonLinearRC, to facilitate implementation and overcome this barrier.

The remainder of the article is organized as follows. Section 2 introduces the notation and modelling framework. Section 3 outlines our proposed method, which generalizes regression calibration to accommodate nonlinear exposure effects, and describes the corresponding variance estimator. In Section 4, we evaluate the finite-sample performance of the proposed procedure through a simulation study. In Section 5, we provide a tutorial on utilizing our R package NonLinearRC to apply our proposed estimation procedure outlined in Section 3. Concluding remarks are provided in Section 6.

 Methodological Framework Notation

Let Z ∈ R denote the true exposure variable of interest. Our primary objective is to evaluate the association between Z and the time to disease onset, denoted by T. In our random sample, referred to as the analysis cohort, the primary statistical challenges arise from right-censoring of T and the presence of measurement error in Z. Specifically, we observe an error-prone surrogate Q in place of Z, subject to both systematic and random error. In addition, we collect baseline covariates V = ( X ′ , R ′ ) ′ , where X includes potential confounders of the exposure–outcome relationship, and R comprises variables that influence the measurement error mechanism, i.e., the relationship between Q and Z. To correct for measurement error in Z, we incorporate a calibration dataset comprising participants drawn either from the analysis cohort or from an external study. For each participant in the calibration sample, we observe Q, V, and an objective biomarker W, which represents the true value of Z plus random noise. Furthermore, replicate measurements of W are available for a subset of participants in the calibration dataset and are used to estimate the measurement error variance in W.

Overall, the observed data can be partitioned into three components: D 1 = { ( U i , δ i , Q i , V i ) : i ∈ I 1 } , D 2 = { ( W i 1 , Q i , V i ) : i ∈ I 2 } , D 3 = { ( W i 1 , … , W i , m i ) : i ∈ I 3 } , where U i = min { T i , C i }, δ i = I ( T i ≤ C i ), C i denotes a non-informative censoring time for the event of interest. The sets I j ( j = 1 , 2 , 3) index participants in the analysis cohort, the calibration dataset, and the replicate dataset respectively, and m i denotes the number of replicate measurements available for participant i ∈ I 3 . We use D 1 to estimate the association between T and Z, and employ D 2 ∪ D 3 to correct for the measurement error in Q.

 Measurement Error Modelling

Adapted from [10], we assume that the observed exposure Q follows a measurement error model subject to systematic bias influenced by subject-specific characteristics V: Q i = a 0 + a 1 Z i + a 2 ′ V i + ε i , where ε i ∼ iid N ( 0 , σ ε 2 ) and ε i ⊥⊥ { Z i , V i } for all i. We additionally suppose that the biomarker W adheres to a classical measurement error model: W i j = Z i + e i j , where e i j ∼ iid N ( 0 , σ e 2 ), and e i j ⊥⊥ { Z i , V i , ε i } for all j.

 Hazard Modelling

To investigate the association between T and Z, we specify the following hazard regression model: λ ( t ; Z i , Q i , V i ) = lim Δ t → 0 + 1 Δ t P ( t ≤ T i < t + Δ t | T i ≥ t , Z i , Q i , V i ) = λ ( t ; Z i , Q i , X i ) = λ ( t ; Z i , X i ) = λ 0 ( t ) exp { g ( Z i ; β ) + γ X ′ X i } , where (2.3) follows upon assuming the measurement error in (2.1) is nondifferential (i.e., T i ⊥⊥ Q i | Z i ) [2, 21]. Here, λ 0 ( t ) denotes an unspecified baseline hazard function, and g ( Z ; β ) is a known function of Z parameterized by β. This formulation allows for a potentially nonlinear effect of Z on the hazard for T, with β characterizing the strength and form of this association. Common choices for g ( Z ; β ) include a degree-p polynomial, or a linear combination of spline basis functions with percentile-based knots [3].

Our objective can thus be formulated as the statistical problem of estimating the regression parameters β and γ X , along with the baseline hazard function λ 0 ( · ), using the available data D 1 ∪ D 2 ∪ D 3 .

 Nonlinear Regression Calibration Nuisance Parameter Estimation

Suppose we are willing to specify Z i | V i ∼ iid N ( μ ( V i ) , S 2 ), where μ ( V i ) is a known function of V i and S 2 is not a function of V; our modelling can be extended to relax this latter assumption. It then follows from (2.1) that Z i | { Q i , V i } ∼ iid N ( α 0 + α 1 Q i + α 2 ′ V i , σ 2 ), where α 1 = a 1 S 2 a 1 2 S 2 + σ ε 2 , α 0 = − α 1 a 0 , α 2 ′ V = μ ( V ) − α 1 a 1 μ ( V ) − α 1 a 2 ′ V , σ 2 = S 2 σ ε 2 a 1 2 S 2 + σ ε 2 . Since α 0 , α 1 , and α 2 are unknown, we proceed to estimate them with data D 2 ∪ D 3 . Given that e i j in (2.2) is assumed to be normally distributed and independent of { Z i , V i , ε i j }, it follows that W i j | { Q i , V i } ∼ iid N ( α 0 + α 1 Q i + α 2 ′ V i , σ 2 + σ e 2 ) . This formulation suggests that the parameters α 0 , α 1 , α 2 , and ( σ 2 + σ e 2 ) can be estimated by fitting a linear regression model with data D 2 , in which W i 1 is specified as the response variable, and Q i and V i are predictors. Let α ˆ 0 , α ˆ 1 , and α ˆ 2 denote the estimated regression coefficients, and let σ 2 + σ e 2 ˆ represent the estimated residual variance from this model.

To separately estimate σ e 2 , we utilize replicate biomarker measurements in D 3 to obtain the following estimator σ ˆ e 2 = 1 n 3 ∑ i ∈ I 3 1 m i − 1 ∑ j = 1 m i ( W i j − W ¯ i ) 2 , where W ¯ i = ∑ j = 1 m i W i j / m i is the within-subject average, and n 3 is the number of participants in I 3 . An estimator for σ 2 is then obtained as σ ˆ 2 = σ 2 + σ e 2 ˆ − σ ˆ e 2 .

 Hazard Regression Estimation

Since Z is unobserved, we cannot directly evaluate the hazard function in (2.4). Instead, we consider the following induced hazard function based on the observed data in D 1 : λ ( t ; Q i , V i ) = ∫ λ ( t ; z , X i ) d F ( z | T i ≥ t , Q i , V i ) , where λ ( t ; Z i , V i ) is shown in (2.4). Evaluation of (3.2) requires specifying Z i | { T i ≥ t , Q i , V i }, which is generally unknown and difficult to model. However, under the rare disease assumption, where the probability of survival up to time t satisfies P ( T i ≥ t | Z i , V i ) ≈ 1 for time points of practical interest, we can approximate the conditional distribution using Bayes’ Theorem: d F ( z | T i ≥ t , Q i , V i ) ≈ d F ( z | Q i , V i ). Substituting this into (3.2), the induced hazard can be approximated as λ ( t ; Q i , V i ) ≈ λ 0 ( t ) exp { γ X ′ X i } H ( β , Q i , V i ; θ ) , where H ( β , Q i , V i ; θ ) = E Z ( exp { g ( Z ; β ) } | Q i , V i ), with the expectation taken with respect to the conditional distribution F Z ( z | Q i , V i ; θ ), and θ denotes the parameters governing this distribution. Under our modelling and parametric assumptions thus far, we have that F Z ( z | Q i , V i ; θ ) is Gaussian with expectation shown in (3.1) and variance σ 2 .

For fixed θ, we propose to estimate β and γ X by solving the following system of estimating equations: Φ ( β , γ X ; θ ) = 1 n 1 ∑ i ∈ I 1 Φ i β ( β , γ X ; θ ) Φ i γ ( β , γ X ; θ ) = 0 , where Φ i β ( β , γ X ; θ ) = δ i ∂ H ( β , Q i , V i ; θ ) / ∂ β H ( β , Q i , V i ; θ ) − S β ( 1 ) ( U i , β , γ X ; θ ) S ( 0 ) ( U i , β , γ X ; θ ) Φ i γ ( β , γ X ; θ ) = δ i X i − S γ ( 1 ) ( U i , β , γ X ; θ ) S ( 0 ) ( U i , β , γ X ; θ ) , with S ( 0 ) ( t ; β , γ X , θ ) = ∑ i ∈ I 1 Y i ( t ) exp { γ X ′ X i } H ( β , Q i , V i ; θ ) , S β ( 1 ) ( t ; β , γ X , θ ) = ∑ i ∈ I 1 Y i ( t ) exp { γ X ′ X i } ∂ H ( β , Q i , V i ; θ ) ∂ β , S γ ( 1 ) ( t ; β , γ X , θ ) = ∑ i ∈ I 1 Y i ( t ) exp { γ X ′ X i } H ( β , Q i , V i ; θ ) X i , Y i ( t ) = I ( U i ≥ t ), and n 1 denotes the number of participants in I 1 .

Since θ is unknown in practice, we replace it with an estimate, θ ˆ (see Section 3.1), and obtain β ˆ and γ ˆ X as the solution to Φ ( β , γ X ; θ ˆ ) = 0. The approach to evaluating H ( β , Q , V ; θ ) depends on the specific form of g ( Z ; β ). For example, if g ( Z ; β ) = β Z, then H ( β , Q , V ; θ ) = exp { β ( α 0 + α 1 Q + α 2 ′ V ) + β 2 σ 2 / 2 }, which follows from the moment-generating function of the normal distribution. In the case where g ( Z ; β ) = β 1 Z + β 2 Z 2 , [5] show that H ( β , Q , V ; θ ) is 1 1 − 2 β 2 σ 2 exp { − 1 2 ( ( α 0 + α 1 Q + α 2 ′ V ) 2 σ 2 − ( β 1 + ( α 0 + α 1 Q + α 2 ′ V ) / σ 2 ) 2 1 / σ 2 − 2 β 2 ) } . This leads to the approximation of the induced hazard in (3.3) as λ 0 ∗ ( t ; β 1 , β 2 , σ 2 ) × exp { β ˜ 1 ( α 0 + α 1 Q i + α 2 ′ V i ) + β ˜ 2 ( α 0 + α 1 Q i + α 2 ′ V i ) 2 + γ X ′ X i } , where λ 0 ∗ ( t ; β 1 , β 2 , σ 2 ) = λ 0 ( t ) 1 1 − 2 β 2 σ 2 exp { σ 2 2 β 1 2 1 − 2 β 2 σ 2 } , β ˜ k = β k 1 − 2 β 2 σ 2 k = 1 , 2 . One can then estimate β ˜ k and Λ 0 ∗ ( t ) = ∫ 0 t λ 0 ∗ ( u ) d u using the maximum partial likelihood estimator and the Breslow estimator, respectively, and subsequently estimate β k and λ 0 ( · ) by inverting the relationships shown above. In general, H ( β , Q , V ; θ ) does not have a closed-form expression when g ( Z ; β ) is a degree-p polynomial with p > 2, or a linear combination of basis functions. However, since the conditional distribution of Z | { Q , V } can be deduced and estimated using the data D 2 ∪ D 3 , we approximate H ( · ) using Monte Carlo integration: H ( β , Q i , V i ; θ ) ≈ 1 C ∑ c = 1 C exp { g ( Z i c ; β ) } , where Z i c | { Q i , V i } ∼ iid N ( α 0 + α 1 Q i + α 2 ′ V i , σ 2 ), and C is sufficiently large. Once θ is estimated by θ ˆ, we can generate samples { Z i c } c = 1 C for all i ∈ I 1 . In practice, the Monte Carlo size C should be chosen large enough that the resulting estimates are numerically stable. A practical strategy is to begin with a moderate value of C and increase it until the parameter estimates change negligibly. In our experience, C = 1000 provided a reasonable starting value in the examples considered here, although the required value may vary depending on the complexity of g ( · ) and the underlying measurement error model.

The proposed approach requires specification of the conditional distribution Z | { Q , V }. In the present paper, we assume this distribution is Gaussian because it yields tractable conditional moments for Taylor approximation and allows straightforward sampling for Monte Carlo approximation. The approximation strategies themselves are not conceptually restricted to the Gaussian setting; more generally, alternative measurement error models could be accommodated provided they permit valid estimation, approximation, or sampling from Z | { Q , V }.

As an alternative to the Monte Carlo approximation for H ( β , Q i , V i ; θ ), one can apply a second-order Taylor series expansion of exp { g ( Z ; β ) } around Z = α 0 + α 1 Q i + α 2 ′ V i to obtain H ( β , Q i , V i ; θ ) ≈ exp { g ( α 0 + α 1 Q i + α 2 ′ V i ; β ) } + σ 2 2 d 2 d Z 2 exp { g ( Z ; β ) } Z = α 0 + α 1 Q i + α 2 ′ V i .

In Section 4, we compare the resulting estimators of ( β ′ , γ X ′ ) ′ obtained using both first- and second-order Taylor series approximations to those derived via the proposed Monte Carlo approximation.

 Baseline Hazard and Survival Probability Estimation

By treating λ 0 ( · ) as a piecewise constant function between uncensored failure times, d Λ 0 ( t ) = λ 0 ( t ) d t becomes a finite-dimensional parameter. Then, for fixed Ω = ( β ′ , γ X ′ , θ ′ ) ′ , the following expression defines an unbiased estimating equation for d Λ 0 ( t ) under the model in (2.4) and the rare disease assumption: ∑ i ∈ I 1 Y i ( t ) [ d N i ( t ) − exp { γ X ′ X i } H ( β , Q i , V i ; θ ) d Λ 0 ( t ) ] = 0 , where N i ( t ) = I ( T i ≤ t ) is the counting process for participant i, and Y i ( t ) = I ( U i ≥ t ) denotes the at-risk indicator. Solving the estimating equation above for d Λ 0 ( t ) yields the following estimator: d Λ ˆ 0 ( t ; β , γ X , θ ) = ∑ i ∈ I 1 Y i ( t ) d N i ( t ) ∑ j ∈ I 1 Y j ( t ) exp { γ X ′ X j } H ( β , Q j 1 , V j ; θ ) . The corresponding estimator for the cumulative baseline hazard function, Λ 0 ( t ) = ∫ 0 t d Λ 0 ( u ), is given by Λ ˆ 0 ( t ; β ˆ , γ ˆ X θ ˆ ) = ∫ 0 t d Λ ˆ 0 ( t ; β ˆ , γ ˆ X , θ ˆ ).

We can thus estimate the survival probability at time t for a subject with exposure Z = z and covariates X = x using P ˆ ( T > t | z , x ; β ˆ , γ ˆ X , θ ˆ ) = exp { − Λ ˆ 0 ( t , β ˆ , γ ˆ X ; θ ˆ ) exp { g ( z ; β ˆ ) + γ ˆ X ′ x } } .

Asymptotics and Variance Estimation

Let Ω = ( β ′ , γ X ′ , θ ′ ) ′ denote the full parameter vector, and let Ω ˆ = ( β ˆ ′ , γ ˆ X ′ , θ ˆ ′ ) ′ denote its estimator that solves the estimating equation Γ ( Ω ) = ∑ i ∈ I 1 ∪ I 2 ∪ I 3 Φ i β ( β , γ X ; θ ) Φ i γ ( β , γ X ; θ ) Ψ i ( θ ) = 0 , where Ψ ( θ ) = ∑ i ∈ I 2 ∪ I 3 Ψ i ( θ ) is the estimating function for the nuisance parameter θ, as defined in Appendix A. We define each Φ i β ( β , γ X ; θ ) = 0 and Φ i γ ( β , γ X ; θ ) = 0 for i ∉ I 1 and all Ω, and similarly, Ψ i ( θ ) = 0 for i ∉ I 2 ∪ I 3 and all θ.

Using results for M-estimators and arguments in Section 2.3 of [1], we have that Ω ˆ → p Ω ∗ as m → ∞, where m = min { n 1 , n 2 , n 3 }, and Ω ∗ is the solution to E ( Γ ( Ω ) ) = 0. Define the empirical derivative matrix A ˆ ( Ω ) = ∑ i ∈ I 1 ∪ I 2 ∪ I 3 ∂ ∂ Ω ′ Γ i ( Ω ) → p A ( Ω ) as m → ∞, where A ( Ω ) = ∂ E ( Γ ( Ω ) ) / ∂ Ω ′ is assumed to be negative definite. Let Ω 0 denote the true parameter value. Then under correct model specification, the consistency of Ω ˆ depends on the approximation used for H ( β , Q i , V i ; θ ):

Monte Carlo Approach: If we apply a Monte Carlo approximation to H ( · ), then Ω ∗ = Ω 0 provided the rare disease assumption holds.

Second Order Taylor Series Approach: If a second-order Taylor approximation is used, then Ω ∗ = Ω 0 provided the rare disease assumption holds and the higher-order terms of the Taylor Series expansion are negligible.

 Following [9], the estimating function can be expressed as Γ i ( Ω ) = Γ i ∗ ( Ω ) + o p ( 1 ), where Γ i ∗ ( Ω ) = Φ i β ∗ ( β , γ X ; θ ) Φ i γ ∗ ( β , γ X ; θ ) Ψ i ( θ ) , Φ i β ∗ ( β , γ X ; θ ) = ∫ 0 ∞ ∂ H ( β , Q i , V i ; θ ) ∂ β − s β ( 1 ) ( t , β , γ X ; θ ) s ( 0 ) ( t , β , γ X ; θ ) d N i ( t ) − ∫ 0 ∞ Y i ( t ) exp { γ X ′ X i } H ( β , Q i , V i ; θ ) s ( 0 ) ( t , β , γ X ; θ ) ∂ H ( β , Q i , V i ; θ ) ∂ β − s β ( 1 ) ( t , β , γ X ; θ ) s ( 0 ) ( t , β , γ X ; θ ) d E ( N ¯ ( t ) ) , Φ i γ ∗ ( β , γ X ; θ ) = ∫ 0 ∞ X i − s γ ( 1 ) ( t , β , γ X ; θ ) s ( 0 ) ( t , β , γ X ; θ ) d N i ( t ) − ∫ 0 ∞ Y i ( t ) exp { γ X ′ X i } H ( β , Q i , V i ; θ ) s ( 0 ) ( t , β , γ X ; θ ) X i − s γ ( 1 ) ( t , β , γ X ; θ ) s ( 0 ) ( t , β , γ X ; θ ) d E ( N ¯ ( t ) ) , with s ( 0 ) ( t , β , γ X ; θ ) = E ( S ( 0 ) ( t , β , γ X ; θ ) ) , s β ( r ) ( t , β , γ X ; θ ) = E ( S β ( r ) ( t , β , γ X ; θ ) ) , s γ ( r ) ( t , β , γ X ; θ ) = E ( S γ ( r ) ( t , β , γ X ; θ ) ) , and N ¯ ( t ) = ∑ i ∈ I 1 N i ( t ). Under standard counting process arguments, these terms behave as approximately independent across individuals, which allows the application of the Central Limit Theorem to conclude: m ( Ω ˆ − Ω ∗ ) → d N ( 0 , V ( Ω ∗ ) ) , where V ( Ω ) = A − 1 ( Ω ) B ( Ω ) ( A − 1 ( Ω ) ) ′ , and B ( Ω ) = E ( Γ 1 ∗ ( Ω ) ⊗ 2 ). The regularity conditions required for consistency and asymptotic normality of Ω ˆ are provided in Appendix B. We also highlight that the asymptotic variance derived here does not explicitly incorporate finite-C Monte Carlo error. Instead, when Monte Carlo integration is used to evaluate H ( · ), the asymptotic theory is developed under the condition that C is sufficiently large that the resulting numerical error is o p ( n − 1 / 2 ), and hence negligible relative to sampling variability.

To estimate V ( Ω ∗ ), we use the plug-in estimator V ˆ ( Ω ˆ ) = A ˆ − 1 ( Ω ˆ ) B ˆ ( Ω ˆ ) ( A ˆ − 1 ( Ω ˆ ) ) ′ . The empirical counterpart B ˆ ( Ω ) is given by B ˆ ( Ω ) = ∑ i ∈ I 1 ∪ I 2 ∪ I 3 Γ ˆ i ∗ ( Ω ) ⊗ 2 , where Γ ˆ i ∗ ( Ω ) is similar to Γ i ∗ ( Ω ), but replaces Φ i β ∗ ( β , γ X ; θ ), and Φ i γ ∗ ( β , γ X ; θ ) with Φ ˆ i β ∗ ( β , γ X ; θ ), Φ ˆ i γ ∗ ( β , γ X ; θ ), respectively, and Φ ˆ i β ∗ ( β , γ X ; θ ) = δ i ∂ H ( β , Q i , V i ; θ ) / ∂ β H ( β , Q i , V i ; θ ) − S β ( 1 ) ( U i , β , γ X ; θ ) S ( 0 ) ( U i , β , γ X ; θ ) − ∑ j ∈ I 1 δ j Y i ( U j ) exp { γ X ′ X i } H ( β , Q i , V i ; θ ) S ( 0 ) ( U j , β , γ X ; θ ) × ∂ H ( β , Q i , V i ; θ ) / ∂ β H ( β , Q i , V i ; θ ) − S β ( 1 ) ( U j , β , γ X ; θ ) S ( 0 ) ( U j , β , γ X ; θ ) , Φ ˆ i γ ∗ ( β , γ X ; θ ) = δ i X i − S γ ( 1 ) ( U i , β , γ X ; θ ) S ( 0 ) ( U i , β , γ X ; θ ) − ∑ j ∈ I 1 δ j Y i ( U j ) exp { γ X ′ X i } H ( β , Q i , V i ; θ ) S ( 0 ) ( U j , β , γ X ; θ ) × X i − S γ ( 1 ) ( U j , β , γ X ; θ ) S ( 0 ) ( U j , β , γ X ; θ ) . That is, Φ ˆ i β ∗ and Φ ˆ i γ ∗ replaces s ( r ) ( t , β , γ X ; θ ) and E ( N ¯ ( t ) ) with their empirical estimators. By applying the law of large numbers and properties of the empirical distribution function, we have A ˆ ( Ω ˆ ) → p A ( Ω ∗ ) and B ˆ ( Ω ˆ ) → p B ( Ω ∗ ) as m → ∞, and implies that V ˆ ( Ω ˆ ) → p V ( Ω ∗ ) as m → ∞.

A ( 1 − ξ ) × 100 % pointwise confidence band for the estimated survival probability P ˆ ( T > t | z , x ; Ω ˆ ) can be constructed using the approach of [11] as follows:

Generate Ω ˆ 1 , … , Ω ˆ M for M ≫ 1, where each Ω ˆ m ∼ N ( Ω ˆ , V ˆ ( Ω ˆ ) ).

For each draw m, compute the survival probability P ˆ ( T > t | z , x ; Ω ˆ m )

Construct the lower and upper bounds of the confidence band by taking the empirical ξ / 2 and 1 − ξ / 2 quantiles of the set { P ˆ ( T > t | z , x ; Ω ˆ m ) : m = 1 , … , M }.

Numerical Study

We conducted simulation studies to evaluate the performance of the proposed generalized regression calibration approach outlined in Section 3, and to assess its performance when we approximated E Z ( exp { g ( Z ; β ) } | Q i , V i ) in (3.3) with a Taylor series approximation. The numerical study parameters were chosen to mimic the numerical study parameter specifications in [5].

We generated the true exposure Z and covariate V ≡ X from a standard bivariate normal distribution with correlation 0.2. The noisy exposure variable Q was generated according to (2.1) with parameters ( a 0 , a 1 , a 2 ′ ) = ( 0 , 0.8 , − 0.8 ) ′ and measurement error variance σ ε 2 = 1. The observed biomarker W was generated under (2.2) with measurement error variance σ e 2 = 0.25. The analysis cohort sample size was set to N 1 = | I 1 | = 10000, which measures Q, V, and an event time subject to random right-censoring. The censoring variable was generated from an exponential distribution with rate parameter 0.2; the generation of event times is described in the following paragraph. The calibration sample size was set to N 2 = | I 2 | = 500, which measures W, Q, and V. To evaluate the proposed method in which the analysis and calibration cohorts share units, we considered three settings in which proportion of units in I 2 also belong to I 1 : 0%, 50%, and 100%. When generating the calibration sample, we generated m i = 2 replicate measurements of the observed biomarker, and randomly selected 20% of subjects to compose the replicate dataset, so that N 3 = | I 3 | = 100.

The event time of interest T was generated under the hazard function shown in (2.4), where λ 0 ( t ) ≡ λ 0 ∈ { 0.02 , 0.04 , 0.08 }, and γ X = 0.1. The purpose of considering different values of λ 0 is to assess the performance of the proposed estimator across different event-rate settings, particularly as the rare-disease assumption becomes increasingly questionable. We considered the following three settings for g ( Z ; β ):

Setting 1: g ( Z ; β ) = β 1 Z + β 2 Z 2 with β = ( β 1 , β 2 ) ′ = ( − 0.2 , 0.05 ) ′ .

Setting 2: g ( Z ; β ) = β 1 Z + β 2 Z 2 + β 3 Z 3 with β = ( β 1 , β 2 , β 3 ) ′ = ( − 0.2 , 0.05 , 0.02 ).

Setting 3: g ( Z ; β ) = β 1 Z / 1 + Z 2 with β = β 1 = 0.2.

 The results were derived from 1000 Monte Carlo simulations.

Results for Quadratic Effects on the Hazard Function

For the setting with g ( Z ; β ) = β 1 Z + β 2 Z 2 , our objective is to estimate β and γ X under (2.4) with data D 1 ∪ D 2 ∪ D 3 . We considered the following estimators:

Ideal: Estimate β and γ X under (2.4) with Z observed and g ( Z ; β ) correctly specified.

Naive: Estimate β and γ X under (2.4) with Q replacing Z and g ( Z ; β ) correctly specified.

Linear Regression Calibration: Estimate β and γ X under (3.3) with g ( Z ; β ) = β 1 Z. That is, β 2 is (incorrectly) estimated with 0 for all data replicates.

Quadratic Regression Calibration: Estimate β and γ X under (3.3) with g ( Z ; β ) = β 1 Z + β 2 Z 2 .

Cubic Regression Calibration: Estimate β and γ X under (3.3) with g ( Z ; β ) = β 1 Z + β 2 Z 2 + β 3 Z 3 , where we applied a Monte Carlo approximation with C = 1000 generated samples.

 Table 1 summarizes results across simulated datasets by reporting (i) the average bias, (ii) the empirical standard error, (iii) the median of the estimated standard error (see Section 3.4), and (iv) the coverage probability under a Wald 95% confidence interval.

Summary of estimated parameters under 1000 data replicates where the event time was generated under g ( Z ; β ) = β 1 Z + β 2 Z 2 . Results are separated based on (i) the analysis and calibration dataset overlap proportion: | I 1 ∩ I 2 | / | I 2 | ∈ { 0 % , 50 % , 100 % }, and (ii) the baseline hazard parameter specification: λ 0 ∈ { 0.2 , 0.4 , 0.8 }.

Analysis and Calibration Dataset Overlap Proportion: | I 1 ∩ I 2 | / | I 2 | = 0 %
Average Event Rate ≈ 9.7 % ( λ 0 = 0.2) Average Event Rate ≈ 17.6 % ( λ 0 = 0.4) Average Event Rate ≈ 29.8 % ( λ 0 = 0.8)
Method β 1 β 2 β 3 γ X β 1 β 2 β 3 γ X β 1 β 2 β 3 γ X
Bias × 100 Ideal -0.271 -0.093 - -0.087 -0.160 -0.025 - -0.087 -0.025 0.028 - -0.034
Naive 9.771 -4.525 - -11.120 9.996 -4.491 - -10.910 10.295 -4.503 - -10.558
LRC -1.832 - - -0.108 -1.390 - - -0.149 -0.694 - - -0.146
QRC -0.131 -0.638 - -0.107 0.280 -0.390 - -0.146 0.836 -0.400 - -0.143
CRC -0.707 -1.806 0.043 -0.115 -0.220 -1.222 0.066 -0.155 0.443 -0.916 0.072 -0.152
Empirical SE Ideal 0.032 0.021 - 0.033 0.025 0.016 - 0.025 0.018 0.012 - 0.019
Naive 0.026 0.011 - 0.037 0.019 0.008 - 0.027 0.015 0.007 - 0.021
LRC 0.056 - - 0.036 0.044 - - 0.028 0.035 - - 0.022
QRC 0.055 0.051 - 0.036 0.042 0.037 - 0.028 0.033 0.030 - 0.022
CRC 0.116 0.056 0.034 0.036 0.100 0.042 0.029 0.028 0.083 0.033 0.024 0.022
Median Estimated SE Ideal 0.032 0.021 - 0.033 0.024 0.016 - 0.024 0.018 0.012 - 0.019
Naive 0.025 0.011 - 0.036 0.019 0.008 - 0.027 0.015 0.006 - 0.021
LRC 0.056 - - 0.035 0.042 - - 0.027 0.034 - - 0.021
QRC 0.054 0.048 - 0.035 0.040 0.036 - 0.027 0.032 0.028 - 0.021
CRC 0.104 0.054 0.028 0.035 0.087 0.041 0.025 0.027 0.074 0.031 0.021 0.021
Coverage Probability Ideal 94.7 95.2 - 95.4 94.3 95.0 - 95.0 95.5 95.0 - 95.9
Naive 3.7 1.1 - 13.1 < 0.1 < 0.1 - 1.7 < 0.1 < 0.1 - 0.2
LRC 93.3 - - 94.0 93.2 - - 94.6 93.9 - - 94.6
QRC 94.2 93.3 - 93.5 93.3 93.6 - 94.2 92.5 93.7 - 94.6
CRC 91.8 92.8 89.5 94.3 91.1 93.1 90.6 94.5 91.9 92.2 90.9 94.7
Analysis and Calibration Dataset Overlap Proportion: | I 1 ∩ I 2 | / | I 2 | = 50 %
Average Event Rate ≈ 9.7 % ( λ 0 = 0.2) Average Event Rate ≈ 17.6 % ( λ 0 = 0.4) Average Event Rate ≈ 29.8 % ( λ 0 = 0.8)
Method β 1 β 2 β 3 γ X β 1 β 2 β 3 γ X β 1 β 2 β 3 γ X
Bias × 100 Ideal -0.265 -0.075 - -0.082 -0.172 -0.006 - -0.057 -0.021 0.024 - -0.037
Naive 9.763 -4.507 - -11.131 10.000 -4.492 - -10.886 10.301 -4.502 - -10.556
LRC -1.870 - - -0.080 -1.359 - - -0.100 -0.685 - - -0.123
QRC -0.163 -0.597 - -0.079 0.296 -0.440 - -0.098 0.844 -0.428 - -0.120
CRC -0.710 -1.772 0.038 -0.088 -0.210 -1.254 0.068 -0.107 0.440 -0.938 0.073 -0.129
Empirical SE Ideal 0.032 0.021 - 0.033 0.025 0.015 - 0.025 0.018 0.012 - 0.018
Naive 0.025 0.011 - 0.037 0.019 0.008 - 0.027 0.015 0.007 - 0.021
LRC 0.056 - - 0.036 0.043 - - 0.027 0.035 - - 0.022
QRC 0.055 0.050 - 0.036 0.042 0.037 - 0.027 0.034 0.029 - 0.022
CRC 0.117 0.056 0.034 0.036 0.099 0.041 0.029 0.027 0.083 0.033 0.024 0.022
Median Estimated SE Ideal 0.032 0.021 - 0.033 0.024 0.016 - 0.024 0.018 0.012 - 0.019
Naive 0.025 0.011 - 0.036 0.019 0.008 - 0.027 0.015 0.006 - 0.021
LRC 0.056 - - 0.035 0.043 - - 0.027 0.034 NA - 0.021
QRC 0.054 0.048 - 0.035 0.040 0.036 - 0.027 0.032 0.028 - 0.021
CRC 0.104 0.054 0.028 0.035 0.087 0.040 0.024 0.027 0.074 0.031 0.021 0.021
Coverage Probability Ideal 94.9 95.2 - 95.0 94.2 95.2 - 94.9 95.5 95.0 - 95.9
Naive 3.4 1.3 - 13.5 < 0.1 < 0.1 - 1.7 < 0.1 < 0.1 - 0.2
LRC 93.2 - - 94.5 93.2 - - 94.4 93.1 - - 94.6
QRC 94.0 93.7 - 94.2 93.6 94.1 - 94.2 92.3 93.6 - 94.5
CRC 92.2 92.9 89.1 94.8 92.3 93.0 90.2 94.3 91.5 92.4 90.9 94.6
Analysis and Calibration Dataset Overlap Proportion: | I 1 ∩ I 2 | / | I 2 | = 100 %
Average Event Rate ≈ 9.7 % ( λ 0 = 0.2) Average Event Rate ≈ 17.6 % ( λ 0 = 0.4) Average Event Rate ≈ 29.8 % ( λ 0 = 0.8)
Method β 1 β 2 β 3 γ X β 1 β 2 β 3 γ X β 1 β 2 β 3 γ X
Bias × 100 Ideal -0.309 -0.080 - -0.108 -0.173 -0.025 - -0.066 -0.029 0.029 - -0.018
Naive 9.739 -4.510 - -11.178 9.983 -4.493 - -10.901 10.311 -4.502 - -10.533
LRC -1.929 - - -0.055 -1.411 - - -0.047 -0.680 - - -0.064
QRC -0.236 -0.620 - -0.054 0.245 -0.456 - -0.046 0.850 -0.430 - -0.060
CRC -0.915 -1.797 0.08 -0.064 -0.193 -1.298 0.045 -0.055 0.466 -0.944 0.066 -0.070
Empirical SE Ideal 0.032 0.021 - 0.034 0.025 0.015 - 0.025 0.018 0.012 - 0.019
Naive 0.026 0.011 - 0.038 0.019 0.008 - 0.027 0.015 0.007 - 0.021
LRC 0.057 - - 0.037 0.043 - - 0.027 0.035 - - 0.021
QRC 0.057 0.050 - 0.037 0.043 0.037 - 0.028 0.034 0.030 - 0.021
CRC 0.119 0.057 0.034 0.037 0.101 0.043 0.029 0.027 0.083 0.033 0.024 0.021
Median Estimated SE Ideal 0.032 0.021 - 0.033 0.024 0.016 - 0.024 0.018 0.012 - 0.019
Naive 0.025 0.011 - 0.036 0.019 0.008 - 0.027 0.015 0.006 - 0.021
LRC 0.056 - - 0.035 0.043 - - 0.027 0.034 - - 0.021
QRC 0.054 0.048 - 0.035 0.040 0.036 - 0.027 0.032 0.028 - 0.021
CRC 0.104 0.054 0.028 0.035 0.087 0.040 0.024 0.027 0.074 0.031 0.021 0.021
Coverage Probability Ideal 95.0 95.1 - 94.8 94.1 95.3 - 94.7 95.4 95.2 - 95.8
Naive 3.8 1.2 - 13.0 < 0.1 < 0.1 - 1.7 < 0.1 < 0.1 - 0.2
LRC 93.4 - - 93.9 93.6 - - 94.7 94.1 - - 95.1
QRC 93.5 93.4 - 93.5 93.3 94.2 - 94.6 91.7 93.8 - 95.1
CRC 91.8 92.8 89.8 93.9 91.1 93.3 90.2 94.6 92.1 91.7 91.2 95.1

Abbreviations: Ave.: Average; CRC: Cubic Regression Calibration; LRC: Linear Regression Calibration; QRC: Quadratic Regression Calibration; SE: Standard error.

As expected, the ideal estimator exhibited negligible bias and close-to-nominal coverage across all settings. In contrast, the naive estimator showed substantial bias and poor coverage for the exposure-effect parameters. Among the practical competitors, the quadratic regression calibration estimator performed best overall, with small bias and coverage probabilities generally close to the nominal 95% level across the different event-rate and overlap settings. These results are consistent with those reported in [5]. The cubic regression calibration estimator also performed well, although it showed slightly larger bias and mild undercoverage in some settings, likely reflecting the use of a more flexible model than necessary when the true exposure effect is quadratic. As expected, linear regression calibration was more sensitive to model misspecification, since it omits the quadratic term in g ( Z ; β ). A notable finding is that the performance of the regression calibration estimators remained relatively stable as the event rate increased. In particular, the quadratic and cubic regression calibration estimators continued to show good finite-sample performance even when the event rate approached 30%, suggesting that the proposed approach remains effective in settings where the rare-disease assumption is less tenable. In addition, performance was broadly similar across the three overlap settings, indicating that the estimator is robust to whether the calibration sample is external, partially internal, or fully internal to the analysis cohort.

 Results for Cubic Effects on the Hazard Function

For the setting with g ( Z ; β ) = β 1 Z + β 2 Z 2 + β 3 Z 3 , our objective is to estimate β and γ X under (2.4) with data D 1 ∪ D 2 ∪ D 3 . In addition to the estimators considered in Section 4.1, we considered

FOTSE: Estimate β and γ X under (3.3) where we approximate H ( β , Q i , V i ; θ ) = E Z ( exp { g ( Z ; β ) | Q i , V i } ) with a first-order Taylor series expansion (see Remark 2 of Section 3.2).

SOTSE: Estimate β and γ X under (3.3) where we approximate H ( β , Q i , V i ; θ ) = E Z ( exp { g ( Z ; β ) | Q i , V i } ) with a second-order Taylor series expansion (see Remark 2 of Section 3.2).

 Table 2 summarizes the resulting parameter estimates. The ideal estimator again displayed negligible bias and coverage probabilities close to the nominal level, while the naive estimator performed poorly due to uncorrected measurement error. More importantly, Table 2 highlights the advantage of using the proposed Monte Carlo-based regression calibration approach when the true exposure effect is more complex. In this setting, the cubic regression calibration estimator yielded substantially improved performance relative to the first- and second-order Taylor-series approximations. Although some mild undercoverage remained for certain components of β, the cubic regression calibration estimator generally showed smaller bias and more reliable coverage than the Taylor-based alternatives, particularly as the event rate increased. The deterioration of the Taylor-series approximations was especially apparent for γ X , for which coverage fell well below the nominal 95% level, particularly at higher event rates. This suggests that low-order Taylor approximations may be inadequate when the nonlinear exposure effect is highly complex. By contrast, the proposed cubic regression calibration estimator remained comparatively stable across the different event-rate settings, indicating that the Monte Carlo approximation offers a more reliable strategy for evaluating the induced hazard function in this setting.

Summary of estimated parameters under 1000 data replicates where the event time was generated under g ( Z ; β ) = β 1 Z + β 2 Z 2 + β 3 Z 3 . Results are separated based on (i) the analysis and calibration dataset overlap proportion: | I 1 ∩ I 2 | / | I 2 | ∈ { 0 % , 50 % , 100 % }, and (ii) the baseline hazard parameter specification: λ 0 ∈ { 0.2 , 0.4 , 0.8 }.

Analysis and Calibration Dataset Overlap Proportion: | I 1 ∩ I 2 | / | I 2 | = 0 %
Average Event Rate ≈ 9.6 % ( λ 0 = 0.2) Average Event Rate ≈ 17.5 % ( λ 0 = 0.4) Average Event Rate ≈ 29.7 % ( λ 0 = 0.8)
Method β 1 β 2 β 3 γ X β 1 β 2 β 3 γ X β 1 β 2 β 3 γ X
Bias × 100 Ideal -0.414 -0.122 0.021 0.008 -0.274 -0.005 0.037 -0.013 -0.150 0.024 0.031 -0.012
Naive 12.494 -4.554 -1.917 -7.522 12.693 -4.497 -1.924 -7.414 12.827 -4.489 -1.926 -7.301
LRC 5.046 - - 0.008 5.335 - - -0.038 5.619 - - -0.071
QRC 6.107 -0.786 - 0.010 6.406 -0.533 - -0.036 6.642 -0.493 - -0.067
FOTSE 3.225 -2.814 -1.032 -7.520 3.649 -2.528 -1.123 -7.413 3.956 -2.457 -1.151 -7.300
SOTSE 7.452 -1.495 -2.222 -7.476 4.345 -1.491 -1.765 -7.394 2.458 -1.498 -1.487 -7.294
CRC 1.288 -1.652 -0.582 0.008 1.353 -1.057 -0.504 -0.040 1.504 -0.732 -0.450 -0.072
Empirical SE Ideal 0.050 0.021 0.012 0.034 0.038 0.016 0.009 0.025 0.028 0.012 0.007 0.019
Naive 0.037 0.011 0.004 0.037 0.027 0.009 0.003 0.027 0.022 0.007 0.003 0.021
LRC 0.055 - - 0.036 0.042 - - 0.027 0.034 - - 0.021
QRC 0.053 0.052 - 0.036 0.040 0.038 - 0.027 0.032 0.030 - 0.021
FOTSE 0.084 0.058 0.052 0.037 0.063 0.044 0.038 0.027 0.050 0.034 0.030 0.021
SOTSE 0.219 0.062 0.033 0.037 0.191 0.047 0.028 0.027 0.165 0.035 0.025 0.021
CRC 0.117 0.055 0.033 0.036 0.098 0.040 0.028 0.027 0.081 0.031 0.023 0.021
Median Estimated SE Ideal 0.050 0.021 0.012 0.033 0.037 0.016 0.009 0.024 0.029 0.012 0.007 0.019
Naive 0.037 0.011 0.005 0.036 0.028 0.008 0.003 0.027 0.021 0.006 0.003 0.021
LRC 0.055 - - 0.035 0.041 - - 0.026 0.032 - - 0.020
QRC 0.052 0.049 - 0.034 0.038 0.036 - 0.026 0.030 0.028 - 0.020
FOTSE 0.084 0.057 0.049 0.036 0.063 0.042 0.037 0.027 0.049 0.033 0.029 0.021
SOTSE 0.235 0.058 0.036 0.037 0.187 0.043 0.029 0.027 0.151 0.034 0.023 0.021
CRC 0.103 0.052 0.027 0.035 0.087 0.039 0.024 0.026 0.073 0.030 0.020 0.020
Coverage Probability Ideal 94.5 95.2 94.8 94.9 94.5 95.2 95.0 95.4 95.1 94.1 95.3 95.8
Naive 8.4 1.4 0.9 45.7 0.2 < 0.1 < 0.1 20.6 < 0.1 < 0.1 < 0.1 5.5
LRC 85.2 - - 94.0 73.6 - - 94.8 58.1 - - 94.7
QRC 78.1 93.0 - 93.4 60.9 94.1 - 94.7 40.7 93.5 - 94.5
FOTSE 93.6 91.8 93.0 46.3 92.0 90.5 93.1 21.8 86.1 87.9 92.6 5.6
SOTSE 96.8 93.2 93.2 46.7 95.8 92.5 91.4 22.0 95.3 92.2 89.5 5.6
CRC 91.4 92.8 88.3 94.0 91.8 93.3 90.6 94.8 92.0 93.1 90.9 94.6
Analysis and Calibration Dataset Overlap Proportion: | I 1 ∩ I 2 | / | I 2 | = 50 %
Average Event Rate ≈ 9.6 % ( λ 0 = 0.2) Average Event Rate ≈ 17.5 % ( λ 0 = 0.4) Average Event Rate ≈ 29.7 % ( λ 0 = 0.8)
Method β 1 β 2 β 3 γ X β 1 β 2 β 3 γ X β 1 β 2 β 3 γ X
Bias × 100 Ideal -0.390 -0.123 0.027 0.031 -0.311 0.002 0.045 -0.043 -0.164 0.026 0.034 -0.015
Naive 12.556 -4.561 -1.921 -7.468 12.683 -4.491 -1.922 -7.443 12.825 -4.489 -1.926 -7.305
LRC 5.123 - - 0.054 5.348 - - -0.058 5.625 - - -0.059
QRC 6.165 -0.819 - 0.056 6.419 -0.539 - -0.056 6.646 -0.519 - -0.055
FOTSE 3.387 -2.814 -1.074 -7.465 3.635 -2.461 -1.094 -7.441 3.966 -2.434 -1.144 -7.304
SOTSE 8.024 -1.544 -2.296 -7.417 4.444 -1.350 -1.789 -7.421 2.563 -1.463 -1.503 -7.295
CRC 1.369 -1.668 -0.600 0.052 1.358 -1.013 -0.501 -0.060 1.499 -0.740 -0.447 -0.060
Empirical SE Ideal 0.050 0.021 0.012 0.034 0.037 0.015 0.009 0.025 0.029 0.012 0.007 0.019
Naive 0.037 0.011 0.005 0.037 0.028 0.008 0.003 0.027 0.022 0.007 0.003 0.021
LRC 0.055 - - 0.036 0.042 - - 0.027 0.034 - - 0.021
QRC 0.053 0.051 - 0.036 0.040 0.038 - 0.027 0.032 0.029 - 0.021
FOTSE 0.084 0.058 0.052 0.037 0.063 0.043 0.038 0.027 0.051 0.034 0.030 0.021
SOTSE 0.221 0.063 0.034 0.037 0.198 0.046 0.030 0.027 0.164 0.035 0.025 0.021
CRC 0.116 0.054 0.033 0.036 0.097 0.039 0.027 0.027 0.081 0.031 0.023 0.021
Median Estimated SE Ideal 0.050 0.021 0.012 0.033 0.037 0.016 0.009 0.024 0.029 0.012 0.007 0.019
Naive 0.037 0.011 0.005 0.036 0.028 0.008 0.003 0.027 0.021 0.006 0.003 0.021
LRC 0.055 - - 0.035 0.041 - - 0.026 0.032 - - 0.020
QRC 0.052 0.049 - 0.034 0.038 0.036 - 0.026 0.030 0.028 - 0.020
FOTSE 0.084 0.057 0.049 0.036 0.063 0.042 0.037 0.027 0.049 0.033 0.029 0.021
SOTSE 0.236 0.058 0.037 0.037 0.188 0.043 0.029 0.027 0.151 0.034 0.023 0.021
CRC 0.102 0.053 0.027 0.035 0.087 0.039 0.023 0.026 0.073 0.030 0.020 0.020
Coverage Probability Ideal 95.0 94.9 94.9 94.8 94.7 95.4 95.2 95.2 95.1 94.1 95.4 96.0
Naive 8.0 1.5 1.1 46.2 0.2 < 0.1 < 0.1 20.3 < 0.1 < 0.1 < 0.1 5.5
LRC 85.0 - - 94.1 74.0 - - 94.5 58.5 - - 94.9
QRC 77.1 93.6 - 93.9 59.5 94.0 - 94.5 41.0 93.5 - 94.9
FOTSE 93.2 91.5 92.8 46.8 92.0 91.0 92.6 21.5 85.8 87.6 92.1 5.6
SOTSE 96.9 92.6 93.5 47.3 95.7 92.6 91.1 21.8 94.9 91.7 89.5 5.6
CRC 91.8 93.3 88.1 94.4 92.0 93.7 90.3 94.5 91.6 93.2 91.0 94.9
Analysis and Calibration Dataset Overlap Proportion: | I 1 ∩ I 2 | / | I 2 | = 100 %
Average Event Rate ≈ 9.6 % ( λ 0 = 0.2) Average Event Rate ≈ 17.5 % ( λ 0 = 0.4) Average Event Rate ≈ 29.7 % ( λ 0 = 0.8)
Method β 1 β 2 β 3 γ X β 1 β 2 β 3 γ X β 1 β 2 β 3 γ X
Bias × 100 Ideal -0.410 -0.115 0.033 -0.027 -0.298 -0.007 0.044 -0.054 -0.168 0.023 0.034 -0.011
Naive 12.551 -4.563 -1.920 -7.520 12.685 -4.500 -1.923 -7.459 12.826 -4.487 -1.926 -7.303
LRC 5.118 - - 0.024 5.311 - - -0.029 5.602 - - -0.025
QRC 6.146 -0.851 - 0.026 6.379 -0.559 - -0.027 6.627 -0.520 - -0.022
FOTSE 3.376 -2.813 -1.053 -7.517 3.648 -2.481 -1.128 -7.457 3.975 -2.436 -1.171 -7.302
SOTSE 7.586 -1.520 -2.216 -7.473 4.221 -1.316 -1.767 -7.438 2.612 -1.506 -1.514 -7.295
CRC 1.367 -1.679 -0.598 0.021 1.279 -1.035 -0.486 -0.031 1.462 -0.743 -0.444 -0.027
Empirical SE Ideal 0.050 0.021 0.012 0.034 0.038 0.015 0.009 0.025 0.029 0.012 0.007 0.019
Naive 0.037 0.011 0.005 0.037 0.028 0.008 0.003 0.027 0.022 0.007 0.003 0.021
LRC 0.055 - - 0.036 0.042 - - 0.027 0.034 - - 0.020
QRC 0.054 0.050 - 0.036 0.041 0.038 - 0.027 0.033 0.029 - 0.020
FOTSE 0.084 0.058 0.051 0.037 0.063 0.044 0.038 0.027 0.050 0.034 0.030 0.021
SOTSE 0.219 0.061 0.033 0.037 0.189 0.046 0.029 0.027 0.159 0.035 0.024 0.021
CRC 0.118 0.054 0.033 0.036 0.099 0.040 0.028 0.027 0.081 0.031 0.023 0.020
Median Estimated SE Ideal 0.050 0.021 0.012 0.033 0.037 0.016 0.009 0.024 0.029 0.012 0.007 0.019
Naive 0.037 0.011 0.005 0.036 0.028 0.008 0.003 0.027 0.021 0.006 0.003 0.021
LRC 0.055 - - 0.035 0.041 - - 0.026 0.032 - - 0.020
QRC 0.052 0.048 - 0.034 0.039 0.036 - 0.026 0.030 0.028 - 0.020
FOTSE 0.084 0.057 0.049 0.036 0.063 0.042 0.037 0.027 0.049 0.033 0.029 0.021
SOTSE 0.237 0.058 0.037 0.037 0.187 0.043 0.028 0.027 0.152 0.034 0.023 0.021
CRC 0.103 0.053 0.027 0.035 0.086 0.039 0.023 0.026 0.073 0.030 0.020 0.020
Coverage Probability Ideal 94.6 95.1 94.8 94.6 94.6 95.2 95.0 95.0 95.1 94.2 95.3 95.9
Naive 8.0 1.4 1.0 45.9 0.2 < 0.1 < 0.1 20.2 < 0.1 < 0.1 < 0.1 5.5
LRC 84.8 - - 94.2 73.7 - - 95.0 59.2 - - 95.8
QRC 77.4 93.8 - 93.6 59.9 94.2 - 94.9 42.0 93.3 - 95.7
FOTSE 93.3 91.3 93.3 46.5 92.2 90.3 93.5 21.4 86.5 87.9 92.6 5.5
SOTSE 97.0 92.8 93.6 47.1 96.0 93.5 90.7 21.7 94.9 92.0 89.4 5.5
CRC 91.4 93.8 88.8 93.9 91.4 93.7 89.4 95.0 91.8 93.0 91.4 95.8

Abbreviations: Ave.: Average; CRC: Cubic Regression Calibration; FOTSE: First Order Taylor Series Expansion; LRC: Linear Regression Calibration; QRC: Quadratic Regression Calibration; SE: Standard error; SOTSE: Second Order Taylor Series Expansion.

As in Table 1, the results were broadly stable across the different overlap proportions. Thus, the main message of Table 2 is not only that the proposed regression calibration approach can accommodate more complex nonlinear exposure effects, but also that its performance is not unduly sensitive to whether the calibration sample is fully external or shares subjects with the analysis cohort.

 Results for Non-Linear Effect on the Hazard Function

For the setting with g ( Z ; β ) = β 1 Z / 1 + Z 2 , we fixed λ 0 = 0.02, which yielded an average uncensored event rate of approximately 9.19% across data replicates. In this scenario, our focus shifts from estimating β and γ X to estimating the conditional survival function. For illustrative purposes, we showcase the estimation of the conditional survival functions with x = z = 0, which corresponds to their marginal expectations. To capture the nonlinear effect of Z on the hazard function, we approximate g ( Z ; β ) with a Bernstein polynomial of degree 4 [17]. We then estimate the survival function by approximating H ( β , Q i , V i ; θ ) = E Z ( exp { β 1 Z / 1 + Z 2 } | Q i , V i ) using three methods:

a first-order Taylor series expansion, as used in Section 4.2

a second-order Taylor series expansion, as used in Section 4.2

the proposed Monte Carlo approximation with C = 1000 generated samples.

 Figure 1 illustrates the simulation results. While the Taylor series expansions offer substantially lower computational cost compared to Monte Carlo integration, they fail to adequately approximate H ( β , Q i , V i ; θ ), as evidenced by the noticeable bias in the resulting estimates. In contrast, the Bernstein polynomial of degree 4 specification with a Monte Carlo approximation successfully captures the nonlinear effect of Z. Alternative basis functions, such as B-splines or natural cubic splines, could also be employed.

Averaged conditional survival functions with x = z = 0 from a Bernstein polynomial of degree 4 and coverage probabilities of the estimated survival functions, where we approximated H ( β , Q i , V i ; θ ) = E Z ( exp { β 1 Z / 1 + Z 2 } | Q i , V i ) with a first-order Taylor series expansion (FOTSE), a second-order Taylor series expansion (SOTSE), and a Monte Carlo approximation with C = 1000 (GRC_MC). The 95% pointwise confidence band for the survival functions were computed with the procedure outlined in Section 3.4.

 R Package Overview

We have developed an R package, NonLinearRC, to implement regression calibration with Monte Carlo integration, as described in Section 3. The package is publicly available on CRAN and GitHub (https://github.com/drtjt/NonLinearRC). We proceed to demonstrate its use with a simulated dataset.

Following the notation introduced in Section 2, we assume that users have datasets D 1 , D 2 , and D 3 available, which we will refer to as the data frames df_Analysis, df_Biomarker, and df_Replicate, respectively. As an illustration, these data frames might appear as follows:

where we assume that units i ∈ I 3 have two replicate observations (i.e., m i ≡ 2 for all i ∈ I 3 ). Note that the unit identification variable, ID, shows that some units are shared between df_Biomarker and df_Replicate, whereas these units are distinct from those in df_Analysis. All other variable names are consistent with the definitions presented in Section 2.

The main function in NonLinearRC is NL_RC. It takes the following arguments as inputs:

df1: Data frame containing right-censored observations of T, used to estimate parameters in (2.4).

df2: Data frame used to estimate parameters in (2.1).

df3: Data frame used to estimate measurement error variance σ e 2 in (2.2).

df1_ID_index = NA: Column number for the identification variable ID in df1. If not specified (default), it is assumed that no units are shared between df1 and df2/df3.

df1_U_index: Column number for the right-censored event time U in df1.

df1_delta_index: Column number for the event time indicator δ in df1.

df1_Q_index: Column number for the observed exposure variable Q in df1.

df1_X_index: Column number(s) in df1 for any additional covariates included in model (2.4). Note that categorical variables need to be transformed to dummy variables.

df2_ID_index = NA: Column number for the identification variable ID in df2. If not specified (default), it is assumed that no units are shared between df2 and df1/df3.

df2_Response_index: Column number for the observed biomarker W in df2 under (3.1).

df2_Covariates_index: Column number(s) in df2 for covariates under (3.1). Note that categorical variables need to be transformed to dummy variables.

df3_ID_index = NA: Column number for the identification variable ID in df3. If not specified (default), it is assumed that no units are shared between df3 and df1/df2.

Function_Type = "Polynomial": Specifies the function form of g ( Z ; β ). Users can select one of the following:

Polynomial: Degree DF polynomial of Z.

Bernstein: Linear combination of DF Bernstein polynomial basis functions.

B-Spline: Linear combination of DF B-spline basis functions.

C-Spline: Linear combination of DF C-spline basis functions.

I-Spline: Linear combination of DF I-spline basis functions.

M-Spline: Linear combination of DF M-spline basis functions.

NC-Spline: Linear combination of DF natural cubic spline basis functions.

DF = 1: Degree of the polynomial for g ( Z ; β ), or number of basis functions for splines.

Method = "Cox": Specifies the estimation procedure under (2.4) or (3.3). Users can select one of the following:

Cox: Estimate parameters under (2.4) with Q replacing Z (default).

Linear: Estimate parameters under (3.3) with g ( Z ; β ) = β Z using regression calibration.

Quadratic: Estimate parameters under (3.3) with g ( Z ; β ) = β 1 Z + β 2 Z 2 , using the regression calibration procedure of [5].

General_1: Estimate parameters under (3.3) using a first-order Taylor series expansion of exp { g ( Z ; β ) } around Z = α 0 + α Q + α 2 ′ V.

General_2: Estimate parameters under (3.3) using a second-order Taylor series expansion of exp { g ( Z ; β ) } around Z = α 0 + α Q + α 2 ′ V.

General_MC: Estimate parameters under (3.3) using a Monte Carlo approximation to E Z ( exp { g ( Z ; β ) } ∣ Q i , V i ).

EstimateNuisance_Logi: Logical indicator specifying whether σ e 2 and parameters in (3.1) should be estimated.

mu_vec = NULL: Vector of E ( Z | Q i , V i ) for each unit in df1; if provided, these values are used instead of estimating them.

Variance_vec = NULL: Vector of V a r ( Z | Q i , V i ) for each unit in df1; if provided, these values are used instead of estimating them.

t_grid: Vector of time points at which to estimate the conditional survival function.

z_grid: Vector of Z values at which to evaluate the conditional survival function.

C = 1000: Number of Monte Carlo samples for estimating β ˆ and/or its standard error.

seed = 235346: Random number seed.

x = 0: Value of X at which the conditional survival function is evaluated.

Spline_Z_vec = NULL: Vector of Z values for generating spline basis functions.

Boundary.knots = NULL: Boundary points for the spline basis functions.

Estimate_CI_Surv = TRUE: Logical indicator specifying whether a 95% pointwise Monte Carlo confidence interval is computed, following [11].

 To illustrate, we use NL_RC to estimate parameters in (3.3) with g ( Z ; β ) = β 1 Z + β 2 Z 2 + β 3 Z 3 . With our datasets, this looks like:

The output of NL_RC is a list containing the following components:

Estimates: Estimates for β and γ X under (2.4).

Nuisance_Estimates: Estimated model parameters of (3.1) and σ ˆ e 2 .

Baseline_Hazard: Data frame containing the estimated cumulative baseline hazard function Λ ˆ 0 ( · ) at each observed uncensored event time.

SurvProb_zgrid_list: Estimated conditional survival function, with each element of the list corresponding to a fixed value of Z in z_grid at times t_grid.

SurvProb_zgrid_CI_list: 95% point-wise confidence interval of the conditional survival function for each value in z_grid at times t_grid (if requested).

SurvProb_zgrid_Var_list: Estimated variance of the conditional survival function for each value in z_grid at times t_grid (if requested).

 In practice, we recommend beginning with relatively simple spline representations, for example natural cubic splines with 3–5 degrees of freedom, in order to balance flexibility and numerical stability. Knot placement may be based on empirical quantiles of the observed exposure surrogate, and the number of knots may be selected using model selection criteria such as AIC or BIC [4] More generally, we recommend fitting a small sequence of candidate basis specifications and selecting the simplest representation that adequately captures the apparent exposure and response relationship, while monitoring numerical stability.

 Conclusion

In this paper, we developed a method for estimating nonlinear covariate effects in a Cox regression model when the covariate is subject to measurement error, assuming a sample of objective biomarkers is available for correcting the covariate measurement error. Our approach generalizes the work of [5] by allowing for arbitrary nonlinear effects through Monte Carlo integration, and we proposed a sandwich variance estimator to quantify the uncertainty of the regression estimates. Simulation studies demonstrate that our method accurately captures nonlinear effects; capabilities that are not offered by existing alternatives. While Monte Carlo integration is computationally intensive, our results indicate that simplifying this step with a Taylor series approximation can lead to biased estimates. To facilitate practical use, we have implemented our method in the R package NonLinearRC.

As shown in (3.2) and (3.3), the rare disease assumption makes the induced-hazard approximation tractable and therefore estimable within the proposed framework. Although our simulation study suggests that the resulting approach remains effective for uncensored event proportions up to roughly 30%, its performance will in general depend on the underlying data-generating mechanism and the complexity of the true exposure effect. When the rare disease assumption is more substantially violated, the regression calibration estimators may become inadequate, and alternative methods that avoid this assumption are required for estimation of the parameters in (2.4). One approach is to adopt risk set regression calibration [19, 14], in which H ( β , Q , V ; θ ) would need to be re-estimated for each risk set. Another approach is the expected estimating equation method proposed by [16] that bases its estimation of parameters in (2.4) with the expected value of the partial score function.

Although the present paper focused on a linear-Gaussian measurement error model, the broader framework is driven by the need to characterize the conditional distribution of the true exposure given covariates, namely Z | { Q , V }. Consequently, the proposed approach is not conceptually limited to the Gaussian setting. Extensions to alternative parametric or semiparametric measurement error models may be possible when this conditional distribution can be validly estimated, approximated, or sampled from. Although the required modelling makes a variety of parametric assumptions, the induced hazard model in (3.3) should be viewed as a working model. More specifically, the assumptions and model specifications in Section 2.2 used to characterize Z | { Q , V } should be assessed, to the extent possible, using the calibration data. Potential violations of linearity for the mean in model (3.1) may be examined by comparing the fitted calibration model with more flexible alternatives, such as models including interaction or higher-order terms in Q and V, and by inspecting residual-versus-fitted plots. Heteroscedasticity may be assessed using residual diagnostics and, when replicate biomarker measurements are available, by examining whether within-subject variability changes with Q or components of V. Departures from conditional normality may be investigated using residual histograms or normal quantile plots. The nondifferential error assumption is less directly testable and should primarily be justified on scientific and design grounds. If V is not sufficiently rich to capture important determinants of the relationship between the true exposure and the observed measurements, then the working model for Z | { Q , V } may be misspecified, which can bias the approximation of H ( β , Q , V ; θ ) and, consequently, the estimated exposure effect parameters.

As noted by an anonymous referee, there has been increasing methodological work on settings involving a main study together with an external validation study (e.g., [20]), in which the conditional distribution Z | { Q , V } may be characterized using data from the external validation study. Our setting differs from this framework, since our auxiliary data consist of biomarker and replicate measurements rather than direct observations of the true exposure in a validation sample. Nevertheless, a related extension of the present framework when external validation data are available may be possible in such settings, provided that the relationship used to characterize Z | { Q , V } is transportable to the analysis cohort.

 Estimating Functions for Nuisance Parameters

We estimate the nuisance parameter vector θ = ( α 0 , α 1 , α 2 , σ 2 + σ e 2 ) ′ as the solution to the system of estimating equations Ψ ( θ ) Ψ α 0 Ψ α 1 Ψ α 2 Ψ σ 2 + σ e 2 = 0 , where Ψ α 0 = 1 n 2 ∑ i ∈ I 2 W i 1 − ( α 0 + α 1 Q i + α 2 ′ V i ) Ψ α 1 = 1 n 2 ∑ i ∈ I 2 ( W i 1 − ( α 0 + α 1 Q i + α 2 ′ V i ) ) Q i Ψ α 2 = 1 n 2 ∑ i ∈ I 2 ( W i 1 − ( α 0 + α 1 Q i + α 2 ′ V i ) ) V i Ψ σ 2 + σ e 2 = 1 n 2 ∑ i ∈ I 2 { ( W i 1 − ( α 0 + α 1 Q i + α 2 ′ V i ) ) 2 − ( n 2 − p ) n 2 ( σ 2 + σ e 2 ) } , with n 2 denoting the number of participants in I 2 and p = dim ( α ).

To estimate σ e 2 , we use the replicate biomarker measurements available in D 3 , and obtain σ ˆ e 2 = 1 n 3 ∑ i ∈ I 3 1 m i − 1 ∑ j = 1 m i ( W i j − W ¯ i ) 2 , where W ¯ i = ∑ j = 1 m i W i j / m i is the within-subject average and n 3 is the number of participants in I 3 . This estimator satisfies the estimating equation Ψ σ e 2 = 1 n 3 ∑ i ∈ I 3 1 m i − 1 ∑ j = 1 m i ( W i j − W ¯ i ) 2 − σ e 2 = 0 . By stacking all estimating functions, we have Ψ ( θ ) = ∑ i ∈ I 2 ∪ I 3 Ψ i ( θ ), where each participant contributes Ψ i ( θ ) = 0 for all θ, if they are not in the relevant subsample (i.e., individuals in I 2 ∖ I 3 do not contribute to estimating σ e 2 , and individuals in I 3 ∖ I 2 do not contribute to the regression parameters in θ). Finally, we obtain an estimate of σ 2 with σ ˆ 2 = σ 2 + σ e 2 ˆ − σ ˆ e 2 .

 Regularity Conditions for Asymptotic Results

Let Ω = ( β ′ , γ X ′ , θ ′ ) ′ ∈ O ⊂ R q denote the full parameter vector, where q is the dimension of Ω, and let Γ ( Ω ) = ∑ i ∈ I 1 ∪ I 2 ∪ I 3 Γ i ( Ω ) = ∑ i ∈ I 1 ∪ I 2 ∪ I 3 Φ i β ( β , γ X ; θ ) Φ i γ ( β , γ X ; θ ) Ψ i ( θ ) denote the stacked estimating function. We assume the following regularity conditions.

(A1): Independent Sampling. The data are independent across distinct study units in I 1 ∪ I 2 ∪ I 3 .

(A2): Parameter Space. The parameter space O is compact, and the target parameter Ω ∗ lies in the interior of O.

(A3): Identifiability. The population estimating equation E ( Γ ( Ω ) ) = 0 has a unique solution at Ω = Ω ∗ .

(A4): Smoothness. For each i, the map Ω ↦ Γ i ( Ω ) is continuously differentiable in an open neighbourhood of Ω ∗ , denoted by N ( Ω ∗ ). Moreover, there exists an integrable random variable M i such that sup Ω ∈ N ( Ω ∗ ) ‖ ∂ Γ i ( Ω ) / ∂ Ω ′ ‖ ≤ M i with E ( M i ) < ∞.

(A5): Finite Moments. For all Ω ∈ N ( Ω ∗ ), E ( ‖ Γ i ∗ ( Ω ) ⊗ 2 ‖ ) < ∞, and, if needed for the Taylor expansion and uniform convergence arguments, higher moments exist in N ( Ω ∗ ).

(A6): Independent Censoring. Conditional on the covariates entering the model, the censoring time C i is independent of the event time T i . In addition, the risk-set quantities are well defined over the time interval of interest.

(A7): Regularity of Cox Risk-set Quantities. For r = 0 , 1 , 2 , …, we assume that sup t ‖ S β ( r ) ( t , β , γ X ; θ ) − s β ( r ) ( t , β , γ X ; θ ) ‖ → p 0 sup t ‖ S γ ( r ) ( t , β , γ X ; θ ) − s γ ( r ) ( t , β , γ X ; θ ) ‖ → p 0 for every Ω ∈ N ( Ω ∗ ). We also require inf t s β ( 0 ) ( t , β , γ X ; θ ) > 0 , inf t s γ ( 0 ) ( t , β , γ X ; θ ) > 0 , for every Ω ∈ N ( Ω ∗ ).

(A8): Regularity of the Induced Hazard Component. The function H ( β , Q , V ; θ ) = E ( exp { g ( Z ; β ) } | Q , V ) is well defined in N ( Ω ∗ ), and its first derivative with respect to Ω exists and is continuous.

(A9): Nonsingularity of the Sensitivity Matrix. A ( Ω 0 ) = ∂ E ( Γ ( Ω ) ) / ∂ Ω ′ | Ω = Ω ∗ exists, is finite, and is nonsingular.

(A10): Positive Definiteness of the Variability Matrix. The matrix B ( Ω ∗ ) = E ( Γ i ∗ ( Ω ∗ ) ⊗ 2 ) exists and is finite.

(A11): Rare Disease Assumption. The rare-disease approximation used to derive the induced hazard is sufficiently accurate so that Ω ∗ is interpreted as the true parameter corresponding to the approximating model.

(A12): Monte Carlo Approximation. If Monte Carlo integration is used to evaluate H ( · ), the number of Monte Carlo draws C is taken sufficiently large that the resulting numerical error is o p ( n − 1 / 2 ).

 References Andersen, P. K. and Gill, R. D. (1982). Cox’s Regression Model for Counting Processes: A Large Sample Study. The Annals of Statistics 10(4) 1100–1120. MR0673646 Carroll, R. J., Ruppert, D., Stefanski, L. A. and Crainiceanu, C. M. (2006) Measurement Error in Nonlinear Models: A Modern Perspective - Second Edition. Chapman & Hall/CRC. MR2243417 Green, P. J. and Silverman, B. W. (1994) Nonparametric Regression and Generalized Linear Models. Chapman and Hall, London. MR1270012 Huang, J. Z. and Liu, L. (2006). Polynomial Spline Estimation and Inference of Proportional Hazards Regression Models with Flexible Relative Risk Form. Biometrics 62 793–802. MR2247208 Huang, Y. and Prentice, R. L. (2025). Biomarker-assisted reporting in nutritional epidemiology: addressing measurement error in exposure-disease associations. Biostatistics 26(1) 014. MR4917684 Huang, Y., van Horn, L., Tinker, L. F., Neuhouser, M. L., Carbone, L., Mossavar-Rahmani, Y., Thomas, F. and Prentice, R. (2013). Measurement error corrected sodium and potassium intake estimation using 24-hour urinary excretion. Hypertension 63 238–244. Huang, Y., Zheng, C., Tinker, L. F., Neuhouser, M. L. and Prentice, R. L. (2022). Biomarker-based methods and study designs to calibrate dietary intake for assessing diet–disease associations. The Journal of Nutrition 152(3) 899–906. Kim, S., Li, Y. and Spiegelman, D. (2014). A semiparametric copula method for Cox models with covariate measurement error. Lifetime Data Analysis 22(1) 1–16. MR3447180 Lin, D. Y. and Wei, L. J. (1989). The Robust Inference for the Cox Proportional Hazards Model. Journal of the American Statistical Association 84(408) 1074–1078. MR1134495 Neuhouser, M. L., Tinker, L., Shaw, P. A., Schoeller, S. A. Bingham, Van Horn, L., Beresford, S. A. A., Caan, B., Thomson, C., Satterfield, S., Kuller, L., Heiss, G., Smit, E., Sarto, G., Ockene, J., Stefanick, M. L., Assaf, A., Runswick, S. and Prentice, R. L. (2008). Use of recovery biomarkers to calibrate nutrient consumption self-reports in the Women’s Health Initiative. American Journal of Epidemiology 167 1247–1259. Preacher, K. J. and Selig, J. P. (2012). Advantages of Monte Carlo confidence intervals for indirect effects. Communication Methods and Measures 6 77–98. Prentice, R. L. (1982). Covariate Measurement Errors and Parameter Estimation in a Failure Time Regression Model. Biometrika 69(2) 331–342. MR0671971 Prentice, R. L., Mosavar-Rahmani, Y., Huang, Y., van Horn, L., Beresford, S. A. A., Caan, B., Tinker, L., Schoeller, D., Bingham, S., Eaton, C. B., Thomson, C., C., J. K., Ockene, J., Sarto, G. and Heiss, G. (2011). Evaluation and comparison of food records, recalls, and frequencies for energy and protein assessment by using recovery biomarkers. American Journal of Epidemiology 174 591–603. Shaw, P. A. and Prentice, R. L. (2012). Hazard Ratio Estimation for Biomarker-Calibrated Dietary Exposures. Biometrics 68 397–407. MR2959606 Sleeper, L. A. and Harrington, D. P. (1990). Regression Splines in the Cox Model With Application to Covariate Effects in Liver Disease. Journal of the American Statistical Association 85(412) 941–949. Wang, C. Y. and Pepe, S. (2000). Expected Estimating Equations to Accommodate Covariate Measurement Error. Journal of the Royal Statistical Society, Series B 62(3) 509–524. MR1772412 Wang, W. and Yan, J. (2021). Shape-Restricted Regression Splines with R Package splines2. https://doi.org/10.6339/21-JDS1020. Women’s Health Initiative Study Group (1998). Design of the Women’s Health Initiative Clinical Trial and Observational Study. Controlled Clinical Trials 19(1) 61–109. https://doi.org/10.1016/S0197-2456(97)00078-0. Xie, S. X., Wang, C. Y. and Prentice, R. L. (2001). A Risk Set Calibration Method for Failure Time Regression by Using a Covariate Reliability Sample. Journal of the Royal Statistical Society. Series B (Statistical Methodology) 63(1) 855–870. MR1872071 Yang, C., Zhang, N., Li, J., Mehta, U. V., Hart, J. E. and Spiegelman, D. L. (2025). Exposure Measurement Error Correction in Longitudinal Studies With Discrete Outcomes. Statistics in Medicine 44 70191. MR4939014 Yi, G. Y. (2017) Statistical Analysis with Measurement Error or Misclassification: Strategy Method and Application. Springer. MR3676914