Two-phase analysis and study design for survival models with error-prone exposures

Let $T$ be a discrete random variable. Denote the $j$-th discrete value of $T$ by $t_{j}$ and write $\lambda_{0j}=\lambda_{0}(t_{j})$, where $\lambda_{0}$ is the baseline hazard function for $T$ defined by $\lambda_{0}(t_{j})=P(T=t_{j}|T\geq t_{j})$ for $j\in J$, where $J$ is the index set for the discrete times with positive mass. We assume the time to event response $T$ is associated with a $d$-dimensional covariate vector $\mathbf{X}=(X_{1},\ldots,X_{d})^{\top}$ such that the conditional hazard function $\lambda(t|\mathbf{x})=P(T=t|T\geq t,\,\mathbf{X}=\mathbf{x})$ is given by

for some coefficient vector $\boldsymbol{\beta}=(\beta_{1},\ldots,\beta_{d})^{\top}$, where $g:[0,1]\rightarrow\mathbf{R}$ is a monotone transformation and $\mathcal{T}=\{t_{j}:j\in J\}$. For example, the odds transformation $g_{1}(u)=\frac{u}{1-u}$ yields the logit hazard model, where $\textrm{logit}\big{(}\lambda_{j}(\mathbf{x})\big{)}=\alpha_{j}+\boldsymbol{\beta}^{\top}\boldsymbol{x}$ for $\alpha_{j}=\textrm{logit}(\lambda_{0j})$ and $\lambda_{j}(\mathbf{x})=\lambda(t_{j}|\mathbf{x})$. It follows that the likelihood function under the logit hazard model has the same representation with the logistic regression model such that the number of events at a time $t_{j}$ represents binomial outcomes with probability $\lambda_{j}(\mathbf{x})$ for each $j\in J$ (Meier et al., 2003). The complementary log transformation $g_{2}(u)=-\log(1-u)$ gives a proportional log-survival model $\log\big{(}S_{j}(\mathbf{x})\big{)}=\log(S_{0j})\exp(\boldsymbol{\beta}^{\top}\mathbf{x})$, where $S_{j}(\mathbf{x})=P(T\geq t_{j}|\mathbf{X}=\mathbf{x})$ and $S_{0j}=S_{0}(t_{j})$ for the baseline survival function $S_{0}$ for $T$. Kalbfleisch and Prentice (2011) provide further details.

We consider that $T$ may be subject to random censoring prior to the finite maximum follow-up time of $\tau<\infty$. For a random censoring time $C$, independent of $T$, let $Y=\min\{T,C\}$ be the observed censored survival time and $\Delta=\mathbb{I}(T\leq C)$ be the event indicator. We assume that $\mathcal{X}_{N}=\{(Y_{i},\Delta_{i},\mathbf{X}_{i}):1\leq i\leq N\}$ are independent and identically distributed and will not be completely observed on all subjects as $N$-random copies of $(Y,\Delta,\mathbf{X})$. Instead, $\mathcal{X}_{\textrm{I},N}=\{(Y_{i},\Delta_{i},\mathbf{Z}_{i}):1\leq i\leq N\}$ are the data available for all subjects in the first phase of a study, where $\mathbf{Z}=(Z_{1},\ldots,Z_{q})^{\top}$ are discrete surrogates or auxiliary variables associated with $\mathbf{X}$. Complete data for the likelihood are observed on the phase two sample, denoted by $\mathcal{X}_{\textrm{II},n}=\{(Y_{i},\Delta_{i},\mathbf{X}_{i}):i\in\mathcal{I}\}$, where $\mathcal{I}$ is a subset of $\{1,\ldots,N\}$ having the cardinality of $n$. For example, the auxiliary variables $\mathbf{Z}$ might be error-prone discrete measures of $\mathbf{X}$ while $\mathbf{X}$ may or may not be discrete. We also allow for settings where some components of $\mathbf{X}$ are available on all subjects, such as sex or other demographic information ascertained in the first phase of the study. For this setting, we may introduce a slight abuse of notation writing $\mathbf{X}=(\mathbf{X}^{C},\mathbf{X}^{I})$ and $\mathbf{Z}=(\mathbf{X}^{C},\mathbf{A})$, where $\mathbf{X}^{C}$ are the components of $\mathbf{X}_{i}$ observed on everyone, $\mathbf{X}^{I}$ are the incomplete components of $\mathbf{X}$ not observed at the phase one, and $\mathbf{A}$ indicates a generic notation for auxiliary variables in the first phase of the study. In this case, $\mathcal{I}$ denotes the set of subject indices for which the true covariates $\mathbf{X}$ are fully observed, with $\mathbf{X}^{I}$ sampled in the second phase of the study, and the auxiliary phase one variables are limited to $\mathbf{A}$ (Reilly and Pepe, 1995; McIsaac and Cook, 2015).

