Semiparametric Bayesian analysis of censored linear regression with errors-in-covariates

Abstract

The accelerated failure time (AFT) model is a well-known alternative to the Cox proportional hazard model for analyzing time-to-event data. In this paper we consider fitting an AFT model to right censored data when a predictor variable is subject to measurement errors. First, without measurement errors, estimation of the model parameters in the AFT model is a challenging task due to the presence of censoring, especially when no specific assumption is made regarding the distribution of the logarithm of the time-to-event. The model complexity increases when a predictor is measured with error. We propose a non-parametric Bayesian method for analyzing such data. The novel component of our approach is to model (1) the distribution of the time-to-event, (2) the distribution of the unobserved true predictor, and (3) the distribution of the measurement errors all non-parametrically using mixtures of the Dirichlet process priors. Along with the parameter estimation we also prescribe how to estimate survival probabilities of the time-to-event. Some operating characteristics of the proposed approach are judged via finite sample simulation studies. We illustrate the proposed method by analyzing a data set from an AIDS clinical trial study.

Keywords

Buckley–James estimator Dirichlet process prior functional approach measurement errors mixture distributions posterior inference

1 Introduction

Right censored time-to-event data are often analyzed by fitting a Cox proportional hazard (CPH) model. Although fitting the CPH model and obtaining the estimate of the relative risk parameters via the partial likelihood method are easy, model parameter interpretation requires the understanding of instantaneous hazard. On the other hand, the accelerated failure time (AFT) model is easy to interpret. In the AFT model, the logarithm of the time-to-failure T is assumed to be a linear function of the covariates and an error term which is assumed to be free from the covariates. That means, for the AFT model,

\log (T) = β_{1}^{T} Z + β 2 X + e

(1)

where the error e is assumed to follow a distribution with finite variance and is assumed to be independent of the covariates

(Z T, X) T

. Here Z is assumed to be a vector of error-free covariates while the continuous scalar covariate X is not observed in the data. Instead, multiple replications of an erroneous unbiased surrogate W for X are observed in the data. By unbiasedness we mean E(W|X) = X, and by surrogate we mean

f (T | W, X, Z) = f (T | X, Z)

(Carroll¹). In the error-free case (i.e. when X is accurately observed) the regression parameters

β = (β_{1}^{T}, β 2) T

are difficult to estimate due to the presence of censoring especially when the distribution of e is left unspecified. There are several choices for fitting an AFT model to right censored data, such as Buckley–James estimating equations,² modified Buckley–James estimating equations proposed by Lai and Ying,³ some more recent approaches proposed by Lin and his co-authors,^4,5 and the empirical likelihood approach of Zhou and Li.⁶ Although our interest is in the semiparametric AFT model where e is left unspecified, the AFT model can be fitted assuming some flexible parametric models (generalized Gamma, log-logistic, splines, etc.) for e (Cox et al.⁷).

In the error-free AFT model context, Christensen and Johnson⁸ first considered the Dirichlet process (DP) prior for nonparametric modeling of the time-to-event, and proposed an elegant semi-Bayesian approach for estimating survival curves and the finite dimensional regression coefficient. Later, Kuo and Mallick⁹ considered a mixture of the DP prior on e, and Walker and Mallick¹⁰ proposed to use the Polya tree prior on e and a noninformative prior on the regression coefficients β. The last two papers considered full Bayesian inferences using the Markov chain Monte Carlo (MCMC) method.

In this paper we consider fitting an AFT model (1) to right censored data when the scalar covariate X is measured with error and with repeated measurements at the baseline. The motivation comes from a clinical study on AIDS. One of the important indicators for the time to AIDS or death of HIV-infected people is the CD4 count at the baseline examination before any treatment starts. The true CD4 count cannot be measured. Therefore, multiple measurements of a surrogate variable for CD4 count at the baseline are considered as the erroneous measurements of the true CD4 count. The goal is to estimate the regression coefficients utilizing the erroneous measurements for CD4 counts. While errors-in-covariate are a common issue in clinical or observational studies, fitting an AFT model when the predictor is measured with error has received little attention from researchers. He et al.¹¹ proposed a simulation and extrapolation (SIMEX) approach for estimating model parameters when the time-to-event data are subject to right censoring and a covariate is measured with error. They assumed that (1) the distribution of e belongs to a known parametric family, and (2) the errors associated with the covariate follow a normal distribution. These assumptions limit the application of their SIMEX method. Another paper in this context is by Ma and Yin.¹² They considered a broader issue by proposing a novel method of handling covariate measurement errors in a semiparametric quantile regression model. However, they require that the censoring mechanism and the actual time-to-event are marginally independent.

In order to circumvent these issues we propose a general method where (1) we do not make any parametric assumption regarding the distribution of e, (2) we do not make any parametric assumption regarding the distribution of the unobserved covariate, (3) we do not make any parametric assumption regarding the distribution of the measurement errors U in W. All of these three issues are handled by a novel application of the non-parametric Bayesian methods. In particular, in a likelihood framework, the distributions of e, X, and U are modeled non-parametrically using a mixture of a finite-dimensional Dirichlet process (FDDP), a special case of stick-breaking prior.¹³ In addition to a non-parametric modeling of e, since our approach does not make any parametric assumption regarding the distributions of X and U, the method can be considered as a functional approach in view of the modern measurement error literature. Since we use a parametric prior on the unknown regression parameter β along with nonparametric prior models for the distributions of e, X, and U, we call the proposed method semiparametric. The novelty of the proposed approach lies in the robustness of the procedure through non-parametric modeling of several nuisance densities. When a distribution is modeled by a DP mixture of kernel densities (we have taken them to be normal kernels), the distribution is essentially modeled by a mixture of infinitely many kernel densities, where the mixing proportions and the parameters of the kernel densities are random. This structure of the prior model for a density, in principle, leads to a posterior that is weakly consistent for the true density (Theorems 5.6.1–5.6.3 of Ghosh and Ramamoorthi¹⁴). This posterior consistency not only holds when the true density is a mixture of normals, but also when the true density has a compact support, such as the uniform distribution. In our set-up, instead of using a DP, for computational convenience we use a FDDP as a close approximation of the DP. Since we are modeling three nuisance distributions non-parametrically, our results are generally robust towards the distributions of e, U, and X. In the simulation studies, we numerically show the robustness of the proposed method by considering different types of distributions for e, X, and U, and comparing with some partly semiparametric approaches. In the partly semiparametric methods, one of three nuisance infinite-dimensional parameters is treated parametrically, and the results show that lack of proper modeling of at least one nuisance parameter may result in biased estimates of the regression parameters.

Previously, Müller and Roeder¹⁵ used a non-parametric Bayesian approach for handling errors in a covariate in case–control studies that do not involve censored subjects. Gustafson et al.¹⁶ considered a parametric Bayesian method for handling errors in a covariate in case–control studies. Further, Sinha et al.¹⁷ considered a non-parametric Bayesian approach for handling errors in a covariate in the logistic regression model while the effect of the covariate was modeled as a non-parametric function. However, to the best of the authors’ knowledge, our current problem is unique in that no one has addressed it before. Overall, our non-parametric Bayesian approach is useful not only for estimating the regression parameters β, but also for estimating the survival probabilities and the quantiles of the failure time distribution.

A brief outline of the remainder of the article is as follows. Section 2 contains basic models and assumptions. Section 3 discusses likelihood and priors. Posterior computation and parameter estimation are given in Section 4. Section 5 outlines some other statistical inferences using the posterior samples. Sections 6 and 7 are devoted to simulation studies and the analysis of a real data set from an AIDS clinical trial study, respectively. Concluding remarks are given in Section 8. The details of the MCMC steps and some further data analysis are relegated to the appendix.

2 Basic models and assumptions

Suppose we observe the data $(V i, Δ i, W ij, j = 1, \dots, m, Z i), i = 1, \dots, n$ , where $V i = min (T i, C i)$ , and the time-to-failure T_i is assumed to be independent of the censoring time C_i conditional on the observed covariates $(W i 1, \dots, W im, Z i)$ and the binary variable $Δ i = I (T i \leq C i)$ denotes the censoring indicator. For non-parametric modeling of the measurement error distribution we require the number of replications to be at least two (i.e. m ≥ 2). Even for handling a more restrictive scenario, such as a symmetric error distribution, one needs m ≥ 2 to identify the error distribution.¹⁸ Without repeated measurements on W, one needs to specify the distribution of the error for any structural or functional approach. We assume that T_i follows model (1), and $e ~ F e$ which is unknown. Furthermore, assume that Z_i is a vector of error-free covariates, and the surrogate variable $W_{i}^{T} = (W i 1, \dots, W im)$ is related to the unobserved latent variable X_i via the classical additive measurement error model

W ij = X i + U ij, for j = 1, \dots, m, m \geq 2

where U_ij are independent and identically distributed (i.i.d.) following a mean zero distribution F_U with a finite variance and are independent of

(V i, Δ i, X i, Z i)

. Furthermore, conditional on Z the unobserved X is assumed to follow a distribution

F X (\cdot | Z)

