What is already known

• • Heterogeneous cutoff values over studies often occur in meta-analysis examining the association of a continuous biomarker with a binary outcome (diagnostic study) and with a time-to-event outcome (prognostic study).

• • For diagnostic meta-analysis, the summary receiver operating characteristics (SROC) curve and its area under the curve are widely recognized as a useful tool and used to handle the issue of heterogeneous cutoff values.

• • For prognostic meta-analysis with time-to-event outcomes, the time-dependent SROC method, which is an extension of the SROC method in diagnostic meta-analysis, has been developed by Hattori and Zhou (2016, Statistics in Medicine 35(26), 4746–4763) and is useful to address the issue of heterogeneous cutoff values.

• • In meta-analysis, large studies or studies with significant results are more likely to be published and observed for meta-analysis; synthesizing only these studies may result in reporting bias, also known as, publication bias or small-study effects.

• • Reporting bias is unavoidable and may also induce biased estimates of the time-dependent SROC method by the bivariate normal model of Hattori and Zhou (2016, Statistics in Medicine 35(26), 4746–4763).

What is new

• • For prognostic studies with time-to-event outcomes, the selective publication process of studies and the selective reporting of the Kaplan–Meier estimates can be influenced by the significance of the hazard ratios, equivalently, the significance of the log-rank tests.

• • To model selective publication processes driven by the log-rank test, inference of trivariate model has been successfully established, which is essentially beyond the sensitivity analysis method for reporting bias in meta-analysis of diagnostic studies.

• • The proposed model was constrained by the marginal selection probability as a sensitivity parameter; the impact of reporting bias on the time-dependent SROC method can be evaluated by varying the marginal selection probability within (0,1).

Potential impact for Research Synthesis Methods readers outside the authors’ field

• • We proposed the first method dealing with reporting bias on the time-dependent SROC method for meta-analysis of prognosis studies with time-to-event outcomes; then, users could evaluate the robustness of the estimates when evaluating the prognostic capacity of the biomarker.

• • By incorporating the ability to address reporting bias, our method enhanced the utility and applicability of the time-dependent SROC method of Hattori and Zhou (2016, Statistics in Medicine 35(26), 4746–4763) for meta-analysis of prognosis studies.

3.1 Notations and data structure

Suppose that S prognosis studies are conducted to evaluate the association between the expressions of biomarker and subjects’ time-to-event outcomes. Among the S studies, data of N prognosis studies are published and used in meta-analysis to summarize the prognostic capacity, while data of $S-N$ studies are missing. In this section, we do not consider the existence of reporting bias, that is, data of N published studies consist of the population or random sample from the S studies. Since we focus on the HZ model to synthesize data of the published studies, we follow the notations in Hattori and ZhouReference Hattori and Zhou ¹⁵ to introduce the latent individual patient data (IPD) in each study and the observable data for meta-analysis as follows.

Suppose that each prognosis study $i\left (i = 1, 2, \dots , N\right )$ includes $n^{(i)}$ subjects, and the subjects are assumed to be random samples from the population of interest. For each subject, let $\tilde X$ be the baseline measurement of biomarker, $\tilde T$ the failure time, $\tilde C$ the right-censored time, and $\tilde Y = \mathrm {min}\left (\tilde T, \tilde C\right )$ the follow-up time. We use tilde to denote these variables for individuals, which are unextractable from the literature for meta-analysis. We assume that the distribution of $\tilde C$ is identical across N studies as well as $\tilde C {\perp \!\!\!\perp } \tilde X$ , and $\tilde T {\perp \!\!\!\perp } \tilde C\mid \tilde X$ , where ${\perp \!\!\!\perp }$ indicates variables are independent. These assumptions lead to $\tilde T {\perp \!\!\!\perp } \tilde C$ .

For each study i, let $v^{(i)}$ denote the study-specific cutoff value, which is not necessarily reported. This cutoff value separates subjects into the low expression group, denoted by $Z=0$ , if $\tilde X \le v^{(i)}$ or the high expression group ( $Z=1$ ) if $\tilde X> v^{(i)}$ . The survival functions of the low and high expression groups are respectively denoted and defined by

$$ \begin{align*} S_0^{(i)}\left(t\right)&=P\left(\tilde T> t \mid \tilde X \le v^{(i)}\right)\\ S_1^{(i)}\left(t\right)&=P\left(\tilde T > t\mid \tilde X > v^{(i)}\right). \end{align*} $$