Unlike previous literature (Lawless et al., 1999; Chatterjee et al., 2003; Chatterjee and Chen, 2007; Tao et al., 2017), we note that the model (2.1) does not include auxiliary variables as predictors, but rather $\mathbf{Z}$ is a surrogate by the Prentice criterion (Prentice, 1989). That means the likelihood function $L_{1}(\boldsymbol{\theta};Y,\Delta,\mathbf{X},\mathbf{Z})$ equals $L_{1}(\boldsymbol{\theta};Y,\Delta,\mathbf{X})$, where $\boldsymbol{\theta}$ is the collection of parameters involved in (2.1). Thus, the complete information for $\boldsymbol{\theta}$ is carried by $(Y,\Delta,\mathbf{X})$, while $\mathrm{E}\{\log L_{1}(\boldsymbol{\theta};Y,\Delta,\mathbf{X})\,|\,\allowbreak Y,\Delta,\mathbf{Z}\}$ may represent extra information of $(Y,\Delta,\mathbf{Z})$ compared to likelihood-based inference using only complete case data for $(Y,\Delta,\mathbf{X})$. In this paper, we consider the mean score method, which is valid when the validation subset $\mathcal{I}$ for the phase two sample is randomly chosen from the full cohort given information obtained from the first phase of a study. Therefore, we use the log-likelihood from all available observations written by

where $h(\mathbf{x}|\mathbf{z})$ denotes the conditional density function of $\mathbf{X}$ given $\mathbf{Z}=\mathbf{z}$ and $\mathcal{I}^{c}=\{1,\ldots,N\}\setminus\mathcal{I}$ indicates a set of $(N-n)$ indices for individuals whose complete covariates are not available. For an individual $i$ and time $j$, we define the observed event indicator ${D}_{ij}=\mathbb{I}(Y_{i}=t_{j},\,\Delta_{i}=1)$ and denote subject $i$’s censored survival time index by $J(i)=\arg\{j\in J:Y_{i}=t_{j}\}$ for $1\leq i\leq n$. Thus, $D_{ik}=0$ for all $k<J(i)$. Then, it follows that the log-likelihood function given $(Y_{i},\Delta_{i},\mathbf{X}_{i})$ can be written by

for each $1\leq i\leq n$. In the above equation, we used the fact that $L_{1}(\boldsymbol{\theta};Y,\Delta,\mathbf{X})=S(Y|\mathbf{X})\lambda(Y|\mathbf{X})^{\Delta}(1-\lambda(Y|\mathbf{X}))^{1-\Delta}$ and the conditional survival function $S_{j}(\mathbf{x})=P(T\geq t_{j}|\mathbf{X}=\mathbf{x})$ was calculated by $\prod_{k=1}^{j-1}(1-\lambda_{k}(\mathbf{x}))$ for $j\geq 2$, together with $S_{1}(\mathbf{x})=1$ by definition.

We apply the mean score method (Reilly and Pepe, 1995; Reilly, 1996) to the conditional hazard model (2.1) when auxiliary variables are discrete. Employing the Expectation-Maximization (EM) technique (Dempster et al., 1977) with (2.2) and (2.3), we may find the maximum likelihood estimator of $\boldsymbol{\theta}$. However, $h(\mathbf{x}|\mathbf{z})$ is generally unknown and a parametric approach for the estimation of $h(\mathbf{x}|\mathbf{z})$ may result in inconsistent inference of $\boldsymbol{\theta}$ in likelihood-based methods. Lawless et al. (1999) introduced a semi-parametric method, estimating the conditional density function $h(\mathbf{x}|\mathbf{z})$ nonparametrically, such that the integration in (2.2) is replaced with a single-step approximation (McIsaac and Cook, 2014). Following the mean score approach, for those not in the phase two subset, we can replace the unobserved score contribution based on $\mathbf{X}$ with its expected value based on the observed phase one data. By replacing the expected value to an empirical mean, the semi-parametric estimation of $\boldsymbol{\theta}$ can then be achieved by maximizing an inverse probability weighted log-likelihood

where $\hat{\pi}(Y_{i},\Delta_{i},\mathbf{Z}_{i})$ is an empirical estimate of ${\pi}(Y_{i},\Delta_{i},\mathbf{Z}_{i})$, the sampling probability of the $i$-th individual selected into the validation subset, which can be consistently estimated by $n(Y_{i},\Delta_{i},\mathbf{Z}_{i})/N(Y_{i},\Delta_{i},\mathbf{Z}_{i})$ as $n$ and $N$ increase. Here, $n(y,\delta,\mathbf{z})$ is the number of subjects in $\mathcal{I}$ who have the same observations with $(y,\delta,\mathbf{z})$ in the first phase study; $N(y,\delta,\mathbf{z})$ is defined similarly to $n(y,\delta,\mathbf{z})$ with replacement of the index set $\mathcal{I}$ with $\{1,\ldots,N\}$ of the full sample.