which is also unknown. It is known that the naive analysis of the data by replacing X_i by

W ¯ i = \sum_{j = 1}^{m} W ij / m

will, in principle, yield a biased estimator of β, and consequently the estimator of the survival function is biased.^12,19

It is worth mentioning that without measurement errors and assuming the response variables are subject to only right censoring, the Buckley–James estimator of β is obtained by solving

S (β; V, X, Z) = \sum_{i = 1}^{n} (Z i X i) {T ~ i (β) - ~ ¯ T (β) - (Z i - Z ¯) T β 1 - (X i - X ¯) β 2} = 0

where

T ~ i (β) = Δ i \log (V i) + (1 - Δ i) (Z_{i}^{T} β 1 + X i β 2) + \frac{\int_{e i (β)}^{\infty} u d F \land e (u, β)}{1 - F \land e (e i (β), β)}, ~ ¯ T (β) = \frac{1}{n} \sum_{i = 1}^{n} T ~ i (β), X ¯ = \frac{1}{n} \sum_{i = 1}^{n} X i Z ¯ = \frac{1}{n} \sum_{i = 1}^{n} Z i, F \land e (t, β) = 1 - \underset{i : e i (β) < t}{Π} [1 - \frac{Δ i}{\sum_{j = 1}^{n} I {e j (β) \geq e i (β)}}]

and

e i (β) = \log (V i) - Z_{i}^{T} β 1 - X i β 2

. The estimating equation

S (β; V, X, Z) = 0

is based upon the normal equations of the least squares method and is then adjusted for censoring (see Buckley and James² for details). The estimating function

S (β; V, X, Z)

involves a non-smooth function

F \land e (\cdot, β)

making the estimating function non-continuous and non-monotone in β. Commonly, in the traditional functional approach of handling covariate measurement errors where unobserved X is treated as an unknown constant, one seeks an estimating function

S * (β; V, W, Z)

such that

E {S * (β; V, W, Z) | V, X, Z} = S (β; V, X, Z)

. However, due to the presence of

F \land e (t, β)

T ~ i (β)

, it is not obvious how to construct such function

S * (β; V, W, Z)

. Alternatively, for this problem with four infinite-dimensional nuisance parameters: (a) the distribution of e; (b) the distribution of the censoring process; (c) the distribution of X given Z; and (d) the distribution of the measurement errors, it would be interesting to investigate the existence and computational feasibility of an efficient estimator along the lines of Ma and Li²⁰ and Ma and Carroll.²¹ The most challenging aspect will be handling the censoring process that may depend on Z. To circumvent these issues we propose a likelihood-based approach with only a few general regularity assumptions on these nuisance distributions, and statistical inferences are made using the MCMC method.

3 Likelihood and priors

In the Bayesian analysis, likelihood function takes a key role. For this purpose we assume that e, X, and U are absolutely continuous random variables and define $S e (\cdot) = 1 - F e (\cdot), f e (e) = dF e (e) / de$ , $f X (x | Z) = d F X (x | Z) / d x$ , and $f U (u) = d F U (u) / d u$ . Then the likelihood of the observed data ignoring the components related to the censoring is

L obs = Π_{i = 1}^{n} \int S_{e}^{1 - Δ i} (\log (V i) - β_{1}^{T} Z i - β 2 X i) f_{e}^{Δ i} (\log (V i 1) - β_{1}^{T} Z i - β 2 X i) f X (X i | Z i) Π_{j = 1}^{m} f U (W ij - X i) d X i

For non-parametric modeling of F_e, F_X(·|Z), F_U, oftentimes a DP mixture model is used that can essentially capture any shape for the distribution of the underlying variable. However, the computation involving a DP prior is time consuming, and it is proportional to the sample size. For efficient computation we shall use a FDDP prior. Before we describe the FDDP, we provide a general definition of the stick-breaking process.

A stick-breaking process is a random probability measure $P$ defined as $P (A) = \sum_{k = 1}^{N} p k δ Y k (A)$ for a measurable set A. Here $δ Y k (\cdot)$ denotes a measure concentrated at Y_k, Y_k are i.i.d. from a distribution H, N is the number of components, and p_k are random probabilities such that $0 \leq p k \leq 1$ and $\sum_{k = 1}^{N} p k = 1$ . Since p_k and Y_k are random, $P (A)$ is also random. The name stick breaking comes due to the structure of random weights p_k, where

p 1 = V 1, p 2 = (1 - V 1) V 2, \dots, p k = (1 - V 1) \dots (1 - V k - 1) V k

and

V k ~ indep Beta (a k, b k)

are assumed to be independent of Y_k. This process allows finite and infinite values of N. As special cases it includes DP, Poisson-DP, Dirichlet-multinomial process, etc.¹³ For a DP, a_k = 1, b_k = α, and N = ∞, and it is denoted by DP(αH). A FDDP, denoted by DP_N(αH), has a finite number of components N and

(p 1, \dots, p N) ~ Dirichlet (α / N, \dots, α / N)

.²² Theorem 2 of Ishwaran and Zarepour²² states that for any real valued measurable integrable function g,

DP N (α H) (g) \to DP (α H) (g)

in distribution as

N \to \infty

. They also described a convenient mechanism of selecting N, the maximum possible cluster size. We shall use N_e, N_u, and N_x to denote N for the FDDPs corresponding to e, U, and X, respectively.

Now we assume that

e i | θ ie ~ indep Normal (θ ie, 1, θ ie, 2), θ ie = (θ ie, 1, θ ie, 2) θ ie | P e ~ i . i . d . P e (arandomdistribution) P e ~ DP N e (α e H 0 e), i = 1, \dots, n

The random probability measure P_e follows a FDDP with the base probability measure H₀ _e on R × R⁺, and H₀ _e is viewed as the center of the random process P_e. The above assumption implies that for any measurable set A, the prior expectation of the probability that

θ e \in A

under the probability measure

DP N e (α e H 0 e)

E {pr P e (θ e \in A)} = pr H 0 e (θ e \in A)

, and

var {pr P e (θ e \in A)} = pr H 0 e (θ e \in A) {1 - pr H 0 e (θ e \in A)} (1 + α e / N e) (1 + α e) - 1

. Thus, larger values of α_e lead to smaller variance of

pr P e (θ e \in A)

. Therefore, α_e can be interpreted as a precision parameter. We assume that under H₀ _e ,

θ e, 2 ~ IG (a σ e, b σ e)

(IG ≡ Inverse Gamma), and conditional on

θ e, 2

θ e, 1 ~ Normal (m e, τ e θ e, 2)

. Using the prior assumption on θ_ie we can now write

f e (e | p 1 e, \dots, p N e e, Y 1 e, \dots, Y N e e) = \sum_{k = 1}^{N e} p ke \frac{1}{\sqrt{2 π θ e, 2}} \exp {- \frac{(e - θ e, 1) 2}{2 θ e, 2}} δ Y ke (θ e)

(p 1 e, \dots, p N e e) ~ Dirichlet (α e / N e, \dots, α e / N e) Y ke ~ i . i . d . H 0 e

In Appendix A1 we give a brief discussion connecting the non-parametric kernel smoothing and DP mixture ideas for density estimation.

We model F_X(·|Z) as a FDDP mixture of normal distributions. In other words, we assume that

X i | Z i, γ 1, θ ix ~ indep Normal (θ ix, 1 + Z_{i}^{T} γ 1, θ ix, 2), θ ix = (θ ix, 1, θ ix, 2) T θ ix ~ i . i . d . P x (arandomdistribution) P x ~ DP N x (α x H 0 x)

where H₀ _x is the base probability measure on R × R⁺ and α_x is the precision parameter. Under H₀ _x , we assume that

θ x, 2 ~ IG (a σ x, b σ x)

, and conditional on θ_x_,2,

θ x, 1 ~ Normal (m x, τ x θ x, 2)

. An alternative statement of the above model is that

f X (x | Z i, γ 1) = \int f X, P x (x | Z i, γ 1) d DP N x (P x; α x H 0 x)

, where

f X, P x (x | Z i, γ 1) = \int (1 / \sqrt{2 π θ ix, 2}) \exp {- (x - Z_{i}^{T} γ 1 - θ ix, 1) / 2 θ ix, 2} dP x (θ ix)

. If we knew the true P_x, say P_x₀, we would model the distribution of X as

f X, P x 0 (\cdot)

without requiring a FDDP prior on P_x. Furthermore,

f X (x | Z, γ 1, p 1 x, \dots, p N x x, Y 1 x, \dots, Y N x x) = \sum_{k = 1}^{N x} p kx \frac{1}{\sqrt{2 π θ x, 2}} \exp {- \frac{(x - θ x, 1 - Z_{i}^{T} γ 1) 2}{2 θ x, 2}} δ Y kx (θ x) (p 1 x, \dots, p N x x) ~ Dirichlet (α x / N x, \dots, α x / N x) Y kx ~ i . i . d . H 0 x

Now we model F_U(·). Although symmetric measurement error is a commonly used assumption,^18,23 for m ≥ 2 one can still identify the distributions of U and X as long as U has mean zero and finite variance.²⁴ Therefore, in our development we use the weaker assumption that the distribution of