The disease status at time t is defined by the counting process $D(t)=1$ if $\tilde T\le t$ or $D(t)=0$ if $\tilde T> t$ , indicating that subjects have an event before time t or survive after t, respectively. To facilitate the understanding, we showed the number of subjects given different disease statuses and expression groups in Table S1 in the Supplementary Material. Analogous to the diagnostic study, sensitivity and specificity at time t are respectively denoted and defined by

$$ \begin{align*} \text{se}\left(x,t\right) &= P\left(\tilde X> x\mid \tilde T \le t\right)\\ \text{sp}\left(x,t\right) &= P\left(\tilde X \le x\mid\tilde T > t\right). \end{align*} $$

For meta-analysis, the following summarized data are required and can be extracted from published literatures. For study $i~(i=1,2,\dots , N)$ , let $n_0^{(i)}$ and $n_1^{(i)}$ denote the number of subjects separated into the low and high expression groups, respectively; the total number of subjects is $n^{(i)} = n_0^{(i)} + n_1^{(i)}$ . Let $\hat S_0^{(i)}\left (t\right )$ and $\hat S_1^{(i)}\left (t\right )$ denote the KM estimates of the low and high expression groups, respectively, and they can be extracted from the reported plots of the KM curves at the partition of time interval $\left [0,t\right ]$ : $0=t_0<t_1<\dots <t_K=t$ . Let $\hat \mu _{\text {lnHR}}^{(i)}$ denote the reported lnHR between high versus low expression groups and $\hat s_{\text {lnHR}}^{(i)}$ the corresponding SE; both are estimated by the Cox model.Reference Cox ²⁶ Among N studies, the sample medians of the follow-up time over the total subjects can also be extracted from some studies. These medians of follow-up time are used to estimate the censoring distribution in estimating the SEs of $\hat {S}_0^{(i)}(t)$ and $\hat {S}_1^{(i)}(t)$ . (See Section 3.2 in Hattori and ZhouReference Hattori and Zhou ¹⁵ for more details.) We take the example of meta-analysis of Ki67 to illustrate the data structure with details explained in Section S1 of the Supplementary Material.

3.2 Bivariate normal model

In this section, we review the structure of HZ model and the definitions of $\text {SROC}(t)$ and $\text {SAUC}(t)$ . The HZ model employs the bivariate normal-normal random-effects model to model the empirical pairs of logit-transformed time-dependent sensitivities and specificities, which are estimated using KM estimates from published literatures, allowing for heterogeneity from different sources (e.g., cutoff values, study designs, population, etc.) in time-dependent sensitivities and specificities across studies. It extended the model of Reitsma et al.Reference Reitsma, Glas, Rutjes, Scholten, Bossuyt and Zwinderman ⁵ which has been widely used in meta-analysis of diagnostic studies.

According to Hattori and Zhou,Reference Hattori and Zhou ¹⁵ with the Bayes rule and simple algebraic manipulations, the true sensitivity and specificity at time t of study i can be respectively re-expressed by

(1) $$ \begin{align} \begin{aligned} \text{se}\left(v^{(i)}, t \right) &= \dfrac {\left\{1- S_1^{(i)}\left(t\right)\right\}q_1^{(i)}} {\left\{1- S_0^{(i)}\left(t\right)\right\}q_0^{(i)} + \left\{1- S_1^{(i)}\left(t\right)\right\}q_1^{(i)}}\\ \text{sp}\left(v^{(i)}, t\right) &= \dfrac { S_0^{(i)}\left(t\right)\cdot q_0^{(i)}} { S_0^{(i)}\left(t\right)\cdot q_0^{(i)} + S_1^{(i)}\left(t\right)\cdot q_1^{(i)}}, \end{aligned} \end{align} $$

where $q_0^{(i)} = P\left (\tilde X\le v^{(i)}\right )$ and $q_1^{(i)} = P\left (\tilde X>v^{(i)}\right )$ . With the data extracted from collected studies, $q_0^{(i)}$ and $q_1^{(i)}$ are estimated by $\hat q_0^{(i)} = n_0^{(i)}/n^{(i)}$ and $\hat q_1^{(i)} = n_1^{(i)}/n^{(i)}$ , respectively; the time-dependent sensitivity and specificity in equation (1) are consistently estimated by substituting $\left (\hat q_0^{(i)}, \hat q_1^{(i)}, \hat S_0^{(i)}, \hat S_1^{(i)}\right )$ for $\left (q_0^{(i)}, q_1^{(i)}, S_0^{(i)}, S_1^{(i)}\right )$ . The resulting consistent estimators are denoted by $\hat {\textrm {se}}\left (v^{(i)}, t \right )$ and $\hat {\textrm {sp}}\left (v^{(i)}, t \right )$ , respectively.