We note that Flanders and Greenland (1991) also proposed a similar estimating equation to (2.4) based on the pseudo-likelihood method, which coincides to the filling in method introduced by Vach and Blettner (1991) when $\mathbf{X}$ is categorical. However, depending on the choice of transformation $g$ in (2.1), different forms of score equations follow from the above weighted log-likelihood (2.4). In Supplementary Material Section A.2, we provide the detailed forms of the mean score equations and the associated Hessian matrices when the logit transformation $g_{1}$ and the complementary log transformation $g_{2}$ are used.

Our expectation that the optimal design for our discrete time proportional hazards model will also be advantageous for the continuous time Cox model is based on the connection between the parameters in these two models. For this, we briefly review this connection. For further discussion, see Kalbfleisch and Prentice (2011).

From the log-likelihood (2.3), we note that the logit hazard model is the canonical form of the discrete-time survival model (2.1) under the odds transformation $g_{1}(u)=\frac{u}{1-u}$, such that $\textrm{logit}(\lambda_{j}(\mathbf{x}))=\alpha_{j}+\boldsymbol{\beta}^{\top}\mathbf{x}$, where $\alpha_{j}=\textrm{logit}(\lambda_{0j})$ is the logit transformation of the baseline hazard. However, we also note that any model with an arbitrary monotone transformation in (2.1) leads to a reparameterization of the logit hazard model. For example, suppose the complementary log transformation $g_{2}(u)=-\log(1-u)$ defines the true survival model. Then it can be easily seen that $\textrm{logit}(\lambda_{j}(\mathbf{x}))=\exp\big{(}e^{\alpha_{j}+\boldsymbol{\beta}^{\top}\mathbf{x}}\big{)}-1$, where $\alpha_{j}=\log(-\log(1-\lambda_{0j}))$ is the complementary log-log transformation of the baseline hazard, and the application of chain rules to (2.3) is followed by likelihood-based estimation of $\boldsymbol{\theta}=(\boldsymbol{\alpha},\boldsymbol{\beta})$, which is also equivalent to reparameterization of the logit hazard model.

In particular, if usual continuous-time Cox proportional hazards are grouped into discrete disjoint intervals, this will lead to the conditional hazard model (2.1) for discrete-time outcomes with the complementary log transformation. To be specific, let $\lambda^{C}(t|\mathbf{x})$ be the conditional hazard function in the Cox model for the continuous-time survival outcome such that $\lambda^{C}(t|\mathbf{x})=\lambda_{0}^{C}(t)\exp(\boldsymbol{\beta}^{\top}\mathbf{x})$, where $\lambda_{0}^{C}$ is the associated baseline hazard function. Suppose we are only able to observe survival outcomes censored by disjoint intervals $(t_{j-1},t_{j}]$ with $t_{0}\equiv 0$, so that statistical inference on conditional hazards by each interval is a goal of the study. Then, we may denote the conditional hazards on the $j$-th interval by

for each $j\in J$. Since the conditional survival function $S^{C}(t|\mathbf{x})=P(T>t\,|\,\mathbf{X}=\mathbf{x})$ is equivalent to $\exp\big{(}-\int_{0}^{t}\lambda^{C}(s|\mathbf{x})\,\mathrm{d}s\big{)}$ in the Cox model, we note that $-\log\big{(}1-\lambda_{j}(\mathbf{x})\big{)}=\int_{t_{j-1}}^{t_{j}}\lambda^{C}(s|\mathbf{x})\,\mathrm{d}s$ represents the cumulative conditional hazard on the interval $[t_{j-1},t_{j})$ and, under the proportional hazards assumption, it can be shown that

where $\lambda_{0j}=1-\exp\big{(}-\int_{t_{j-1}}^{t_{j}}\lambda_{0}^{C}(s)\,\mathrm{d}s\big{)}$ is the associated cumulative baseline hazard on $[t_{j-1},t_{j})$. Thus, the cumulative hazard model (2.6) is directly connected with the discrete-time survival model (2.1) under the complementary log transformation, and the two models have the same regression coefficient $\boldsymbol{\beta}$.

Here, we evaluate the empirical performance of the proposed mean score estimator for the discrete survival time setting via a simulation study. We also evaluate the degree to which the proposed adaptive validation design improves estimation performance compared to simple random sampling and a balanced design for several scenarios.