U ij = W ij - X i

has mean zero and finite variance and model it as a finite-dimensional centered Dirichlet process (CDP) mixture of a normal kernel.²⁵ For notational convenience we shall use the index l = m(i − 1) + j to identify duplex (i, j) for

i = 1, \dots, n, j = 1, \dots, m

. Thus, we consider

U ij | θ lu ~ indep \frac{1}{\sqrt{2 π θ lu, 2}} \exp {- \frac{(U ij - θ lu, 1) 2}{2 θ lu, 2}}, θ lu = (θ lu, 1, θ lu, 2) T θ lu | P u ~ i . i . d . P u (arandomdistributionwithmeanzero) P u ~ CDP N u (α u H 0 u)

More specifically, the zero mean of U_ij is ensured by the fact that

E (θ lu | P u) = 0

when P_u is randomly drawn from the

CDP N u (α u H 0 u)

. Now we write

f U (W - X | p 1 u, \dots, p N u u, Y 1 u, \dots, Y N u u) = \sum_{k = 1}^{N u} p ku \frac{1}{\sqrt{2 π θ u, 2}} \exp {- \frac{(W - X - θ u, 1) 2}{2 θ u, 2}} δ Y ku (θ u)

(p 1 u, \dots, p N u u) ~ Dirichlet (α u / N u, \dots, α u / N u) Y ku = Y_{ku}^{*} - μ_{u}^{*} μ_{u}^{*} = \sum_{k = 1}^{N u} p ku Y_{ku}^{*} Y_{ku}^{*} ~ i . i . d . H 0 u

Here H₀ _u is the base probability measure on R × R⁺. Under H₀ _u we assume that the second component of

Y_{ku}^{*}

Y_{ku, 2}^{*} ~ IG (a σ u, b σ u)

, and conditional on

Y_{ku, 2}^{*}

, the first component of

Y_{ku}^{*}

Y_{ku, 1}^{*} ~ Normal (m u, τ u Y_{ku, 2}^{*})

. We further assume that a priori

β 1 ~ Normal (μ β 1, Σ β 1)

β 2 ~ Normal (μ β 2, σ_{β 2}^{2})

, and

γ 1 ~ Normal (μ γ 1, Σ γ 1)

. On α_e, α_u, and α_x we put

Gamma (a α e, b α e)

Gamma (a α u,

b α u)

, and

Gamma (a α x, b α x)

priors, respectively. Also, we assume that a priori

τ e ~ IG (g e, h e)

τ u ~ IG (g u, h u)

, and

τ x ~ IG (g x, h x)

. We use

IG (η e, ζ e), IG (η u, ζ u)

, and

IG (η x, ζ x)

priors on

b σ e, b σ u

, and

b σ x

, respectively.

Further notation is needed for posterior computation. Define $Θ_{e}^{T} = (θ 1 e, \dots, θ ne)$ , $Θ_{x}^{T} = (θ 1 x, \dots, θ nx)$ , and $Θ_{u}^{T} = (θ 1 u, \dots, θ Mu)$ , where M = n × m. Let $φ e$ be an N_e × 2 matrix that contains N_e distinct elements of Θ _e . Similarly define $φ x$ and $φ u$ . For updating random elements of Θ _e , define configuration indicators $s_{e}^{T} = (s 1 e, \dots, s ne)$ such that s_ie = j if $θ ie = φ je, j = 1, \dots, N e, i = 1, \dots, n$ . Also define the size of the jth cluster $n_{j}^{e} = \sum_{i = 1}^{n} I (s ie = j)$ , for $j = 1, \dots, N e$ . Thus, $0 \leq n_{j}^{e} \leq n$ and $\sum_{j = 1}^{N e} n_{j}^{e} = n$ . Similarly, define $n_{j}^{x} = \sum_{i = 1}^{n} I (s ix = j)$ that satisfies $0 \leq n_{j}^{x} \leq n$ and $\sum_{j = 1}^{N x} n_{j}^{x} = n$ , and $n_{j}^{u} = \sum_{l = 1}^{M} I (s lu = j)$ with $0 \leq n_{j}^{u} \leq M$ and $\sum_{j = 1}^{N u} n_{j}^{u} = M$ .

Since knowing s_e and $φ e$ is equivalent to knowing Θ _e , in the MCMC method Θ _e is updated via resampling s_e and $φ e$ . Similarly, s_x, s_u can be defined, and Θ _x is updated by resampling s_x and $φ x$ and Θ _u is updated by resampling s_u and $φ u$ . From now on, we shall write θ_ie as $φ_{s ie e}^{T} = (φ s ie e, 1, φ s ie e, 2)$ . Similarly, we shall use $φ s ix x$ and $φ s lu u$ instead of θ_ix and θ_lu.

4 Posterior computation and parameter estimation

Inference regarding the parameters are made from the respective posterior distribution. Using the MCMC method we draw random numbers from the posterior distribution.

Define $T_{i}^{*} = \log (T i)$ . When Δ _i = 0 the value of $T_{i}^{*}$ is unknown. Then it will be treated as an unknown parameter in our Bayesian computation and resampled conditional on the observed data and the other parameters. The important feature of the following MCMC technique is that all of the conditional distributions except the one related to α_e, α_u, and α_x are in the form of standard well-known distributions. We follow Ishwaran and James¹³ for updating the parameters related to the stick-breaking priors.

In the MCMC method we repeat the steps 1–8 (given in Appendix A2) for a large number (e.g. 20,000) of iterations. Along with the unknown parameters and hyperparameters we shall resample all X_i for $i = 1, \dots, n$ , and $T_{i}^{*}$ for those i where Δ _i = 0.

After discarding the first few thousand samples (e.g. 5,000) as burn-in (see, e.g., Cowles and Carlin²⁷), we shall consider the remaining MCMC samples as the random numbers from the joint posterior distribution of the parameters. These sampled observations will be used for calculating parameter estimates and other statistics.

5 Other statistical inferences based on posterior samples

5.1 Estimation of survival probabilities

In addition to the estimation of β in the AFT model (1), another key objective in this context is to estimate the survival probability $pr (T \geq t 0 | X 0, Z 0, Θ)$ for given t₀, X₀, and Z₀, where Θ denotes the set of all parameters. Let $π (Θ | D)$ be the generic notation for the posterior distribution of Θ given the observed data $D$ . Then a random number from the posterior distribution of the survival probability can be obtained by computing $pr (T \geq t 0 | X 0, Z 0, Θ)$ when Θ is randomly drawn from $π (Θ | D)$ . A Bayes estimator of this survival probability is the posterior mean $pr \land (T \geq t 0 | X 0, Z 0, Θ) = \int pr (T \geq t 0 | X 0, Z 0, Θ) π (Θ | D) d Θ$ that can be estimated by taking the Monte Carlo average of

1 - \sum_{k = 1}^{N e} p ke Φ {\frac{\log (t 0) - β_{1}^{T} Z 0 - β 2 X 0 - θ ke, 1}{\sqrt{θ ke, 2}}}

over B (e.g. B = 10,000 or more) MCMC samples of

(β_{1}^{T}, β 2, θ_{ke}^{T}, p ke)

drawn from their joint posterior distribution.

5.2 Model selection

In clinical studies we are also interested in testing hypotheses, such as H₀: β₂ = 0 versus H₁: β₂ ≠ 0. In the Bayesian set-up one can conduct hypothesis testing by calculating the Bayes factor $BF = pr (D | H 1) pr (H 0) / {pr (D | H 0) pr (H 1)}$ , where pr(H₀) and pr(H₁) are the prior probabilities of H₀ and H₁, $pr (D | H k) = \int pr (D | Θ k) π (Θ k)$ dΘ _k for k = 0, 1 with Θ _k being the finite- and infinite-dimensional parameter under the hypothesis H_k, and $π (Θ k)$ is the corresponding prior distribution. Usually BF larger than 10 indicates a strong evidence for the alternative model specified by H₁. Following Newton and Raftery²⁸ we shall calculate the marginal probability or likelihood $pr (D | H k)$ using the harmonic mean:

pr (D | H k) = [E π (Θ k | D) {\frac{1}{L (D | Θ k)}}] - 1

that can be estimated by

[B - 1 \sum_{b = 1}^{B} L - 1 (D | Θ_{k}^{(b)})] - 1,

where

(Θ_{k}^{(1)}, \dots, Θ_{k}^{(B)})

are B MCMC samples from the posterior distribution

π (Θ k | D)

. Since under the Bayesian set-up unobserved X_i is also considered as an unknown parameter, the likelihood is, under H₀,