The HZ model, following Reitsma et al.,Reference Reitsma, Glas, Rutjes, Scholten, Bossuyt and Zwinderman ⁵ employs a logit-transformation to map the time-dependent sensitivity and specificity to the real number line $(-\infty , \infty )$ . While other transformations are also applicable,Reference Doebler, Holling and Böhning ²⁷ logit-transformation is widely accepted. Let $\boldsymbol \mu ^{(i)} = \left (\mu _{\text {se}}^{(i)},\mu _{\text {sp}}^{(i)}\right )^\top $ denote the true time-dependent sensitivity and specificity pair of study i on the logit scale, that is,

$$ \begin{align*} \begin{aligned} \mu_{\text{se}}^{(i)}=\text{logit}\left\{\text{se}\left(v^{(i)},t\right)\right\}\text{ and } \mu_{\text{sp}}^{(i)}=\text{logit}\left\{\text{sp}\left(v^{(i)},t\right)\right\}, \end{aligned} \end{align*} $$

where $\text {logit}=\log x-\log (1-x)$ . At the between-study level, it is assumed that

(2) $$ \begin{align}{\boldsymbol \mu}^{(i)} \sim N_2 \left (\boldsymbol \mu, \boldsymbol \Omega \right ) \text{ with } \boldsymbol \Omega= \begin{bmatrix} \tau_{\text{se}}^2 & \tau_{\mathrm{se,sp}}\\ \tau_{\mathrm{se,sp}} & \tau_{\text{sp}}^2 \end{bmatrix}, \end{align} $$

where $N_2$ denotes the bivariate normal distribution, and $\boldsymbol \mu = \left (\mu _{\text {se}}, \mu _{\text {sp}}\right )^\top $ is the overall mean at time t across multiple prognosis studies; $\boldsymbol \Omega $ is the between-study variance-covariance matrix of ${\boldsymbol \mu }^{(i)}$ , where $\tau _{\text {se}}^2$ and $\tau _{\text {sp}}^2$ are the variance of $\mu _{\text {se}}^{(i)}$ and $\mu _{\text {sp}}^{(i)}$ , respectively, and $\tau _{\mathrm {se,sp}}$ the covariance between them.

Let $\hat {\boldsymbol \mu }^{(i)} = \left (\text {logit}\left \{\hat {\textrm {se}}\left (v^{(i)},t\right )\right \}, \text {logit}\left \{\hat {\textrm {sp}}\left (v^{(i)},t\right )\right \}\right )^\top $ denote the consistent estimates of $\boldsymbol \mu ^{(i)}$ in each study. Hattori and ZhouReference Hattori and Zhou ¹⁵ showed that $\hat {\boldsymbol \mu }^{(i)}$ has the following asymptotic distribution at the within-study level:

(3) $$ \begin{align} \hat{\boldsymbol \mu}^{(i)}\mid \boldsymbol \mu^{(i)} \sim N_2 \left (\boldsymbol \mu^{(i)}, \left.{ {\mathbf H}^{(i)}}\right/{n^{(i)}} \right ) \text{ with } {\mathbf H}^{(i)} = \begin{bmatrix} \left\{\sigma_{\text{se}}^{(i)}\right\}^2 & \sigma_{\mathrm{se,sp}}^{(i)}\\ \sigma_{\mathrm{se,sp}}^{(i)} & \left\{\sigma_{\text{sp}}^{(i)}\right\}^2 \end{bmatrix}, \end{align} $$

where ${\mathbf H}^{(i)}$ indicates the within-study asymptotic variance-covariance matrix. The detailed expression of ${\mathbf H}^{(i)}$ is presented in equation (S1) of the Supplementary Material. The parameter ${\mathbf H}^{(i)}$ in model (3) can be replaced with its consistent estimator $\hat {\mathbf H}^{(i)}$ estimated by the Greenwood formula with median follow-up time, which follows the convention that the within-study variance–covariance matrix is known in meta-analysis. (See Section 3.2 in Hattori and ZhouReference Hattori and Zhou ¹⁵ for more details.) Combining models (2) and (3) induces the HZ model:

(4) $$ \begin{align} \hat{\boldsymbol \mu}^{(i)} \sim N_2 \left (\boldsymbol \mu, \left. \boldsymbol \Omega + {\hat{\mathbf H}^{(i)}}\right/{n^{(i)}} \right ), \end{align} $$

where $(\boldsymbol \mu , \boldsymbol \Omega )$ are the unknown parameters at time t, and they can be estimated by the maximum likelihood (ML) method; we denote their ML estimators as $(\hat {\boldsymbol \mu }, \hat {\boldsymbol \Omega })$ .

Based on the HZ model, $\text {SROC}(t)$ is derived by taking the conditional expectation of $\mu _{\text {se}}^{(i)}$ given $\mu _{\text {sp}}^{(i)}$ in model (2). Let x denote $1-\text {sp}\left (x,t\right )$ , $\text {SROC}(t)$ is defined by the following time-dependent function:

(5) $$ \begin{align} \text{SROC}(t) = \text{SROC}\left(x, t; \boldsymbol \mu, \boldsymbol \Omega \right) = \text{logit}^{-1}\left[ \mu_{\text{se}} - \dfrac{\tau_{\mathrm{se,sp}}}{\tau_{\text{sp}}^2}\left\{ \text{logit}(x)+\mu_{\text{sp}}\right\}\right]. \end{align} $$

Accordingly, $\text {SAUC}(t)$ is defined by

(6) $$ \begin{align} \text{SAUC}(t) = \text{SAUC}\left(t; \boldsymbol \mu, \boldsymbol \Omega \right) = \int_0^1 \text{SROC}\left(x, t; \boldsymbol \mu, \boldsymbol \Omega\right) dx. \end{align} $$

$\text {SROC}(t)$ and $\text {SAUC}(t)$ can be estimated by replacing the unknown parameters $(\boldsymbol \mu , \boldsymbol \Omega )$ in (5) and (6).

4.1 Trivariate model incorporating the lnHR

In intervention studies, whether or not a study is published is influenced by the significance of the result, such as, the P-value of the t-statistic. Similarly, a prognosis study with time-to-event outcomes is more likely to be published when its log-rank test (or equivalently, the t-statistic of the lnHR without covariates) is significant, and vice versa. The phenomenon of selective publication causes N published studies to be biased sample from the S studies and may induce reporting bias in the estimates of $\text {SROC}(t)$ and $\text {SAUC}(t)$ . To quantify reporting bias on the estimates, we aim to extend the methods of CopasReference Copas ¹⁹ and Zhou et al.Reference Zhou, Huang and Hattori ²⁴ and introduce the selection function of the log-rank test (equivalently, the t-statistic of the lnHR) to model the publication mechanism of prognosis studies; then, the inference of reporting bias is made based on the conditional likelihood of the HZ model given the published studies. However, the HZ model (4) does not involve the lnHR. Thus, before introducing the proposed method, we need to expand the HZ model (4) to correlate the time-dependent sensitivity and specificity with the lnHR for constructing the conditional likelihood in the later section.

To distinguish from some notations for the HZ model (4) in Section 3.2, we let $\hat {\mathbf y}^{(i)}=\left (\hat \mu _{\text {se}}^{(i)},\hat \mu _{\text {sp}}^{(i)}, \hat \mu _{\text {lnHR}}^{(i)}\right )^\top $ denote the consistent estimators of $\boldsymbol \theta ^{(i)} = \left ( \mu _{\text {se}}^{(i)},\mu _{\text {sp}}^{(i)}, \mu _{\text {lnHR}}^{(i)}\right )^\top $ at time t, where $\mu _{\text {lnHR}}^{(i)}$ denotes the lnHR and $\hat \mu _{\text {lnHR}}^{(i)}$ the empirical lnHR estimated by the Cox model from the published prognosis studies. At the between-study level, it is assumed that

(7) $$ \begin{align} {\boldsymbol\theta}^{(i)} \sim N_3\left ({\boldsymbol\theta}, \boldsymbol \Psi \right ) ~\text{with}~ \boldsymbol \Psi = \begin{bmatrix} \psi_{\text{se}}^2 & \psi_{\mathrm{se,sp}} & \psi_{\mathrm{se,lnHR}} \\ \psi_{\mathrm{se,sp}} & \psi_{\text{sp}}^2 & \psi_{\mathrm{sp,lnHR}} \\ \psi_{\mathrm{se,lnHR}} & \psi_{\mathrm{sp,lnHR}} & \psi_{\text{lnHR}}^2 \end{bmatrix}, \end{align} $$