We consider the conditional hazard model (2.1) with the complementary log transformation $g_{2}(u)=-\log(1-u)$, where $\boldsymbol{\beta}=(\log(1.5),\log(0.7),\log(1.3),-\log(1.3))^{\top}$. We assume survival status is observed at discrete times $0<t_{1}<t_{2}<\cdots<t_{10}<\infty$. As previously mentioned in Section 2.3, this model will estimate the same $\boldsymbol{\beta}$ as in the underlying continuous-time Cox proportional hazards model. We first generate a four dimensional covariate vector $\mathbf{X}=(X_{1},\ldots,X_{4})^{\top}$, which consists of both continuous and binary variables.

We simulate correlated covariates with a unit scale between $0$ and $1$, by first considering a multivariate normal random vector $\mathbf{W}=(W_{1},\ldots,W_{4})^{\top}$ with zero mean and $\textrm{Cov}(W_{j},W_{k})=0.3^{|j-k|}$, so that we put $X_{j}=Q_{j}(\Phi(W_{j}))$ for $j=1,2$ and $X_{j}\sim\textrm{Bernoulli}(\Phi(W_{j}))$ for $j=3,4$, where $\Phi$ is the cumulative standard normal distribution function, and $Q_{1}$ and $Q_{2}$ are quantile functions of the beta distribution with pairs of the shape and rate parameters $(2,1.5)$ and $(3,3)$, respectively. By doing this, all $X_{j}$’s are correlated with each other, and particularly $X_{1}$ and $X_{2}$ are marginally beta random variables with $\textrm{Corr}(X_{1},X_{2})\approx 0.290$. Since continuous covariates are generally bounded in practice, the simulated continuous covariates $(X_{1},X_{2})^{\top}$ represent standardized covariates over the range of observations into unit intervals.

For discrete survival outcomes, we note that $P(T=t_{j}\,|\,\mathbf{X}=\mathbf{x})=\lambda_{j}(\mathbf{x})\prod_{k=1}^{j-1}(1-\lambda_{k}(\mathbf{x}))$ enables us to generate the discrete survival outcomes as multinomial random variables associated with covariates. To simplify the censoring mechanism in our simulation, we set a fixed censoring time $C=t_{6}$ as the maximum follow-up for all subjects, so that we observe a truncated survival time $Y=\min\{T,t_{6}\}$ and $\Delta=\mathbb{I}(T\leq t_{6})$, It is worth mentioning that the proposed method also allows a random censoring time $C$. An application of the method with a more complex censoring mechanism will be considered in Section 4.2.

Next, we generate $N$-random copies $\mathcal{X}_{N}=\{(Y_{i},\Delta_{i},\mathbf{X}_{i}):1\leq i\leq N\}$ of $(Y,\Delta,\mathbf{X})$, which will inform the benchmark full cohort analysis. We consider full cohort sizes $N=2000,4000,8000$ in the simulation. By choosing different baseline hazards, we considered three scenarios of overall censoring rates with $30\%$, $50\%$ and $70\%$. Here, we mainly refer to simulation results when the censoring rate is $50\%$, since we found that this setting fundamentally provides similar lessons with the other scenarios which can be found in the Supplementary Material Tables B.1 and B.2. Figure 1 illustrates the conditional hazard, survival and probability functions in our simulation settings.

We assume that direct observations of $X_{1}$ are not available in the phase one sample of the cohort, but instead one observes a discretized and error-prone exposure $Z$ such that

where $X_{1}^{\ast}=X_{1}+\varepsilon$ is perturbation of $X_{1}$ with an independent measurement error $\varepsilon\sim N(0,0.1^{2})$. Let $Z^{o}$ be the true discretization of $X_{1}$ defined similarly to $Z$ by replacing $X_{1}^{\ast}$ with $X_{1}$. For the assumed parameter values, the discordance or misclassified rate between $Z^{o}$ and $Z$ was $P(Z^{o}\neq Z)\approx 0.284$ and this shows that $Z$ is not only a discrete but also an error-prone surrogate of $X_{1}$, and potentially associated with $(X_{2},X_{3},X_{4})$. Finally, we set $\mathcal{X}_{\textrm{I},N}=\{(Y_{i},\Delta_{i},Z_{i}):1\leq i\leq N\}$ to be the phase one sample available on all subjects.