L (D | Θ 0) = Π_{i = 1}^{n} (\frac{1}{\sqrt{2 π φ s ie e, 2}} \exp [- \frac{{\log (V i) - β_{1}^{T} Z i - φ s ie e, 1} 2}{2 φ s ie e, 2}]) Δ i \times {1 - Φ (\frac{\log (V i) - β_{1}^{T} Z i - φ s ie e, 1}{\sqrt{φ s ie e, 2}})} 1 - Δ i \times Π_{j = 1}^{m} \frac{1}{\sqrt{2 π φ s lu u, 2}} \exp {- \frac{(W ij - X i - φ s lu u, 1) 2}{2 φ s lu u, 2}} \frac{1}{\sqrt{2 π φ s ix x, 2}} \exp {- \frac{(X i - φ s ix x, 1 - Z_{i}^{T} γ 1) 2}{2 φ s ix x, 2}}

and similarly under H₁,

L (D | Θ 1) = Π_{i = 1}^{n} (\frac{1}{\sqrt{2 π φ s ie e, 2}} \exp [- \frac{{\log (V i) - β_{1}^{T} Z i - β 2 X i - φ s ie e, 1} 2}{2 φ s ie e, 2}]) Δ i \times {1 - Φ (\frac{\log (V i) - β_{1}^{T} Z i - β 2 X i - φ s ie e, 1}{\sqrt{φ s ie e, 2}})} 1 - Δ i \times Π_{j = 1}^{m} \frac{1}{\sqrt{2 π φ s lu u, 2}} \exp {- \frac{(W ij - X i - φ s lu u, 1) 2}{2 φ s lu u, 2}} \frac{1}{\sqrt{2 π φ s ix x, 2}} \exp {- \frac{(X i - φ s ix x, 1 - Z_{i}^{T} γ 1) 2}{2 φ s ix x, 2}}

In the real data analysis this Bayes factor approach will also be used for model comparisons where we compute marginal probability of

D

under a given model. One numerical problem in calculating

\sum_{b = 1}^{B} 1 / L (D | Θ_{k}^{(b)})

is that oftentimes

l_{k}^{(b)} = \log {L (D | Θ_{k}^{(b)})}

is a large positive or negative number in the order of 1000, making it impossible to calculate the quantity. Thus, we adopt the following approximation using the Taylor series expansion:

\frac{B}{\sum_{b = 1}^{B} \exp (- l_{k}^{(b)})} \approx \frac{B}{\sum_{b = 1}^{B} {\exp (- μ l k) - (l_{k}^{(b)} - μ l k) \exp (- μ l k) + 0.5 (l_{k}^{(b)} - μ l k) 2 \exp (- μ l k)}} = \exp (μ l k) (1 + 0.5 σ_{*, k}^{2}) - 1

where

μ l k = \sum_{b = 1}^{B} l_{k}^{(b)} / B

and

σ_{*, k}^{2} = \sum_{b = 1}^{B} (l_{k}^{(b)} - μ l k) 2 / B

. Hence, based on this approximation

\log {pr (D | H k)} \approx μ l k - \log (1 + 0.5 σ_{*, k}^{2})

6 Simulation studies

6.1 Simulation design

While a violation of model assumptions may lead to biased estimates of the parameters, the amount of bias depends on the degree of violation, and intricate interplay among the several model assumptions and their violations. We conducted simulation experiments with several scenarios, but due to limited space we shall discuss mainly two scenarios that clearly show the advantage of the proposed method in terms of bias whereas for the other scenarios the semiparametric and partly semiparametric (we shall discuss it in the next paragraph) approaches are comparable. We point out that inconsistency of partly semiparametric methods are manifested via large bias in the parameter estimates. We simulated a cohort of size n = 200 and 300, by simulating $Z ~ Normal (0, 1)$ , and then X and e in the following scenarios. Finally, we obtained T by setting $\log (T) = 1 + Z + 2 X + e$ . Two (m = 2) erroneous measurements W_i₁ and W_i₂ were obtained by adding U_i₁ and U_i₂ with X_i, for $i = 1, \dots, n$ . For scenario 1, $e ~ Exponential (1), X ~ 0.2 Z + (1 / 3) Normal (0, 0.72) + (2 / 3) Normal (2, 0.32)$ , and $U ~ Gamma (1, 1) - 1$ . To create approximately 25% and 50% censored data the censoring variable was simulated as $C = 0.5 X 2 + Unif (0, 2, 000)$ and $C = 0.5 X 2 + Unif (0, 400)$ , respectively. For scenario 2, $e ~ t 3, X ~ {Gamma (6, 0.5) - 3} / 1.22$ , $U ~ Normal (0, 0.712)$ , and C followed two distributions: $C = 0.5 Z 2 + 0.5 X 2 + Unif (0, 40)$ and $C = 0.5 Z 2 + 0.5 X 2 + Unif (0, 5)$ , for 25% and 50% censoring, respectively. For both scenarios we took $var (U) / var (X) \times 100 % = 50 %$ to closely match with the noise-to-signal ratio of the real data. Note that in these scenarios C violates our assumption by making it depend on unobserved X variable. The results when C does not depend on X are similar, thus is omitted. Also, we have intentionally taken non-normal distribution for e, U, and X, to show the robustness of our approach. For completeness, we also ran additional simulations with normal distributions for X, e, and U. The results indicate that SP, SPPE, SPPU, and SPPX worked equally well in this case. The details are omitted.

6.2 Methods for the analyses

The observed data were $(V i, Δ i, Z i, W ij, j = 1, 2, i = 1, \dots, n)$ , and X was no longer used in the analysis stage. The first method is the naive method, where we used $W ¯ i = \sum_{j = 1}^{2} W ij / 2$ in place of X_i in the Buckley–James method and used an existing program to compute the estimates (bj within the R package rms), and this approach will be referred to as the naive method. Next, we analyzed the data using the regression calibration (RC) approach. Here we assume that $W ¯ i | X i ~ Normal (X i, σ_{w | x}^{2})$ and $X i | Z i ~ Normal (γ 0 + γ 1 Z i, σ_{x | z}^{2})$ which imply $X i | W ¯ i, Z i ~ Normal [(σ_{x | z}^{- 2} + σ_{w | x}^{- 2}) - 1 {(γ 0 + γ 1 Z i) σ_{x | z}^{- 2} + W ¯ i σ_{w | x}^{- 2}}, (σ_{x | z}^{- 2} + σ_{w | x}^{- 2}) - 1]$ . We then analyzed the data with X_i being replaced by $X \land i = (1 / {σ \land}_{x | z}^{2} + 1 / {σ \land}_{w | x}^{2}) - 1 {(γ \land 0 + γ \land 1 Z i) / {σ \land}_{x | z}^{2} + W ¯ i / {σ \land}_{w | x}^{2}}$ in the Buckley–James method. Here $γ \land 0$ and $γ \land 1$ are the estimated coefficients obtained by regressing $W ¯ i$ on Z_i, $i = 1, \dots, n, {σ \land}_{w | x}^{2} = (2 n) - 1 \sum_{i = 1}^{n} (W i 1 - W i 2) 2$ and ${σ \land}_{x | z}^{2} = (n - 2) - 1 \sum_{i = 1}^{n} (W ¯ i - γ \land 0 - γ \land 1 Z i) 2 - {σ \land}_{w | x}^{2}$ . Next, we analyzed the data using the proposed method which is referred to as the semiparametric method (SP) where we treated all three infinite-dimensional nuisance parameters non-parametrically.

One may analyze these data sets using several parametric and partly semiparametric approaches. In principle, these approaches may produce biased results when the parametric assumptions are violated. For the sake of comparisons, here we also analyzed the data sets using three partly semiparametric approaches denoted by SPPE, SPPU, SPPX, where two of the three nuisance parameters were treated non-parametrically while the third was treated parametrically. The SPPE model is the same as SP model except that e is modeled parametrically as $e ~ Normal (θ e, 1, θ e, 2), θ e, 2 ~ IG (a σ e, b σ e), θ e, 1 | θ e, 2 ~ Normal (m e, τ e θ e, 2)$ , $τ e ~ IG (g e, h e), b σ e ~ IG (1, 1)$ . The SPPU model is the same as the SP model except that U is modeled parametrically as $U = W - X ~ Normal (0, θ u)$ , $θ u ~ IG (a σ u, b σ u), b σ u ~ IG (1, 1)$ . The SPPX model is the same as the SP except that X given Z is modeled parametrically as $X | Z ~ Normal (θ x, 1 + γ_{1}^{T} Z i, θ x, 2), θ x, 2 ~ IG (a σ x, b σ x)$ , $θ x, 1 | θ x, 2 ~ Normal (m x, τ x θ x, 2)$ , $b σ x ~ IG (1, 1), τ x ~ IG (g x, h x)$ .

