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.

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;\boldsymbol{\beta })\}|{Q_{i}},{\boldsymbol{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 $\boldsymbol{V}\equiv \boldsymbol{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}},{\boldsymbol{a}^{\prime }_{2}})={(0,0.8,-0.8)^{\prime }}$ and measurement error variance ${\sigma _{\varepsilon }^{2}}=1$. The observed biomarker W was generated under (2.2) with measurement error variance ${\sigma _{e}^{2}}=0.25$. The analysis cohort sample size was set to ${N_{1}}=|{\mathcal{I}_{1}}|=$ 10000, which measures Q, $\boldsymbol{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}}=|{\mathcal{I}_{2}}|=500$, which measures W, Q, and $\boldsymbol{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 ${\mathcal{I}_{2}}$ also belong to ${\mathcal{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}}=|{\mathcal{I}_{3}}|=100$.

The event time of interest T was generated under the hazard function shown in (2.4), where ${\lambda _{0}}(t)\equiv {\lambda _{0}}\in \{0.02,0.04,0.08\}$, and ${\gamma _{X}}=0.1$. The purpose of considering different values of ${\lambda _{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;\boldsymbol{\beta })$: The results were derived from 1000 Monte Carlo simulations.

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 ${\mathcal{D}_{1}}$, ${\mathcal{D}_{2}}$, and ${\mathcal{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\in {\mathcal{I}_{3}}$ have two replicate observations (i.e., ${m_{i}}\equiv 2$ for all $i\in {\mathcal{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: To illustrate, we use NL_RC to estimate parameters in (3.3) with $g(Z;\boldsymbol{\beta })={\beta _{1}}Z+{\beta _{2}}{Z^{2}}+{\beta _{3}}{Z^{3}}$. With our datasets, this looks like:

The output of NL_RC is a list containing the following components: 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.

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(\boldsymbol{\beta },Q,\boldsymbol{V};\boldsymbol{\theta })$ 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,\boldsymbol{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,\boldsymbol{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 $\boldsymbol{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 $\boldsymbol{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 $\boldsymbol{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,\boldsymbol{V}\}$ may be misspecified, which can bias the approximation of $H(\boldsymbol{\beta },Q,\boldsymbol{V};\boldsymbol{\theta })$ 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,\boldsymbol{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,\boldsymbol{V}\}$ is transportable to the analysis cohort.