Efficient estimation of the regression coefficient of $X_{1}$ is of interest and the optimal sampling allocations are designed to minimize the variance of the $\beta_{1}$ estimate in our simulation study. Since the optimal sampling design (3.1) depends on nuisance parameters defined by the population structure, namely $I_{V}$ and $\mathrm{Var}\big{(}U_{1}(\boldsymbol{\theta})|y,\delta,z\big{)}$ in Theorem 1, we approximate their true values with empirical estimates obtained from an externally generated large sample of size $N_{0}=10^{4}$ that was independent from the full cohort $\mathcal{X}_{N}$. We define the optimal sampling strategy based on these values as the oracle procedure and refer to Supplementary Material Section A.2 for some technical details for the derived optimal allocation. In practice, however, the oracle procedure is infeasible and so we are interested in evaluating the adaptive sampling design as described in Section 3. For this, we first sample a pilot validation subset $\mathcal{X}_{\textrm{II},n}^{\textrm{Pilot}}$ and estimate the optimal sampling design together with the phase one sample $\mathcal{X}_{\textrm{I},N}$. To accommodate all possible strata information in the first phase sample, we employ balanced sampling for the pilot study, with $n^{\textrm{Pilot}}(Y_{i},\Delta_{i},Z_{i})$ equal for all $\{Y_{i},\Delta_{i},Z_{i}\}$. We then estimate $I_{V}$ and ${\mathrm{Var}}\big{(}U_{1}(\boldsymbol{\theta})|y,\delta,z\big{)}$, as outlined at the end of the previous section. Next, we draw an additional validation subset $\mathcal{X}_{\textrm{II},n}^{\textrm{Adapt}}$ for each stratum following (3.2), and the mean score method is finally applied to the two-stage analysis for the phase one sample $\mathcal{X}_{\textrm{I},N}$ and the adaptive phase two sample $\mathcal{X}_{\textrm{II},n}=\mathcal{X}_{\textrm{II},n}^{\textrm{Pilot}}\cup\mathcal{X}_{\textrm{II},n}^{\textrm{Adapt}}$. We considered validation subset sizes of $n=200,400,800$ and took equal proportions for the pilot and adaptive samples.

Due to saturated strata, there might be some remaining validation size to be allocated, for example $n/2-\sum_{(y,\delta,z)}n^{\textrm{Pilot}}(y,\delta,z)>0$ in the pilot study. In this case, we randomly select unvalidated subjects for the remaining allocation, which is also similarly applied to the adaptive sample. We note that the proportion of the pilot sample size to the adaptive entails a trade-off between precision of the nuisance parameters needed for optimal sampling design and efficiency gains from the adaptive validation when the final phase two sample size is fixed. Some preliminary simulations showed that the adaptive sampling design with nearly equal sizes of the pilot and adaptive samples usually produced robust and efficient estimates (data not shown), which is similar to observations made by McIsaac and Cook (2015) for the two-stage analysis with binary outcomes.

We compare our proposed adaptive sampling method, which we refer to here as mean score adaptive (MS-A), to the mean score method using the oracle design (MS-O) and to some other standard estimation methods for two-phase designs. The complete case analysis of fully randomly selected $n$-validation sample (CC-SRS) will give unbiased results in likelihood-based inference. Since CC-SRS does not use auxiliary information of the first stage sample, we examine if the optimal design for the mean score method improves estimation performance of CC-SRS and evaluate efficiency gains from the two-stage analysis (MS-SRS). We also conduct design-based estimation with balanced sampling such that validation size is equally allocated to each $(y,\delta,z)$-stratum of the first stage sample. Similarly to the proposed adaptive procedure, if there are some remaining individuals to be allocated after balanced sampling, due to saturated strata, we randomly sample the remaining from the unvalidated subjects to yield a final total phase two sample of $n$ individuals. The design-based estimation with balanced sampling is given by a Horvitz-Thompson type estimator (MS-BAL), where sampling proportions of validation within strata are used for inverse probability weights as in (2.4). For our setting, we note that the inverse probability weighted (IPW) estimator is technically the same with the design-based mean score estimator which incorporates the balanced sampling weights pre-specified in the two-phase analysis (McIsaac and Cook, 2014). We implement the proposed mean score estimator with adaptive sampling described above, estimating the necessary nuisance parameters from the validation subset (MS-A) and the when using the oracle design (MS-O). Finally, we consider the full cohort analysis (Full-CC) based on fully observed data, as the benchmark performance, which empirically gives the upper bound of efficiency for the two-phase analysis, since the Full-CC uses the complete covariate information on all subjects.