Although the general modeling technique and inference method of SP are described in Sections 3 and 4, some necessary details are described in this paragraph. Before each analysis we re-centered the W values by subtracting the sample mean of all of the W from W itself. For the Bayesian methods (SP, SPPE, SPPU, SPPX), posterior inference was made through the MCMC method with 20,000 iterations using the following priors and hyperparameters. We took the RC estimates of β₁ and β₂ as the prior mean of β₁ and β₂, and used five as the prior variance for β₁ and β₂. We used $γ \land 1$ as the prior mean for γ₁ and used 2 times the corresponding standard error as the prior standard deviation of γ₁. We set $a α e = a α u = a α x = b α e = b α u = b α x = 1$ that lead to an Exponential(1) prior for the precision parameters that covers a wide range of plausible values. Then we set $η e = ζ e = η u = ζ u = η x = ζ x = 1, g e = h e = g u = h u = g x = h x = 1$ . By setting these parameters of the inverse gamma distribution to both be one, we allow very large finite variances of the distributions of the hyperparameters that in turn result in estimates that are less affected by the prior choice. In addition, we set $m e = m u = m x = 0$ that are involved in the base probability measure of the FDDPs. Consider the case of $θ e = (θ e, 1, θ e, 2) T$ . Conditional on θ_e_,2 and τ_e, under the base measure of the FDDP, we assumed $θ e, 1 ~ Normal (m e, τ e θ e, 2)$ . For any choice of m_e, this normal distribution can cover a wide spectrum of values for appropriate choice of τ_e and θ_e_,2. Thus, even a completely wrong choice of m_e is compensated by flexible values of τ_e and θ_e_,2 supported by their almost non-informative prior distributions. Thus, following the same analogy for m_u and m_x as well, the choice of m_e, m_u, m_x does not need to be perfect. We initialize $τ e = τ u = τ x = 1, α e = α u = α x = 2, X i = X \land i$ , and set β to the RC estimates of β and γ₁ to $γ \land 1$ .

The “error” in approximating a DP by a FDDP can be measured via the L₁ difference between two marginal probabilities, one corresponding to the FDDP and the other corresponding to the DP. Had we observed $e 1, \dots, e n$ , then based on our model assumption on the distribution of e of the AFT model we could write $f N e (e 1, \dots, e n) \equiv \int {\int Π_{i = 1}^{n} f (e i | θ ie) dP e (θ ie)} d DP N e (P e; α e H 0 e),$ and had we assumed that $P e ~ DP (α e H 0 e)$ , then $f \infty (e 1, \dots, e n) \equiv \int {\int Π_{i = 1}^{n} f (e i θ ie) dP e (θ ie)} d DP (P e; α e H 0 e) .$ Ishwaran and Zarepour²² showed that $\int | f N e (e 1, \dots, e n) - f \infty (e 1, \dots, e n) | de 1 \dots de n \approx 4 n \exp {- (N e - 1) / α e}$ . With N_e = 50, α_e = 1, and n = 1, 036 (the sample size for the data) the error is $2.173 \times 10 - 18$ , and for α_e = 2, it is $9.489 \times 10 - 08$ . Our numerical experience shows that all precision parameters α_e, α_u, and α_x are usually smaller than 1, and as long as the error ( $= 4 n \exp {- (N e - 1) / α e}$ ) is in the order of 10^– ⁵, the results do not vary much with the choice of N_e. We have used the same analogy in choosing N_u and N_x both equal to 50.

6.3 Results

We report the estimated bias (Bias), empirical standard deviations of the estimates (SD), and mean squared error (MSE) based on 500 replications. Tables 1 and 2 contain results for scenarios 1 and 2, respectively. All methods show finite sample bias, and for the naive, RC, and partly semiparametric methods usually bias (see the estimates of β₂) increases with the censoring percentage. However, the naive estimates are much more biased than any other methods. Although the RC method is generally inconsistent, for small values of β₂

(β 2 < 1)

, and relatively small measurement errors, the RC estimates are pretty satisfactory (not presented here). However, in scenarios 1 and 2, RC estimates are quite biased. Considering the bias of all methods in different scenarios, the SP method becomes the most robust approach. In Table 1, other than the naive and RC methods, the largest bias is seen in SPPX and then SPPU. However, SPPE where a normal model is used for e turned out to be comparable with the our proposed model SP. An intuitive explanation is that larger values of T involving larger values of e are likely to get censored more often. The tail probabilities of e, needed for handling censored data, are moderately well-approximated by tail probabilities of a normal distribution when the true distribution of e is Exponential (1). Table 2 clearly shows that if the model assumption regarding e is grossly violated that could affect adversely the parameter estimates, which is evident in large bias in the SPPE method. The second largest bias is shown in SPPX. SPPU performs as good as SP, as the normal distribution assumption on U holds true in this scenario.

Table 1.

Results of the simulation study where $\log (T) = Z + 2 X + 1 + e$ , $e ~ Exponential (1), X ~ 0.2 Z + (1 / 3) Normal (0, 0.72) + (2 / 3) Normal (2, 0.32)$ , $C = 0.5 X 2 + Unif (0, 2, 000)$ for 25% censoring, $C = 0.5 X 2 + Unif (0, 400)$ for 50% censoring, and $U ~ Gamma (1, 1) - 1$ . Here $N e = N u = N x = 50$ .

		Bias	SD	MSE	Bias	SD	MSE
Method	Parameter	n = 200 & 25% censoring			n = 200 & 50% censoring
Naive	β ₁	0.152	0.122	0.038	0.127	0.138	0.035
	β ₂	–0.529	0.102	0.291	0.590	0.122	0.183
RC	β ₁	0.012	0.142	0.020	–0.024	0.162	0.026
	β ₂	0.153	0.209	0.067	0.327	0.245	0.167
SP	β ₁	0.001	0.098	0.009	0.011	0.108	0.012
	β ₂	0.007	0.091	0.008	0.025	0.099	0.011
SPPE	β ₁	–0.001	0.098	0.009	0.011	0.109	0.012
	β ₂	0.009	0.093	0.009	0.027	0.100	0.011
SPPU	β ₁	0.003	0.132	0.017	0.011	0.144	0.021
	β ₂	0.035	0.122	0.016	0.074	0.137	0.024
SPPX	β ₁	0.012	0.104	0.011	0.004	0.117	0.013
	β ₂	0.076	0.117	0.019	0.159	0.128	0.042
		n = 300 & 25% censoring			n = 300 & 50% censoring
Naive	β ₁	0.156	0.095	0.034	0.134	0.112	0.034
	β ₂	–0.529	0.083	0.287	0.409	0.095	0.078
RC	β ₁	0.018	0.106	0.011	–0.019	0.121	0.021
	β ₂	0.140	0.161	0.045	0.316	0.184	0.142
SP	β ₁	0.008	0.076	0.006	0.019	0.089	0.012
	β ₂	0.009	0.068	0.005	0.027	0.075	0.009
SPPE	β ₁	0.004	0.081	0.007	0.019	0.091	0.012
	β ₂	0.012	0.072	0.005	0.028	0.077	0.010
SPPU	β ₁	0.003	0.100	0.010	0.014	0.110	0.016
	β ₂	0.034	0.097	0.010	0.065	0.106	0.022
SPPX	β ₁	0.017	0.084	0.007	0.011	0.095	0.013
	β ₂	0.070	0.088	0.013	0.149	0.095	0.032

Table 2.

Results of the simulation study where $\log (T) = Z + 2 X + 1 + e$ , and $e ~ t 3, X ~ {Gamma (6, 0.5) - 3} / 1.22$ , $C = 0.5 Z 2 + 0.5 X 2 + Unif (0, 40)$ for 25% censoring, $C = 0.5 Z 2 + 0.5 X 2 + Unif (0, 5)$ for 50% censoring, and $U ~ Normal (0, 0.712)$ . Here $N e = N u = N x = 50$ .

		Bias	SD	MSE	Bias	SD	MSE
Method	Parameter	n = 200 & 25% censoring			n = 200 & 50% censoring
Naive	β ₁	–0.008	0.134	0.018	–0.007	0.159	0.025
	β ₂	–0.466	0.142	0.238	–0.488	0.189	0.275
RC	β ₁	–0.008	0.140	0.019	–0.007	0.164	0.027
	β ₂	–0.069	0.197	0.044	–0.097	0.253	0.074
SP	β ₁	0.002	0.133	0.017	0.009	0.152	0.023
	β ₂	0.023	0.210	0.045	0.034	0.271	0.074
SPPE	β ₁	0.009	0.143	0.021	0.033	0.172	0.031
	β ₂	0.085	0.232	0.061	0.129	0.309	0.113
SPPU	β ₁	0.003	0.134	0.018	0.008	0.153	0.023
	β ₂	0.039	0.208	0.045	0.056	0.281	0.082
SPPX	β ₁	–0.001	0.132	0.017	0.005	0.148	0.022
	β ₂	–0.036	0.192	0.038	–0.052	0.238	0.059
		n = 300 & 25% censoring			n = 300 & 50% censoring
Naive	β ₁	–0.009	0.111	0.012	–0.016	0.133	0.018
	β ₂	–0.467	0.111	0.231	–0.496	0.156	0.271
RC	β ₁	–0.010	0.115	0.013	–0.017	0.137	0.019
	β ₂	–0.076	0.150	0.028	–0.116	0.209	0.057
SP	β ₁	–0.002	0.106	0.011	–0.010	0.126	0.016
	β ₂	–0.013	0.139	0.019	–0.022	0.195	0.038
SPPE	β ₁	0.006	0.115	0.013	0.022	0.141	0.023
	β ₂	0.071	0.173	0.035	0.118	0.251	0.076