where $\boldsymbol \theta = \left (\theta _{\text {se}}, \theta _{\text {sp}}, \theta _{\text {lnHR}}\right )^\top $ are the common means, and $\boldsymbol \Psi $ indicates the between-study variance-covariance matrix of ${\boldsymbol \theta }^{(i)}$ ; the diagonal elements are the corresponding variances of $\theta _{\text {se}}, \theta _{\text {sp}}$ , and $\theta _{\text {lnHR}}$ and the others are the covariances between each two of them.

At the within-study level, we can prove that

(8) $$ \begin{align} \hat{\mathbf y}^{(i)}\mid \boldsymbol\theta^{(i)} \sim N_3 \left ( \left.{\boldsymbol\theta^{(i)}}, {{\boldsymbol \Sigma}^{(i)}}\right/n^{(i)} \right ) \text{ with } {\boldsymbol \Sigma}^{(i)} = \begin{bmatrix} \mathbf H^{(i)} & \boldsymbol \Sigma_{12}^{(i)}\\ {\boldsymbol \Sigma_{12}^{(i)}}^\top & \left\{ \sigma_{\text{lnHR}}^{(i)}\right\}^2 \end{bmatrix}, \end{align} $$

where $N_3$ denotes the trivariate normal distribution, $\boldsymbol \Sigma ^{(i)}$ the within-study asymptotic variance–covariance matrix, ${\mathbf H}^{(i)}$ the variance–covariance matrix in model (3); let $\boldsymbol \Sigma _{12}^{(i)}=\left (\sigma _{\mathrm {se,lnHR}}^{(i)}, \sigma _{\mathrm {sp,lnHR}}^{(i)}\right )^\top $ , where $\sigma _{\mathrm {se,lnHR}}^{(i)}$ is the covariance between $\hat \mu _{\text {se}}^{(i)}$ and $\hat \mu _{\text {lnHR}}^{(i)}$ and $\sigma _{\mathrm {sp,lnHR}}^{(i)}$ the covariance between $\hat \mu _{\text {sp}}^{(i)}$ and $\hat \mu _{\text {lnHR}}^{(i)}$ , and $\left \{ \sigma _{\text {lnHR}}^{(i)}\right \}^2$ the variance of $\hat \mu _{\text {lnHR}}^{(i)}$ . The proof of the asymptotic distribution of $\hat {\mathbf y}^{(i)}$ (8) is presented in Section S2 of the Supplementary Material. In (8), $\boldsymbol \Sigma ^{(i)}$ can be replaced by its consistent estimator $\hat {\boldsymbol \Sigma }^{(i)}$ , according to the convention in meta-analysis that the within-study variance–covariance matrix is known. Combining models (8) and (7), we derives the marginal distribution of $\hat {\mathbf y}^{(i)}$ :

(9) $$ \begin{align} \hat{\mathbf y}^{(i)} \sim N_3 \left (\boldsymbol \theta, \left. \boldsymbol \Psi + {\hat{\boldsymbol \Sigma}^{(i)}}\right/{n^{(i)}} \right ). \end{align} $$

4.2 Selection functions on the significance of the lnHR

As aforementioned, we conjectured that the published N studies are subject to selective publication in that study with significant lnHR (or small P-value of the log-rank test) is more likely to be published. The selective publication process is modeled by the probability of a study being selected given its t-statistic of the lnHR, denoted by $t_{\text {HR}}^{(i)}$ , as shown in the following selection function:

(10) $$ \begin{align} P\left(\text{select}\mid \hat{\mathbf y}^{(i)}, \hat {\boldsymbol \Sigma}^{(i)}\right) = a\left({\hat{\mathbf y}^{(i)}, \hat {\boldsymbol \Sigma}^{(i)}} \right)=a\left( {t_{\text{HR}}^{(i)}}\right), \end{align} $$

where the function $a(\cdot )$ is a non-decreasing function of ${t_{\text {HR}}^{(i)}}$ , that is, the t-statistic of the lnHR. To simplify the inference procedure, we, following Copas,Reference Copas ¹⁹ employ the probit model to $a(\cdot )$ . Thus, equation (10) is defined by:

(11) $$ \begin{align} a\left({t_{\text{HR}}^{(i)}}\right)=\Phi\left(\alpha + \beta\cdot {t_{\text{HR}}^{(i)}}\right), \end{align} $$