We investigated two different aspects of varying sample sizes and conducted six simulation scenarios; (i) increasing the phase two sample size from $n=200,400,800$ with a fixed full cohort size $N=4000$, (ii) increasing the full cohort size from $N=2000,4000,8000$ when the phase two sample size is fixed by $n=200$. Table 1 summarizes the relative performance for estimation of $\beta_{1}$. Compared to CC-SRS, MS-SRS had a reduction of the variance, which demonstrates efficiency gains of the mean score method from employing auxiliary information of the phase one sample. For $n=400$ and 800, MS-A was more efficient than MS-SRS, whereas for $n=200$, MS-SRS was more efficient. This suggests that $n=200$ is too small in this setting to gain efficiency using adaptive sampling; this is reasonable as the optimal sampling allocation is based on estimates derived from a pilot sample half that size. MS-BAL had relatively inferior performance compared to MS-SRS and MS-A, and this was true across all scenarios studied. The MS-O had the smallest mean squared error (MSE) in all scenarios, as expected, and it turned out that the superior performance of MS-O mostly came from variance reduction. As the cohort sample size increased, the associated adaptive sampling design (MS-A) behaved closely to MS-O, which indicates that the proposed adaptive approach consistently approximates the oracle procedure. In the scenarios studied, increasing the phase two sample size was much more beneficial for efficiency gains for MS-A and MS-O, compared to the sample size increment of the phase one study. For example, MS-A achieved $60.7\%(\approx 1-0.147/0.374)$ variance reduction while the validation rate increased four-fold from $n=200$ to $n=800$ given $N=4000$ in the left three columns of Table 1. On the contrary, increasing the phase one sample size from $N=2000$ to $N=8000$ only produced $11.2\%(\approx 1-0.190/0.214)$ improvement in the right three columns when $n=400$ is fixed. This suggests that the performance of the proposed MS-A procedure is sensitive to the size of the phase two study, more so than the total cohort size for a fixed $n$.

Simulation results for all regression parameters and various combinations of $n=200,400,800$ and $N=2000,4000,8000$ can be found in Supplementary Material Tables B.3–B.7. Comparing the relative performance of the methods for the other regression parameters, including those for the discrete baseline hazards, we found that the MS-O outperformed the other estimators for all model parameters. The oracle sampling design is infeasible in most practical situations; however, the two mean score estimators consistently showed the best performance for the target parameter among all other practical competitors. MS-A generally achieved the minimum mean squared errors for all parameters when the phase two sample size exceeded 200, compared to the other practical estimators, for the scenarios studied (Supplementary Material Tables B.3–B.7). For the setting where $n=200$, the MS-SRS performed the best amongst the practical estimators. These results suggest for that MS-A will perform well with respect to MSE for all model parameters, particularly the target and discrete proportional hazard parameters, for robust phase two sample size. For small phase two samples and similar censoring rates, the MS-SRS may be preferred.

Wilms tumor is a rare renal cancer occurring in children, where tumor histology and the disease stage at diagnosis are two important risk factors for relapse and death. We consider data, reported by Kulich and Lin (2004), on 3915 subjects from two randomized clinical trials from the National Wilms Tumor Study (NWTS) (D’angio et al., 1989; Green et al., 1998). There are two measures of tumor histology, which is classified as either favorable (FH) or unfavorable (UH), one by a local pathologist and the other by an expert pathologist from a central facility. Because of the rarity of disease, local pathologists may be less familiar with Wilms Tumor and their assessment is subject to misclassification. The central assessment is considered to be the gold standard, or true histology in our analysis, and the local evaluations are considered surrogate observations for the central evaluation. Since histology for all subjects was validated by the central laboratory, NWTS data has been widely used to evaluate two-phase sampling methodology, such as by Breslow and Chatterjee (1999), Kulich and Lin (2004) and Lumley (2011).

We demonstrate our proposed mean score method for discrete-time survival in an analysis of relapse in this NWTS cohort. We assume that the local histology is available for all subjects in a first phase sample and that only a sub-cohort was sampled for evaluation by the central pathologist in a second phase sample. Specifically, we are interested in the proportional hazards discrete-time survival model (2.1) for time to relapse, under the complementary log transformation $g_{2}(u)=-\log(1-u)$, in order to evaluate the risk associated with unfavorable central histology, late (III/IV) disease stage versus the early (I/II) stage, age at diagnosis (year), and tumor diameter (cm). For this model, we also include an interaction between histology and stage of disease.

In this cohort, $90\%$ of the $669$ events occurred within the first three years of diagnosis, while less than $5\%$ of non-relapsed subjects were censored in the same period. Based on this observation, we first define a modified, or reduced, cohort to include only patients who had an event or were fully followed up in the first 3 years, so that censoring only occurs at the third year ($82.2\%$), the assumed end of the study. This modified cohort included $N=3757$ subjects and is used in our NWTS data analysis. We consider the regression coefficients from the discrete-time survival analysis of the modified full cohort as the reference values. We took this approach first out of concern that the large number of nuisance parameters introduced by the small number of individuals with intermittent censored outcomes, due to the added strata for a binary outcome at each failure time interval, might adversely affect the mean score method. We will evaluate efficiency gains of the optimal mean score design, based on the discretized survival outcomes, for the continuous-time analysis in Section 4.3. Finally, we also conducted an analysis of all individuals, regardless of the time of censoring, and discuss design issues regarding how to handle intermittent censoring in our framework in Section 4.4.