SPPU	β ₁	–0.001	0.106	0.011	–0.010	0.127	0.016
	β ₂	0.006	0.143	0.020	–0.000	0.201	0.040
SPPX	β ₁	–0.007	0.104	0.011	–0.009	0.125	0.015
	β ₂	–0.072	0.134	0.023	–0.097	0.191	0.045

Table 3.

		Bias	SD	MSE	Bias	SD	MSE
Method	Parameter	n = 200 & 25% censoring			n = 200 & 50% censoring
Naive	β ₁	0.143	0.122	0.035	0.119	0.135	0.032
	β ₂	–0.521	0.109	0.282	–0.400	0.118	0.174
RC	β ₁	0.005	0.145	0.021	–0.036	0.162	0.027
	β ₂	0.152	0.193	0.060	0.329	0.223	0.158
SP	β ₁	0.013	0.097	0.009	0.017	0.104	0.011
	β ₂	–0.016	0.146	0.022	0.012	0.124	0.015
SPPE	β ₁	–0.001	0.108	0.012	0.009	0.113	0.013
	β ₂	0.007	0.132	0.017	0.029	0.112	0.013
SPPU	β ₁	0.016	0.127	0.016	0.019	0.136	0.019
	β ₂	–0.012	0.108	0.012	0.048	0.120	0.017
SPPX	β ₁	0.011	0.110	0.012	–0.002	0.123	0.015
	β ₂	0.077	0.127	0.022	0.177	0.140	0.051

Table 4.

Results of the simulation study where $\log (T) = Z + 2 X + 1 + e$ , $e ~ Normal (0, 1), X ~ 0.2 Z + (1 / 3) Normal (0, 0.72) + (2 / 3) Normal (2, 0.32)$ , $C = 0.5 X 2 + Unif (0, 500)$ for 25% censoring. Here $N e = N u = N x = 50$ .

		Bias	SD	MSE	Bias	SD	MSE
Method	Parameter	$σ_{u}^{2} = 0.5$			$σ_{u}^{2} = 1$
$U ~ Normal (0, σ_{u}^{2})$
Naive	β ₁	0.095	0.118	0.023	0.159	0.135	0.043
	β ₂	–0.291	0.098	0.094	–0.515	0.103	0.277
RC	β ₁	0.014	0.124	0.016	0.019	0.155	0.024
	β ₂	0.109	0.144	0.033	0.184	0.221	0.083
SIMEX1	β ₁	0.045	0.144	0.022	0.122	0.165	0.043
	β ₂	–0.095	0.141	0.029	–0.322	0.153	0.127
SIMEX2	β ₁	0.011	0.122	0.015	0.046	0.144	0.023
	β ₂	–0.021	0.121	0.015	–0.091	0.148	0.030
SP	β ₁	0.016	0.132	0.017	0.030	0.155	0.025
	β ₂	0.020	0.123	0.015	0.011	0.153	0.023
$U ~ σ u {Gamma$ $(1, 1) - 1}$
Naive	β ₁	0.085	0.119	0.022	0.141	0.137	0.039
	β ₂	–0.295	0.100	0.097	–0.517	0.113	0.280
RC	β ₁	0.004	0.125	0.015	0.001	0.155	0.024
	β ₂	0.105	0.136	0.029	0.180	0.209	0.076
SIMEX1	β ₁	0.007	0.143	0.020	0.055	0.169	0.031
	β ₂	–0.092	0.169	0.037	–0.300	0.221	0.139
SIMEX2	β ₁	–0.003	0.122	0.015	0.013	0.144	0.021
	β ₂	–0.027	0.114	0.014	–0.094	0.144	0.029
SP	β ₁	0.004	0.117	0.013	0.009	0.132	0.017
	β ₂	0.022	0.105	0.011	0.015	0.157	0.024
$U ~ σ u {0.5 Normal$ $(0.9, 0.52) + 0.5 Normal (- 0.9, 0.52)}$
Naive	β ₁	0.098	0.119	0.023	0.165	0.137	0.046
	β ₂	–0.305	0.093	0.102	–0.538	0.095	0.299
RC	β ₁	0.013	0.123	0.015	0.018	0.152	0.023
	β ₂	0.116	0.146	0.035	0.195	0.233	0.093
SIMEX1	β ₁	0.052	0.142	0.022	0.132	0.165	0.044
	β ₂	–0.104	0.127	0.027	–0.349	0.129	0.139
SIMEX2	β ₁	0.013	0.123	0.015	0.058	0.143	0.022
	β ₂	–0.022	0.119	0.015	–0.100	0.144	0.031
SP	β ₁	0.008	0.125	0.015	0.015	0.138	0.019
	β ₂	0.032	0.111	0.013	0.037	0.126	0.017

Table 5.

Results for the ACTG AIDS clinical trial data. For the naive Buckley–James method the 95% interval refers to the Wald-type confidence interval whereas for the RC method the 95% interval refers to the percentile interval based on 1000 bootstrap samples. For the Bayesian methods the 95% intervals refer to the equal tail credible intervals. For the Bayesian methods we present the posterior mean of the parameters as the estimates. Here Z, Z + D, Z + Z, and D stand for zidovudine, zidovudine plus didanosine, zidovudine plus zalcitabine, and didanosine, respectively.

Method	Z + D (Ref: Z)	Z + Z (Ref: Z)	D (Ref: Z)	$\log (CD 4)$
Accelerated failure time model
NV	0.333	0.411	0.263	1.001
	(–0.042, 0.708)	(0.013, 0.808)	(–0.101, 0.625)	(0.566, 1.433)
RC	0.332	0.407	0.265	1.217
	(0.044, 0.645)	(0.080, 0.755)	(–0.049, 0.554)	(0.563, 2.094)
SP	0.406	0.512	0.350	1.360
	(0.054, 0.771)	(0.164, 0.889)	(0.019, 0.702)	(0.864, 1.897)
SPPE	0.414	0.514	0.355	1.355
	(0.066, 0.777)	(0.161, 0.897)	(0.023, 0.701)	(0.866, 1.898)
SPPU	0.410	0.507	0.352	1.417
	(0.059, 0.776)	(0.139, 0.891)	(0.016, 0.693)	(0.902, 1.985)
SPPX	0.408	0.513	0.356	1.394
	(0.071, 0.780)	(0.156, 0.888)	(0.020, 0.700)	(0.879, 1.952)
Piecewise exponential model
NVPE	–0.801	–1.010	–0.778	–1.927
	(–1.365, –0.262)	(–1.617, –0.443)	(–1.326, –0.237)	(–2.575, –1.284)
RCPE	–0.801	–1.005	–0.783	–2.32
	(–1.372, –0.267)	(–1.607, –0.432)	(–1.321, –0.246)	(–3.093, –1.542)
PCPE	–0.809	–1.014	–0.795	–2.528
	(–1.388, –0.264)	(–1.629, –0.430)	(–1.334, –0.252)	(–3.413, –1.651)

Table 6.

Results for the ACTG AIDS clinical trial data. Here Z, Z + D, Z + Z, and D stand for zidovudine, zidovudine plus didanosine, zidovudine plus zalcitabine, and didanosine, respectively, and AFT stands for accelerated failure time. The 95% Wald-type confidence intervals are given in parentheses right beneath the estimates. The bootstrap method was used to compute the standard error of the regression calibration and the SIMEX methods.

Method	Z + D (Ref: Z)	Z + Z (Ref: Z)	D (Ref: Z)	$\log (CD 4)$
Semiparametric AFT model, e is nonparametric
SIMEX1	0.322	0.417	0.252	1.092
	(–0.044, 0.708)	(–0.004, 0.838)	(–0.122, 0.626)	(0.384, 1.799)
SIMEX2	0.306	0.396	0.248	1.133
	(–0.058, 0.671)	(–0.025, 0.817)	(–0.141, 0.636)	(0.484, 1.782)
Parametric AFT model, e is Generalized Gamma
NV	0.339	0.443	0.338	0.972
	(0.052, 0.625)	(0.135, 0.750)	(0.064, 0.612)	(0.594, 1.350)
RC	0.339	0.440	0.341	1.181
	(0.053, 0.625)	(0.132, 0.748)	(0.063, 0.619)	(0.740, 1.622)
SIMEX1	0.349	0.435	0.334	1.094
	(0.026, 0.672)	(0.139, 0.731)	(0.024, 0.643)	(0.578, 1.609)
SIMEX2	0.333	0.435	0.332	1.121
	(0.015, 0.656)	(0.141, 0.729)	(0.026, 0.637)	(0.623, 1.618)

Of course, as a price for the robustness, the SD of the SP method is often slightly larger than the competing methods but the MSEs are relatively comparable. The SD of the estimates decreases with sample size, and it increases with the percentage of censoring. The difference between the SP and other partly semiparametric approaches is not as large as the difference between the SP and RC methods. The reason lies in the fact that in other partly semiparametric methods two of the three nuisance parameters are non-parametrically modeled that reduces the degree of model violations. Of course, the bias reduction achieved in the SP compared with other partly semiparametric approaches indicates the supremacy of flexible modeling of all three nuisance parameters.