where ${t_{\text {HR}}^{(i)}} = \hat \mu _{\text {lnHR}}\left /\hat s^{(i)}_{\text {lnHR}}\right .$ , parameters $\alpha $ and $\beta $ control the probability of selective publication, and $\hat s^{(i)}_{\text {lnHR}}$ denotes the reported SE of the lnHR with $\hat s^{(i)}_{\text {lnHR}} = \left .\hat \sigma _{\text {lnHR}}^{(i)}\right /\sqrt {n^{(i)}}$ . The monotonic property of the probit model links two cases of randomly selective publication: (1) when $\beta =0$ and suppose $\alpha =\Phi ^{-1}\left (p_0\right )$ , the probability of selective publication is independent of the t-statistic, and each study is randomly published from the population with selection probability $p_0$ ; (2) when $\beta \rightarrow \infty $ , each study is published with probability 1.

According to the definition of the probit model, the selection function (11) can be represented by

(12) $$ \begin{align} a\left({t_{\text{HR}}^{(i)}}\right) =\Phi\left(z^{(i)} < \alpha + \beta\cdot {t_{\text{HR}}^{(i)}}\right)= \Phi\left\{ \alpha + \beta\cdot \left(\left.{\hat \mu^{(i)}_{\text{lnHR}}}\right/{\hat s^{(i)}_{\text{lnHR}}}\right) \right\}, \end{align} $$

where $z^{(i)}$ is the standard normal random variable independent of ${t_{\text {HR}}^{(i)}}$ . Based on the trivariate model (9), the marginal distribution of $\hat \mu ^{(i)}_{\text {lnHR}}$ is

$$ \begin{align*} \hat \mu^{(i)}_{\text{lnHR}}\sim N\left({\theta_{\text{lnHR}}}, \psi_{\text{lnHR}}^2 + \left\{\hat s^{(i)}_{\text{lnHR}}\right\}^2 \right). \end{align*} $$

Then, the distribution of ${t_{\text {HR}}^{(i)}}$ is derived by

$$ \begin{align*} N \left( \dfrac{\theta_{\text{lnHR}}}{\hat s^{(i)}_{\text{lnHR}}}, 1 + \left\{\dfrac{\psi_{\text{lnHR}}}{\hat s^{(i)}_{\text{lnHR}}}\right\}^2\right ), \end{align*} $$

and the selection function $a\left ({t_{\text {HR}}^{(i)}}\right )$ in equation (12) can be written into the following selection function:

(13) $$ \begin{align} P\left(\text{select}\mid {\hat {\boldsymbol \Sigma}^{(i)}}\right) = b\left(\hat{\boldsymbol \Sigma}^{(i)}\right) = \Phi\left\{ \dfrac{\alpha + \beta\cdot \left(\theta_{\text{lnHR}}\left/\hat s^{(i)}_{\text{lnHR}}\right.\right) } {\sqrt{1 + \beta^2\cdot \left\{1+ \left({\psi_{\text{lnHR}}}\left/{\hat s^{(i)}_{\text{lnHR}}}\right.\right)^2\right\}}} \right\}. \end{align} $$

4.3 Likelihood based sensitivity analysis

We estimate the parameters in $\text {SROC}(t)$ and $\text {SAUC}(t)$ by maximizing the loglikelihood subject to a certain value of marginal selection probability, denoted by $P(\text {select})$ and defined by

(14) $$ \begin{align} P(\text{select}) = E_P\left\{a\left({t_{\text{HR}}^{(i)}}\right)\right\} = E_P\left\{b\left(\hat{\boldsymbol \Sigma}^{(i)}\right)\right\}. \end{align} $$

We regard $P(\text {select})=p$ as the sensitivity parameter, implying the the expected proportion of the published from the population studies. With various values of p, the changes in the estimates of $\text {SROC}(t)$ and $\text {SAUC}(t)$ imply the impact of reporting bias on them. In this section, we derive the loglikelihood function conditional on the published studies at a fixed value of p.