We first discretize the continuous event time into six 6-month intervals, so that we model the hazard of relapse during the first three years after diagnosis. As in Section 4.1, we consider four different sampling scenarios for the phase two subsample: simple random sampling, balanced allocation across strata, the adaptive and optimal (oracle) mean score designs, where the last three employ stratified sampling based on the phase one sample. To evaluate the efficiency gains of the proposed mean score approach, we performed the two-stage analysis, with a phase two sample of $n=400$, 1000 times. For implementation of the optimal design, we estimate parameters in (3.1) using the oracle procedure by using the central histology records of the full cohort. We refer to Section 4.1 for further details of implementation. Bias and efficiency are calculated using the estimates of the full cohort analysis with the central histology as the reference.

Table 2 demonstrates the performance of the different methods, where efficient estimation of the interaction between unfavorable histology and late stage of the disease is of main interest. MS-A outperformed the other practical competitors: CC-SRS, MS-SRS and MS-BAL. For example, MS-A had a $38.7\%(\approx 1-0.425/0.693)$ variance reduction compared to CC-SRS for estimating the interaction term. We also note that performance of MS-A was pretty close to the oracle procedure, MS-O; MS-A was only $7\%$ less efficient than MS-O for estimating the regression coefficient of interest. Furthermore, across all parameters, MS-A achieved the smallest mean squared error among the practical competitors, not only for baseline hazards estimation but also regression coefficients. Overall, the proposed method performed better than the other two-phase estimators, with both an efficient design and by incorporating auxiliary information from the phase one sample.

In Figure 2, we examine the relative performance of the different estimators for each of the regression parameters and phase two sample sizes $n=200,400,800$. In the top-left panel (a), we show the estimation results for the regression coefficient for unfavorable histology (UH) over 1000 repetitions of subsampling, when the adaptive and optimal mean score allocations were designed for the efficient estimation of UH. The rest of the panels are similar, showing results for the other regression coefficients they were set as targets, namely when the MS-A and MS-O designs were for efficient estimation of (b) late stage of the disease, (c) age of diagnosis (year), (d) tumor diameter (cm) and (e) the interaction between UH and stage of the disease, respectively. MS-O consistently showed superior performance for all phase two ample sizes; however, MS-A again had the smallest mean squared error among the practical methods. Further, the performance of MS-A tended to be close to MS-O as the validation size increased.

In many practical settings, the continuous-time Cox model will be the analysis of primary interest. We further investigate benefits of the proposed mean score sampling design, when the phase two estimating equation employs the continuous-time Cox proportional hazards model. In Section 2.3, we discussed the direct connection between the the continuous-time Cox model (2.1) and the analogous grouped discrete-time model with the complementary log transformation, in the sense that the two models have the same regression coefficient $\boldsymbol{\beta}$. We can conduct the two-phase analysis of the Cox model with the proposed optimal mean score sampling method derived for the parameter of interest in the discretized model.

Table 3 demonstrates performance comparison of different sampling methods when they were applied to the two-phase analysis of the usual Cox model based on 1000 repeated phase two samples in the NWTS data. Unlike the previous analysis in Table 2 or Figure 2, we assumed that the continuous relapse times were observed for the full cohort in the phase one study and employed the previously developed optimal and adaptive mean score design for the discrete-time survival model (2.1) only to design the allocation of the phase two sample. That is, we first calculated the sampling probabilities and associated inverse probability weights for the phase two sample, and then the design-based estimates of the continuous-time Cox model were reported in Table 3. We assumed efficient estimation of the interaction effect between unfavorable histology and disease stage was of primary interest in the adaptive and oracle sampling designs. We used the survival package in R (Therneau and Lumley, 2019) and applied the inverse probability weights of the phase two sample with the weights option of the coxph function. As shown in Table 3, the adaptive sampling design for the discrete-time analysis also provided efficiency gains for the continuous-time analysis. For example, the variance reduction was about $26\%$ compared to both the simple random sampling and the balanced design. In this example, the adaptive design happened to slightly outperform the oracle one, but we note that these two allocations were for the optimal design for the discrete-time and not the continuous survival model (2.1). We found that the proposed adaptive sampling design provided $20\%$–$25\%$ of efficiency gains in each case varying the target parameter to be (a) unfavorable histology (UH), (b) late stage of the disease, (c) age of diagnosis (year) or (d) tumor diameter (cm) (data not shown).