Prompted by a referee’s comment, we re-ran the simulation study for scenario 1 with $N e = N u = N x = 100$ . The results are presented in Table 3. A close comparison between Tables 1 and 3 reveals that there are no qualitative differences in these results.

In addition to the simulation studies described above, we conducted a small-scale simulation to compare the performance of the proposed method with the SIMEX approach. Here we took the semiparametric AFT model where e was left unspecified. In the first SIMEX approach (refer to as SIMEX1) we considered one of the two measurements of W as the erroneous measurement and estimated the measurement error variance with $\sum_{i = 1}^{n} (W i 1 - W i 2) 2 / (2 n)$ . SIMEX1 does not use all the available data. It is thus likely to lead to more bias. In the second SIMEX (refer to as SIMEX2) we used $W ¯ i = (W i 1 + W i 2) / 2$ as the erroneous measurement and estimated the measurement error variance with $\sum_{i = 1}^{n} (W i 1 - W i 2) 2 / (4 n)$ . We used symmetric and asymmetric measurement error distributions, and used two different values for the measurement error variance var(U) = 0.5 and var(U) = 1. For comparisons we also present the naive and the RC along with our SP method. The results are given in Table 4. For obvious reasons SIMEX1 is worse than SIMEX2. When var(U) = 0.5 the performance of SIMEX2 and SP are similar. However, for var(U) = 1 the bias in SIMEX2 is much larger than SP. Large finite sample bias is a reflection of possible inconsistency of SIMEX. Although the bias of the SIMEX estimates seem to be not much affected by the non-normal measurement error, SIMEX is shown to be consistent only in a handful of cases with normal measurement errors and correct extrapolating function.

It is seen that while SP greatly reduces the bias, it is also accompanied by larger variance compared with the naive or the RC approach. This phenomenon is expected since any method that takes into account measurement errors in a covariate would generally result in larger uncertainty in the parameter estimators than the methods that fail to consider the measurement error issue in the analysis. Also, the variances of the estimators increase with the measurement error variance (see Table 4), and so are the MSEs. In our simulation, the bias of the inconsistent methods (such the naive and RC) is overwhelmingly larger than their corresponding variances, resulting in larger MSEs for the naive and RC methods compared with SP. Of course, there is no guarantee that SP would always have smaller MSE than the inconsistent approaches; eventually it all depends on whether the bias over weights the estimation variance (uncertainty) that in turn depends on the measurement error variance. Lastly, we would like to point out that although MSEs are presented in all of the tables, MSE may not be a good measure to compare consistent and inconsistent estimators.

7 Analysis of the data from an AIDS clinical trial study

This data set comes from a randomized, double-blind trial on AIDS known as ACTG 175 study. One of the study aims is to understand the effect of several antiretoviral drugs on human immunodeficiency virus-1 (HIV-1) infected people who had no history of an AIDS-defining illness other than minimal mucocutaneous Kaposi’s sarcoma (see Hammert et al.²⁹ for details).

Subjects were randomly assigned to one of the four therapies, 600 mg of zidovudine, 600 mg of zidovudine plus 400 mg of didanosine, 600 mg of zidovudine plus 2.25 mg of zalcitabine, and 400 mg of didanosine. In our analysis, the event is the development of AIDS or death, and T is defined as the time (in days) from the start of the treatment to the occurrence of the event. According to the ACTG 175 study, for this group of patients, AIDS and death were considered as the primary end points as both were related to at least 50% decline in CD4 counts.²⁹

For our analysis we considered only n = 1, 036 subjects who did not have antiretroviral treatment before this trial. Among them 262, 257, 260, 257 subjects received the above 4 treatments, respectively. These subjects had two (i.e. m = 2) baseline measurements of CD4 counts prior to the start of their treatment. Among the 1036 subjects 85 experienced the above event and the median and average follow-up time were approximately 27 and 32 months, respectively.

We fit model (1) to this data set, where the logarithm of the actual CD4 count at the baseline minus 5.89 is considered as X. The choice of 5.89 is to make the distribution centered around 0. Note that the exact CD4 count in the blood is impossible to measure mainly due to constant movement of these cells between blood and tissues. In addition, within a short time span (a few days) small changes may occur in the CD4 count due to physical activity, stress, good night’s sleep, etc. Therefore, the two baseline measurements are considered to be two erroneous measurements W_i₁ and W_i₂ for X_i, $i = 1, \dots, 1036$ . The estimated noise-to-signal ratio is $var (U) \land / var (X) \land \times 100 % = 42 %$ .

The three dummy variables corresponding to the four treatments are considered to be error-free covariates Z, where 600 mg of zidovudine was considered as the reference category. We analyzed the data using NV, RC, SP, SPPE, SPPU, and SPPX, and the results are presented in Table 5. For SP we used $N e = N u = N x = 50$ . For RC, the 95% confidence intervals were calculated based on 1000 bootstrap samples.

Based on the 95% confidence intervals and credible intervals all methods indicate that $\log (CD 4)$ has a statistically significant effect on the time-to-event. More importantly, after adjusting for the measurement errors, the estimate of the coefficient for CD4 counts, β₂, is quite different in the SP method from the naive estimate. Clearly, for β₂, the naive estimate is closer to zero than the estimates from the other methods, a trend that is also observed in the simulation studies. The results of SP also indicate that compared with zidovudine, the other three therapies have statistically significant effect on delaying the time-to-event. This result is consistent with the findings of Hammer et al.²⁹ Figure 1 shows the estimated survival probabilities and 95% pointwise credible intervals based on SP for each treatment group when CD4 counts were 232 and 539, the approximate 10th and 90th quantiles of baseline CD4 measurements of the subjects. We found that between the competing hypotheses, H₀: β₂ = 0 versus H₁: $β 2 \neq 0$ , the data unequivocally support H₁ as the Bayes factor was much larger than 10. Note that the analysis of this data set by SP with 60,000 MCMC iterations took approximately 3 minutes on a 2.8 GHz Intel Xenon X5560 processor.

We also analyzed the data using a piecewise exponential (PE) model, i.e. we assumed that the hazard of the time-to-event is $λ (t | X, Z) = λ s \exp (β_{1}^{T} Z + β 2 X)$ for $t s - 1 \leq t < t s$ , $s = 1, \dots, q$ , where the time axis is partitioned into $[t 0, t 1), \dots, [t q - 2, t q - 1)$ , $[t q - 1, t q)$ , with t₀ = 0 and t_q = ∞. First we carried out a naive analysis (NVPE) where we fit the PE model by replacing X_i by $W ¯ i$ . Then we conducted a RC analysis (RCPE) by fitting the PE model where X_i was replaced by $X \land i$ defined in Section 6.2. Third, we fit the PE model with a parametric correction for the measurement errors (PCPE) where we assumed that X given Z followed a $Normal (γ 0 + γ_{1}^{T} Z, σ_{x}^{2})$ and $U ~ Normal (0, σ 2)$ . The likelihood of the data is

L = Π_{i = 1}^{n} \int {\sum_{s = 1}^{q} λ s I (t s - 1 \leq V i < t s) \exp (β_{1}^{T} Z i + β 2 X i)}} Δ i \exp {- \sum_{s = 1}^{q} λ s I (t s - 1 \leq V i) \exp (β_{1}^{T} Z i + β 2 X i)} \times \frac{1}{\sqrt{2 π σ_{x}^{2}}} \exp {- \frac{(X i - γ 0 - γ_{1}^{T} Z i) 2}{2 σ_{x}^{2}}} \times \frac{1}{(2 π σ_{u}^{2}) m / 2} \exp {- \sum_{j = 1}^{m} \frac{(W ij - X i) 2}{2 σ_{u}^{2}}} d X i

In all three methods the parameters were estimated in a Bayesian framework using the MCMC method. Here we assumed a priori

λ s ~ Gamma (a λ, b λ), β 1 ~ Normal (μ β 1, Σ β 1)

β 2 ~ Normal (μ β 2, σ_{β 2}^{2})

γ 0 ~ Normal (μ γ 0, σ_{γ 0}^{2})

γ 1 ~ Normal (μ γ 1, Σ γ 1)

σ_{x}^{2} ~ IG (a σ x, b σ x)

, and

σ_{u}^{2} ~ IG (a σ u, b σ u)

. The MCMC details are given in Appendix A3. In particular, we took q = 5 and partitioned the time axis as

[0, t 0.2), [t 0.2, t 0.4), [t 0.4, t 0.6), [t 0.6, t 0.8), [t 0.8, \infty)

, where t_r denotes the rth quantile of the observed failure times. The prior means of β were the RC estimates of the CPH model and 5 was used as the prior variance for all parameters. We set

μ γ 0 = γ \land 0

μ γ 1 = γ \land 1

, where

γ \land 0

and

γ \land 1

denote the estimated intercept and partial slopes for the linear regression of

W ¯ i

on Z_i. Two times the square of the corresponding standard errors were used as the prior variance. Finally, we used

a σ x = a σ u = b σ x = b σ u = 1