We let $f_P\left (\hat {\mathbf y}^{(i)}\mid \hat {\boldsymbol \Sigma }^{(i)}\right )$ denote the marginal distribution of $\hat {\mathbf y}^{(i)}$ and $f_P\left (\hat {\boldsymbol \Sigma }^{(i)}\right )$ the distribution of $\hat {\boldsymbol \Sigma }^{(i)}$ over the population studies. In meta-analysis without taking into account reporting bias, the empirical data $\left (\hat {\mathbf y}^{(i)}, \hat {\boldsymbol \Sigma }^{(i)}\right )$ are regarded as random sample from the population studies with the joint distribution $f_P\left (\hat {\mathbf y}^{(i)}, \hat {\boldsymbol \Sigma }^{(i)}\right )$ . In the presence of selective publication, the published studies can be biased sample from the population; we let $f_O$ denote the distribution of the selectively published samples. Given a fixed p, the distribution of $\hat {\boldsymbol \Sigma }^{(i)}$ in the published studies is derived by

$$ \begin{align*} f_O\left(\hat{\boldsymbol\Sigma}^{(i)}\right) = f\left(\left.\hat{\boldsymbol\Sigma}^{(i)}\right|\text{select}\right) = \dfrac{P\left(\text{select}\mid\hat{\boldsymbol\Sigma}^{(i)}\right)f_P\left(\hat{\boldsymbol\Sigma}^{(i)}\right)}{P\left(\text{select}\right)} = \dfrac{b\left(\hat{\boldsymbol \Sigma}^{(i)}\right)f_P\left(\hat{\boldsymbol\Sigma}^{(i)}\right)}{p}, \end{align*} $$

which gives

(15) $$ \begin{align} f_P\left(\hat{\boldsymbol\Sigma}^{(i)}\right)= p\cdot \left\{b\left( \hat{\boldsymbol \Sigma}^{(i)}\right)\right\}^{-1}f_O\left(\hat{\boldsymbol\Sigma}^{(i)}\right). \end{align} $$

The joint distribution of the empirical data is then derived by

$$ \begin{align*} f_O\left(\hat{\mathbf y}^{(i)}, \hat{\boldsymbol{\Sigma}}^{(i)}\right) & = f\left(\hat{\mathbf y}^{(i)}, \hat{\boldsymbol{\Sigma}}^{(i)}\mid \text{select}\right) \\ & = \dfrac{P\left(\text{select}\mid \hat{\mathbf y}^{(i)}, \hat{\boldsymbol{\Sigma}}^{(i)} \right) f_P\left(\hat{\mathbf y}^{(i)}, \hat{\boldsymbol{\Sigma}}^{(i)}\right)}{p} \\ & = \dfrac {a\left({t_{\text{HR}}^{(i)}}\right) f_P\left(\hat{\mathbf y}^{(i)} \mid \hat{\boldsymbol{\Sigma}}^{(i)}\right) f_O\left(\hat{\boldsymbol{\Sigma}}^{(i)}\right)} {b \left(\hat{\boldsymbol \Sigma}^{(i)}\right) }. \end{align*} $$

This joint distribution allows us to derive the loglikelihood of published studies, that is,

(16) $$ \begin{align} \begin{aligned} \ell_O\left(\boldsymbol\theta, \boldsymbol{\Psi}, \alpha, \beta\right) &= \log\prod_{i=1}^N f_O\left(\hat{\boldsymbol{\mu}}^{(i)}, \hat{\boldsymbol{\Sigma}}^{(i)}\right) \\ &= \sum_{i=1}^{N}\log f_P\left(\hat{\boldsymbol{\mu}}^{(i)}\left| \hat{\boldsymbol{\Sigma}}^{(i)}\right.\right) + \sum_{i=1}^{N}\log a\left({t_{\text{HR}}^{(i)}}\right) - \sum_{i=1}^{N}\log b\left(\hat{\boldsymbol{\Sigma}}^{(i)}\right) + c, \end{aligned} \end{align} $$

where the second and the third terms are used for correcting reporting bias, and $c = \sum _{i=1}^{N}\log \left \{ {f_O\left (\hat {\boldsymbol {\Sigma }}^{(i)}\right )}\right \}$ is constant. If a random selection of studies holds, then either $\beta =0$ or $\beta \rightarrow \infty $ holds as well. Consequently, these two terms cancel each other, making the likelihood reduce to that without accounting for reporting bias.

Noting that by taking the integral of both sides in equation (15), we can derive

(17) $$ \begin{align} p = 1\left/E_O\left\{b\left(\hat{\boldsymbol{\Sigma}}^{(i)}\right)^{-1}\right\}\right. \simeq N\left/\sum_{i=1}^N b\left(\hat{\boldsymbol{\Sigma}}^{(i)}\right)^{-1}\right.. \end{align} $$