The proposed method in Section 3 can be applied without modification for the setting where there is random intermittent censoring, as well as censoring at the end of the study; the only difference being that there will be more strata at the intermittent event times. The mean score estimation depends on nonparametric estimation of the probability distribution for the phase two variables conditional on each discrete value of the phase one surrogate. We note that $\hat{\pi}(Y_{i},\Delta_{i},\mathbf{Z}_{i})$ in Section 2.2 is the empirical estimate of the sampling probability of the $i$-th individual selected into the validation subset, which may suffer from the curse of dimensionality as the number of phase one strata $(Y_{i},\Delta_{i},\mathbf{Z}_{i})$ increases. One could consider an increasing number of discretized intervals to approximate the continuous-time points, but the number of unique continuous-time survival outcomes observed will increase as the sample size increases. For this reason, in our data example we first studied the two-phase analyses of the discrete-time survival models in the previous sections where individuals were right-censored only at the sixth time point, an induced end of follow-up, which is equivalent to a fixed Type I censoring when individuals who were intermittently censored were excluded from analysis. Even in this simple setting in Section 4.1, there were $(7\times 4)$-strata for the phase two sampling, resulting from a discrete surrogate $Z$ having four categories. If we further considered random right-censoring for this setting, we would have to estimate sampling probabilities for $(6\times 2\times 4)$-strata. The larger number of associated nuisance parameters in this case for MS-A estimator may require larger phase two samples to achieve the expected efficiency gains due to unstable nuisance parameter estimation in small strata. Indeed, in our simulation study with Type I censoring only, we saw that when $n=200$, the MS-A did not outperform MS-SRS.

We suggest a simple strategy for the proposed adaptive mean score design aimed at under-sampling less informative strata in the pilot study, which may preserve efficiency gains by providing more precise nuisance parameters for the more informative strata. For example, under the random right-censoring assumption, individuals censored before the end of the follow-up period should generally be be less informative than a non-censored individuals with similar covariates for that same period. For a fixed phase two sample size, such under-sampling may enable us to re-allocate the pilot sample for the MS-A so that relatively more informative groups can be up-weighted. In fact, we anticipate that the optimal adaptive design would allocate very few individuals in these strata, so we seek to avoid putting too many individuals into these strata in the pilot. Supplementary Material Tables B.8 and B.9 provide analogous results for the discrete-time analysis in Table 2 and continuous time analysis of Table 3 for the NWTS, but we analyzed the full cohort ($N=3915$) and under-sampled the censored individuals not relapsing before the end of the follow-up period in the balanced and the pilot samples. Specifically, we allocated a small number of the pilot sample size for individuals censored before the end of the follow-up period (i.e., $y<3$ and $\delta=0$) at each level of histology, and the remaining allocation was equally distributed to the other strata in the pilot sample. In this particular example, we set the under-sampling size to 4, which was approximately half of the balanced sampling size for each stratum in the pilot study with $n=400$. Indeed, we found that the naive balanced sampling of all strata in the pilot often over-sampled censored-groups compared to the oracle design (data not shown), and consequently this simple remedy made the modified adaptive sampling design be closer to the oracle design than did the balanced pilot sample with the proposed adaptive design. As seen in Supplemental Table B.8, the modified approach produced a $14\%$ variance reduction compared to both the fully balanced and the MS-A with a balanced pilot sample. Additionally, the precision for all estimates, relative to the analogous estimates in Tables 2 and 3, benefited from including the intermittently censored $158$ ($=3915-3757$) individuals who had been excluded in Section 4.2, with more gains seen for the MS-BAL and MS-A estimators. We hypothesize that this modification will allow for more robust efficiency gains for the proposed adaptive design in other settings with intermittent censoring. Supportive simulation studies that examine the performance of MS-A given the anticipated phase two sample size and other study parameters may provide useful insights to guide refinements of this strategy for a given setting.

We assume there exists $J\geq 1$ such that $P(Y\leq t_{J})$ with probability one. Let $\boldsymbol{\theta}_{0}\in\Theta_{0}$ be the true parameter for some open $\Theta_{0}\subset\mathbf{R}^{J+d}$ and further assume the following conditions:

• A1.

  $\mathrm{E}|L_{1}(\boldsymbol{\theta};Y,\Delta,\mathbf{X})|<\infty$ at $\boldsymbol{\theta}_{0}$.

• A2.

  $L_{1}(\boldsymbol{\theta};Y,\Delta,\mathbf{X})$ is three-times continuously differentiable and bounded away from 0 near $\boldsymbol{\theta}_{0}$ (a.s.).

• A3.

  $I_{V}=-\mathrm{E}\big{[}\frac{\partial^{2}}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\top}}\log L_{1}(\boldsymbol{\theta};Y,\Delta,\mathbf{X})\big{]}$ is positive definite at $\boldsymbol{\theta}_{0}$.

• A4.

  $\pi(y,\delta,\mathbf{z})>0$ for all $(y,\delta,\mathbf{z})$.