. The results are given in the lower panel of Table 5. These results, like the SP analysis, also indicate that the treatments are significantly associated (based on the 95% credible interval) with the time-to-event. In particular, compared with zidovudine, the other treatments reduce the hazard of the event. In addition, log(CD4) is significantly negatively associated with the hazard of the time-to-event, a consistent finding with the SP analysis.

Since SP, SPPE, SPPU, SPPX, and PCPE are all Bayesian methods, we were able to compare these approaches using marginal probabilities $pr (D | M)$ where $M$ stands for a generic model. Figure 2 shows the boxplot of $l_{M}^{(b)}$ based on the MCMC samples after discarding the first 10,000 burn-in samples. The range of its values clearly indicates that a straightforward estimate of $pr (D | M)$ using the harmonic mean is not possible. Therefore, we adopted the approximation given in Section 5.2 and obtained $\log {pr (D | M)}$ as 3306.63, 3292.64, 3022.74, 2937.14, and –171.93 for SP, SPPE, SPPU, SPPX, and PCPE, respectively. The above results along with equal prior probability for each of the model in Bayes factor calculations indicate that the SP model is the best and is closely followed by the SPPE model. This explains why the β₂ estimates in these two methods are so close.

8 Conclusions

In this paper we proposed a non-parametric Bayesian method for fitting the AFT model to a right censored data when a covariate is measured with error. We believe that our approach is the first attempt to solve this problem in a non-parametric framework. While we non-parametrically treat the three components (the stochastic noise of the model for the time-to-event, the distribution of the latent unobserved true covariate, and the distribution of the zero mean measurement errors), the computation is simple and easily programmable due to the novel application of the stick-breaking priors. Due to non-parametric modeling of all nuisance distributions the proposed method outperforms the naive, RC, SIMEX, and other partly semiparametric methods in our simulation studies.

In this work, for notational simplicity we have assumed the number of the replications of the surrogate for the true covariate to be the same across all of the subjects (i.e. m_i = m). This assumption can be relaxed by using some more intense notations with some general regularity conditions on m_i. Furthermore, in principle, the proposed method can be extended to handle interval censored data.³⁰ For this purpose, in the posterior computations one should generate $T_{i}^{*}$ from a normal distribution which is truncated on both sides. In this work we used the Monte Carlo estimates for estimating survival probabilities. We believe that using the importance sampling method with proper importance weights one may improve the efficiency of the estimator. We have focused on time-invariant covariate. However, the measurement error issue may arise in a time-varying covariate; see, for example, Veronesi et al.³¹ who considered the RC and SIMEX methods for handling such a covariate in a Cox regression model. Another interesting paper in this area is by Crowther et al.³² who considered joint modeling of survival outcome using a semiparametric Cox model and longitudinally measured prognostic biomarkers using a linear mixed model. It is worth investigating how the nonparametric approaches considered in this paper could be implemented in these settings. The computational code can be obtained from the authors upon request. We are also creating an R-package for practitioners to use that will be available through our website.

Finally, we briefly discuss the efficiency property of the estimator. A semiparametrically efficient estimator is defined in the class of regular asymptotic linear (RAL) estimators that achieves the efficiency bound (Tsiatis,³³ p. 27). The bound is defined as the supremum of the most efficient RAL estimators of the parametric submodels that is a subset of the semiparametric class of models, and the true model belongs to the parametric submodels. For a parametric model where standard regularity conditions hold, the MLE produces the most efficient RAL estimator. By construction our estimator is not a RAL estimator. Therefore, it is generally difficult to compare it with the corresponding efficiency bound. However, like the Cramer–Rao lower bound, there exists a Bayesian Cramer–Rao lower bound.³⁴ It is worth exploring how a Bayesian minimax estimator can be constructed with that lower bound.

Footnotes

Acknowledgment

The authors gratefully acknowledge the constructive comments from the referees that led to a significant improvement of the manuscript.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was partially supported by the NIH (grant number R03CA176760).

Appendix

References

Carroll

. Measurement error in epidemiologic studies. Encyclopedia of Biostatistics, New York: John Wiley & Sons, 1997.

Buckley

James

. Linear regression with censored data. Biometrika 1979; 66: 429–436.

Lai

Ying

. Large sample theory of a modified Buckley–James estimator for regression analysis with censored data. Ann Statist 1991; 19: 1370–1402.

Zeng

Lin

. Efficient estimation for the accelerated failure time model. J Am Statist Assoc 2007; 102: 1387–1396.

Jin

Lin

Wei

Ying

. Rank-based inference for the accelerated failure time model. Biometrika 2003; 90: 341–353.

Zhou

. Empirical likelihood analysis of the Buckley–James estimator. J Multivar Anal 2008; 99: 649–664.

Cox

Chu

Schneider

Muñoz

. Parametric survival analysis and taxonomy of hazard functions for the generalized gamma distribution. Statist Med 2007; 26: 4252–4374.

Christensen

Johnson

. Modelling accelerated failure time with a Dirichlet process. Biometrika 1988; 75: 693–704.

Kuo

Mallick

. Bayesian semiparametric inference for the accelerated failure-time model. Can J Statist 1997; 25: 457–472.

10.

Walker

Mallick

. A Bayesian semiparametric accelerated failure time model. Biometrics 1999; 55: 477–483.

11.

Xiong

. Accelerated failure time models with covariates subject to measurement Error. Statist Med 2007; 26: 4817–4832.

12.

Yin

. Censored quantile regression with covariate measurement errors. Statistica Sinica 2011; 21: 949–971.

13.

Ishwaran

James

. Gibbs sampling methods for stick breaking priors. J Am Statist Assoc 2001; 96: 161–173.

14.

Ghosh

Ramamoorthi

. Bayesian Nonparametrics, New York: Springer, 2003.

15.

Müller

Roeder

. A Bayesian semiparametric model for case–control studies with errors in variables. Biometrika 1997; 84: 523–537.

16.

Gustafson

Vallee

. A Bayesian approach to case–control studies with errors in covariables. Biostatistics 2002; 3: 229–243.

17.

Sinha

Mallick

Kipnis

Carroll

. Semiparametric Bayesian analysis of nutritional epidemiology data in the presence of measurement error. Biometrics 2010; 66: 444–454.

18.

Delaigle

Hall

Meister

. On deconvolution with repeated measurements. Ann Statist 2008; 36: 665–685.

19.

Carroll

Ruppert

Stefanski

Crainiceanu

. Measurement Error in Nonlinear Models: A Modern Perspective, 2nd edn. Boca Raton, FL: Chapman and Hall/CRC Press, 2006.

20.

. Variable selection in measurement error models. Bernoulli 2010; 16: 274–300.

21.

Carroll

. Locally efficient estimators for semiparametric models with measurement error. J Am Statist Assoc 2006; 101: 1465–1474.

22.

Ishwaran

Zarepour

. Exact and approximate sum representations for the Dirichlet process. Can J Statist 2002; 30: 269–283.

23.

Hall

. Measurement error models with unknown error structure. J R Statist Soc Ser B 2007; 69: 429–446.

24.

Vuong

. Nonparametric estimation of the measurement error model using multiple indicators. J Multivar Anal 1998; 65: 139–165.

25.

Yang

Dunson

Baird

. Semiparametric Bayes hierarchical models with mean and variance constraints. Computat Statist Data Anal 2010; 54: 2172–2186.

26.

Escobar

West

. Bayesian density estimation and inference using mixtures. J Am Statist Assoc 1995; 90: 577–588.

27.

Cowles

Carlin

. Markov chain monte carlo convergence diagnostics: A comparative review. J Am Statist Assoc 1996; 91: 883–904.

28.

Newton

Raftery

. Approximate Bayesian inference with the weighted likelihood bootstrap. J R Statist Soc Ser B 1994; 56: 3–26.

29.

Hammer

Katzenstein

Hughes

. A trial comparing nucleoside monotherapy with combination therapy in HIV-infected adults with CD4 cell counts from 200 to 500 per cubic millimeter. N Engl J Med 1996; 335: 1081–1090.

30.

Huang J and Wellner J. Interval censored survival data: a review of recent progress. In Proceedings of the first Seattle symposium in biostatistics: survival analysis (Lecture Notes in Statistics, vol. 123). Berlin: Springer, 1997; 123–169.

31.

Veronesi

Ferrario

Chambless

. Comparing measurement error correction methods for rate-of-change exposure variables in survival analysis. Statist Meth Med Res 2011; 22: 583–597.

32.

Crowther

Lambert

Abrams

. Adjusting for measurement error in baseline prognostic biomarkers included in a time-to-event analysis: a joint modelling approach. BMC Med Res Methodol 2013; 13: 146–146.

33.

Tsiatis

. Semiparametric Theory and Missing Data, New York: Springer, 2006.

34.

Gill

Levit

. Applications of the van Trees inequality: a Bayesian Cramér–Rao bound. Bernoulli 1995; 1: 59–79.

35.

Gelfand

Hills

Racine-Poon

Smith

. Illustration of Bayesian inference in normal data models using Gibbs sampling. J Am Statist Assoc 1990; 85: 972–985.