According to the definition of $b(\cdot )$ (13), equation (17) is monotonic with respect to $\alpha $ ; thus, the parameter $\alpha $ can be represented by the function of $(\boldsymbol \theta , \boldsymbol {\Psi }, \beta )$ given a value of p. We denote this by $\alpha _p = \alpha _p(\boldsymbol \theta , \boldsymbol {\Psi }, \beta )$ . By replacing the $\alpha $ with $\alpha _p$ , we derive the conditional loglikelihood given the published studies from the loglikelihood (16):

(18) $$ \begin{align} \begin{aligned} &\ell_O\left(\boldsymbol\theta, \boldsymbol{\Psi}, \beta; p\right) \\ \propto& -\dfrac{1}{2}\sum_{i=1}^{N}\left\{ \left(\hat{\mathbf y}^{(i)} - \boldsymbol\theta\right)^\top \left(\boldsymbol{\Psi}+\left.\hat{\boldsymbol{\Sigma}}^{(i)}\right/n^{(i)}\right)^{-1} \left(\hat{\mathbf y}^{(i)}- \boldsymbol{\theta}\right) + \log\left|\boldsymbol{\Psi}+\left.\hat{\boldsymbol{\Sigma}}^{(i)}\right/n^{(i)}\right|\right \} \\ +& \sum_{i=1}^{N}\log \Phi\left\{ \alpha_p + \beta\cdot \left({\hat \mu^{(i)}_{\text{lnHR}}}\left/{\hat s^{(i)}_{\text{lnHR}}}\right.\right) \right\} \\ -& \sum_{i=1}^{N}\log \Phi\left\{ \dfrac{\alpha_p + \beta\cdot \left(\theta_{\text{lnHR}}\left/\hat s^{(i)}_{\text{lnHR}}\right.\right) } {\sqrt{1 + \beta^2\cdot \left\{1+ \left({\psi_{\text{lnHR}}}\left/{\hat s^{(i)}_{\text{lnHR}}}\right)^2\right.\right\}}} \right\}. \end{aligned} \end{align} $$

Differently from the loglikelihood (16), with $\alpha _p$ , this loglikelihood is constraint by the marginal selection probability. The parameters $\left (\boldsymbol \theta , \boldsymbol {\Psi }, \beta \right )$ can be estimated by maximizing the conditional loglikelihood given a specified value of p. We denote these ML estimates as $(\hat {\boldsymbol \theta }, \hat {\boldsymbol \Psi }, \hat \beta )$ . With these MLEs, $\text {SROC}(t)$ in equation (5) and $\text {SAUC}(t)$ in equation (6) can be estimated accordingly, denoted by $\mathrm {SR}\hat {\mathrm O}\mathrm C(t) =\text {SROC}(x, t; \hat {\boldsymbol \theta }, \hat {\boldsymbol \Psi })$ and $\mathrm {SA}\hat {\mathrm U}\mathrm C(t) = \text {SAUC}(t; \hat {\boldsymbol \theta }, \hat {\boldsymbol \Psi })$ , respectively. The asymptotic normality of $(\hat {\boldsymbol \theta }, \hat {\boldsymbol \Psi }, \hat \beta )$ follows the general theory of the ML estimation under the assumptions that the number of studies and subjects are large. The asymptotic variance-covariance matrix of $(\hat {\boldsymbol \theta }, \hat {\boldsymbol \Psi }, \hat \beta )$ can be consistently estimated by the inverse of the empirical Fisher information following the ML theory. The two-tailed CI for the SAUC can be constructed by the delta method, and the detailed derivation is presented in Section S3 of the Supplementary Material.

In practice, the true value of p is unknown. As a sensitivity analysis for reporting bias, the value of p should be taken within the range $(0,1)$ in general. It is recommended to specify a decreasing series for p, such as $p=0.9,0.8,\dots , 0.1$ , to thoroughly examine the changes of $\mathrm {SR}\hat {\mathrm O}\mathrm C(t)$ or $\mathrm {SA}\hat {\mathrm U}\mathrm C(t)$ . Specifying a large value, such as $p=0.9$ , implies that the published studies account for 90% of the population, and thereby small reporting bias will be corrected. In contrast, specifying a small value, such as $p=0.1$ , indicates that the published studies represent 10% of the population, resulting in possibly substantial reporting bias to be corrected by the proposed method. If the range for the probability of publication can be reasonably conjectured, it is preferable to specify a decreasing series for p within that narrower range. In the following section, we illustrate the practical implementation of the proposed method.