Highlights
What is already known?
-
• Incidence proportions are often used in meta-analyses of adverse events.
-
• Heterogeneity in follow-up duration often occurs in meta-analyses of adverse events in randomized clinical trials, particularly in oncology, and this makes it difficult to interpret pooled incidence proportions using standard meta-analytic methods.
-
• Although survival analysis techniques would be useful for accounting for heterogeneous follow-up durations, individual studies rarely report summary measures based on survival analysis.
-
• No commonly accepted methods are available for literature-based meta-analyses of adverse events.
What is new?
-
• Novel methods are proposed to estimate cumulative incidence functions for adverse events using incidence proportions reported in individual studies, while accounting for heterogeneous follow-up durations.
Potential impact for RSM readers
-
• The proposed method enables fair comparisons of adverse event occurrence by appropriately adjusting for heterogeneity in follow-up duration across treatments.
-
• This approach would be beneficial for evaluating the risk–benefit balance of treatments.
1 Introduction
In clinical trials, it is often of primary interest to investigate efficacies of some drugs. Even if the primary endpoint is an efficacy parameter, safety assessment is a crucial aspect in all clinical trials. Details of all the adverse events (AEs) and their frequencies observed in the trial should be reported. For each AE, the incidence proportion, which is the number of subjects who experienced the AE divided by the number of subjects, is usually reported. Most randomized clinical trials are primarily designed to compare efficacies. Unless a randomized controlled clinical trial is prospectively designed with a specific safety outcome as its primary endpoint and sized accordingly, the trial may lack sufficient sample size to reliably assess whether there is an increased risk of AEs. In other words, a single clinical trial may have limited capacity to evaluate drug safety and then synthesizing information of AEs over multiple studies would provide valuable information.
Thus, meta-analysis would offer a valuable framework to understand safety profiles of drugs. In meta-analysis of AEs, the incidence proportion is often of primary interest, and the empirical incidence proportions are aggregated over studies and compared across treatment groups. For example, Zhu et al.Reference Zhu, Chen, Liu, Xiao and Liu 1 conducted a meta-analysis to examine whether bevacizumab increased cerebrovascular event risks. Cortes et al.Reference Cortes, Calvo and Ramirez-Merino 2 reported a meta-analysis of the AEs for breast cancer patients treated by bevacizumab and Gyawail et al.Reference Gyawail, Shimokawa, Ando, Honda and Ando 3 investigated serious and fatal AE risks of sorafenib with solid tumors. In the standard meta-analysis, the fixed-effect and the mixed-effect models are often used; the former assumes that the treatment effect is common, and the latter allows for heterogeneity over studies. To estimate the overall treatment effect, the inverse variance weighted average estimator, the normal–normal model, and the generalized linear mixed model are often used. In general, the incidence proportions of AEs are likely to be very low and may be zero for some studies. Thus, the standard methods for meta-analysis may not work in such situations with sparse data and several alternative methods have been developed taking account for such sparsity of the frequencies.Reference Cai, Parast and Ryan 4 – Reference Hong, Wang and Rosner 6
In addition to the sparsity issue, various issues specific to meta-analysis of AEs were discussed by a regulatory guidance document 7 and Cochrane Handbook (see Chapter 19 of Higgins and ThomasReference Higgins and Thomas 8 ). Heterogeneity in follow-up duration for evaluation of AEs is one of them (see Section III-D of FDA guidance 7 ). The duration of follow-up may differ among the treatment groups and thus should be factored into safety meta-analysis. Varying follow-up durations between the treatment groups due to discontinuing assigned drugs or dropouts may occur and impact the comparability among groups with respect to safety outcomes. If this is the case, aggregating the incidence proportions may be misleading and less informative.
When subject-level data are available, survival analysis techniques can be employed to account for differences in follow-up duration among studies. Stegherr et al.Reference Stegherr, Beyersmann and Jehl 9 , Reference Stegherr, Schmoor, Lübbert, Friede and Beyersmann 10 recommended use of survival analysis techniques in the analysis of AEs in clinical trials. However, such meticulous analyses are seldom executed in individual trials, and conducting individual participants data (IPD) meta-analysis is also challenging with difficulty to collect all the IPDs. Consequently, meta-analysis based on study-level data (literature-based meta-analysis) necessitates approaches that can accommodate disparities in follow-up duration. To our best knowledge, no formal methods have been proposed to handle heterogeneous follow-up duration in aggregated meta-analysis.
In clinical trials in oncology, the primary endpoint is often a time-to-event endpoint, such as the overall survival (OS) and the progression-free-survival (PFS). If the primary endpoint is the OS, the incidence proportion of the AEs is only for the period from the date of study initiation (typically the date of randomization) to the OS. Then, if the OSs are substantially longer in the experimental treatment than in the control treatment, it may be misleading to compare the incidence proportions between the treatment groups. This happens even with the PFS as the primary endpoint. Ethical care for subjects often leads to the treatment discontinuation and the subject is off the study once disease progression is observed. More exactly, in practice, with an extra follow-up duration of AEs such as one month, AEs are followed up to certain period after the PFS (e.g., PFS+1 month). If this is the case, interpreting the incidence proportion becomes difficult. Note that the follow-up scheme on the PFS is usually defined in the study protocol in oncology clinical trials. Thus, the issue of the heterogeneous follow-up duration is likely to arise common in oncology trials.
On the other hand, the statistical analysis for the time-to-event primary endpoint (e.g., PFS and OS) are usually conducted with the Kaplan–Meier method and since showing benefits in efficacy is of primary interest, the Kaplan–Meier estimates of the survival functions of the treatment groups are highlighted graphically in almost all the oncology studies. Kaplan–Meier estimates can provide valuable information regarding follow-up durations for AEs. In this article, we propose novel safety meta-analysis methods that address difference in follow-up duration for AEs among the treatment groups. Our idea is to analyze the AE as a time-to-event outcome by considering the cumulative incidence function as the estimand and propose estimation procedures for it utilizing the Kaplan–Meier estimates of the time-to-event primary endpoint, which form the distribution of the follow-up duration of AEs. At first, we introduce our proposed method assuming that the follow-up duration (e.g., PFS and OS) is independent of the time-to-AE. In practice, this independence assumption may be questionable. Thus, we also introduce an inference procedure accounting for dependence between the time-to-AE and the follow-up duration with copulas. As far as we know, our approach is the first safety meta-analysis method that can address difference in the follow-up duration.
The remainder of this article is structured as follows. Section 2 provides a motivating example of the meta-analysis of AEs in randomized clinical trials for bevacizumab and introduce the setting and notations. In Section 3, we introduce our proposed methods. We illustrate our proposed method by re-analyzing the meta-analysis of bevacizumab in Section 4 and examine the performance of the method by simulation studies in Section 5. Finally, we conclude this manuscript by describing advantages and limitations of the proposed method and potential future work in Section 6.
2 Preliminary
2.1 Motivating example: Meta-analysis of adverse events with bevacizumab
Our study was motivated by a meta-analysis of cerebrovascular events of anti-neoplastic agent bevacizumab.Reference Zhu, Chen, Liu, Xiao and Liu 1 This meta-analysis included 12,917 patients from 17 randomized clinical trials and the primary interest was to evaluate cerebrovascular events when administrated by bevacizumab in cancer patients. Basic characteristics of the studies are summarized in Table 1 of Zhu et al.Reference Zhu, Chen, Liu, Xiao and Liu 1 Studies of various cancer types were included with bevacizumab as an experimental treatment group. The primary outcome measure of the meta-analysis was the risk ratio (RR) of the bevacizumab group to the control, which was calculated as the ratio of the incidence proportion for occurrence of cerebrovascular events of the bevacizumab group to that of the control treatment. The RRs of the incidence proportion were aggregated by using the Mantel–Haenszel method over the studies. They reported the aggregated RR 3.28 (95%CI:1.97–5.48), concluding that bevacizumab had a significantly higher cerebrovascular events risk.
Summary of availability of information required by the proposed method. n: the number of subjects for safety analysis; MFT: median follow-up time (month); IP: incidence proportion (bev.: bevacizumab group and cnt.: control group); KM: availability of Kaplan–Meier estimates of the primary endpoint (open circle: reported and cross: missing); Follow-up: definition of the end of-follow-up for AEs in the protocol (Month); At-risk: availability of at-risk information for the primary endpoint typically attached below the figure for Kaplan–Meier estimates (open circle: reported and cross: missing); and NA: not available
Table 1 of the present article summarizes some additional characteristics of the enrolled studies, as well as important information from Table 1 of Zhu et al.Reference Zhu, Chen, Liu, Xiao and Liu 1 Among the 17 studies, 15 studies were designed to demonstrate efficacy of bevacizumab for a time-to-event endpoint. We searched descriptions on the definition of the end of follow-up in each individual publication. As presented in Table 1, some studies with the PFS as the primary endpoint discontinued the study follow-up after certain duration such as one month after the PFS was determined. In many studies, efficacy of bevacizumab was established with statistical significance. For example, Aghajanian et al.,Reference Aghajanian, Sill and Darcy 11 which is the study 1 in Table 1, conducted a doube blind, randomized, and placebo-controlled clinical trial for endometrial cancer patients of bevacizumab. The primary endpoint was the PFS and they reported a statistically significant prolongation of the PFS with bevacizumab; the median PFS was 12.4 month in bevacizumab group and that in the control was 8.4 month. In the study protocol of Aghajanian et al.,Reference Aghajanian, Sill and Darcy 11 the study was designed to discontinue follow-up for a patient at one month after a progression event. Therefore, the follow-up duration of the AE was longer in the bevacizumab group than the control group. Aghajanian et al.Reference Aghajanian, Sill and Darcy 11 reported the incidence proportion of the bevacizumab group as 2/247=0.008 and that of the control group as 1/233=0.004. Since the event rates were very low, it was hard to compare these two incidence proportions. Furthermore, heterogeneous follow-up durations between the groups complicated the comparison. Aggregation over studies may solve the issue of rare events, but cannot solve the heterogeneous follow-up issue. Then, it is unclear what the aggregated RR 3.28Reference Zhu, Chen, Liu, Xiao and Liu 1 estimated.
2.2 Setting and notation
Suppose we are interested in conducting a meta-analysis of AEs for randomized clinical trials in oncology. We consider a literature-based meta-analysis based on summary statistics (not IPD) reported in individual studies. Let the number of studies in the meta-analysis be denoted by S. Each individual study was designed to compare a time-to-event efficacy endpoint, such as the OS and the PFS between two treatment groups (experimental/control). Let k= 1 imply the experimental treatment and k=0 the control group. The number of subjects in the kth group of the sth study is denoted by
$n_{k}^{(s)}$
. To describe our proposed method, we need to prepare many time-to-event variables for individual participants, corresponding cumulative distribution functions and survival functions, and observed data. In Table 2, we summarize key notations. We begin with preparing notations for individuals in each study, which are not observed in the literature-based meta-analysis. As introduced in Section 2.1, the follow-up is discontinued once the primary time-to-event endpoint is confirmed; for the OS, the follow-up is until the date to die and for the PFS, after some extra duration defined in the study protocol, say one month, from the PFS, the follow-up is discontinued due to ethical consideration. We consider several random variables for time-to-events. We set the date of randomization as the origin of them. Let the time until the primary endpoint of the ith subject of the kth group in the sth study be denoted by
$T_{k,i}^{(s)}$
and it is right-censored by
$C_{k,i}^{(s)}$
, which is the end of the follow-up of the efficacy endpoint. Let
$\tilde {T}_{k,i}^{(s)}=min(T_{k,i}^{(s)}, C_{k,i}^{(s)})$
and
$\Delta _{k,i}^{(s)}=I(T_{k,i}^{(s)} \le C_{k,i}^{(s)})$
. The pair
$(\tilde {T}_{k,i}^{(s)}, \Delta _{k,i}^{(s)})$
is the standard data for survival analysis. We suppose the standard assumption for survival analysis;
$T_{k,i}^{(s)}$
and
$C_{k,i}^{(s)}$
are independent, which ensures validity of the Kaplan–Meier estimator. Let the survival and cumulative distribution functions of
$T_{k,i}^{(s)}$
be denoted by
$S_{k}^{(s)}(t)=P(T_{k,i}^{(s)}>t)$
and
$F_{k}^{(s)}(t)=1-S_{k}^{(s)}(t)$
, respectively. The Kaplan–Meier estimate for
$S_{k}^{(s)}(t)$
based on the observations
$(\tilde {T}_{k,i}^{(s)}, \Delta _{k,i}^{(s)})$
(
$i=1,2,\ldots , n_k^{(s)}$
) is denoted by
$\hat {S}_{k}^{(s)}(t)$
. Denote
$G_{k}^{(s)}(t)=P(C_{k,i}^{(s)}>t)$
. Let the time until the end of follow-up be denoted by
$T_{FU:k,i}^{(s)}$
; for the OS,
$T_{FU: k,i}^{(s)}$
agrees with
$T_{k,i}^{(s)}$
and for the PFS, extra follow-up duration, say one month, is added according to the procedure defined in the study protocol. The survival and cumulative distribution functions for
$T_{FU: k,i}^{(s)}$
are denoted by
$S_{FU: k}^{(s)}(t)$
and
$F_{FU: k}^{(s)}(t)$
, respectively. Let
$\tilde {T}_{FU: k,i}^{(s)}=min(T_{FU: k,i}^{(s)}, C_{k,i}^{(s)})$
and its survival and cumulative distribution functions be denoted by
$\tilde {S}_{FU: k}^{(s)}(t)$
and
$\tilde {F}_{FU: k}^{(s)}(t)$
, respectively. Furthermore, we define the time-to-AE,
$T_{AE:k,i}^{(s)}$
by the duration of the day of the first occurrence of the AE from the day of randomization. We assume that the distribution of
$T_{AE:k,i}^{(s)}$
is common over the studies, and then denote its cumulative distribution function and its survival function by
$F_{AE:k}(t)$
and
$S_{AE:k}(t)$
, respectively, dropping the superscript for dependence to the study.
List of notations
Recall that we focus on literature-based meta-analysis, in which only information extracted from published papers is available. Thus, individual participant data are not available. We summarize data, which are supposed to be available in conducting the proposed method. Let
$X_k^{(s)}$
be the number of subjects in the kth groups of the sth study, who experienced AEs during the study. We suppose
$n_k^{(s)}$
and
$X_k^{(s)}$
are available for all the studies. In addition, we assume that all the studies report the Kaplan–Meier estimate
$\hat {S}_{k}^{(s)}(t)$
of the two treatment groups graphically, whose values are extracted from the graphical plot using some graph scan software such as the Digitizelt software (https://www.digitizeit.xyz/ja/). With them, we suppose that the survival function of the time to the end of the follow-up for AEs,
$S_{FU: k}^{(s)}(t)$
are well estimated. Specifically, for the OS,
$T_{FU: k,i}^{(s)}=T_{k,i}^{(s)}$
and then
$\hat {S}_{k}^{(s)}(t)$
is used. For PFS, some additional follow-up duration
$\tau ^{(s)}$
(say, one month) is usually defined in the protocol and then
$S_{FU: k}^{(s)}(t)$
is estimated by
$\hat {S}_{k}^{(s)}(t-\tau ^{(s)})$
. The resulting estimator for
$S_{FU: k}^{(s)}(t)$
is denoted by
$\hat {S}_{FU: k}^{(s)}(t)$
;
$\hat {S}_{FU: k}^{(s)}(t)=\hat {S}_{k}^{(s)}(t)$
for the OS and
$\hat {S}_{FU: k}^{(s)}(t)=\hat {S}_{k}^{(s)}(t-\tau ^{(s)})$
for the PFS. Let
$\hat {F}_{FU: k}^{(s)}(t)=1-\hat {S}_{FU: k}^{(s)}(t)$
. To estimate
$G_{k}^{(s)}(t)$
, the survival function of
$C_{k,i}^{(s)}$
, we suppose that for each study, the number of subjects at risk for the efficacy endpoint is given for some time points. Let
$a_k^{(s)}(t)$
be the number of subjects at-risk at t. Typically, it is given below the graphical plot of the Kaplan–Meier estimates for the efficacy endpoint. It is assumed that for the kth group of the sth study,
$a_k^{(s)}(t)$
are available at
$t=t_1^{(s)}, t_2^{(s)}, \ldots , t_{K^{(s)}}^{(s)}$
, where
$K^{(s)}$
is the number of time points for the sth study. For studies without
$a_k^{(s)}(t)$
available, we suppose that the median follow-up time
$m^{(s)}$
is available, which is defined as the sample median of
$\tilde {T}_{k,i}^{(s)}=min(T_{k,i}^{(s)}, C_{k,i}^{(s)})$
over all the subjects of the two groups in each individual study and sometimes is reported in literature to describe the status of the follow-up of the time-to-event endpoint as presented in Table 1.
3 Proposed methods
3.1 Methods under the independent follow-up duration
The idea of our proposed method is to summarize the occurrence of AEs by estimating the cumulative distribution function
$F_{AE: k}(t)$
with information of the follow-up duration from the Kaplan–Meier estimates of the efficacy outcomes as the end of the follow-up of AEs. We begin with introducing our proposed method under a simpler setting of the independent follow-up duration. We assume that
$\tilde {T}_{FU: k,i}^{(s)}$
and
$T_{AE: k,i}^{(s)}$
are independent. Recall that
$\tilde {T}_{FU: k,i}^{(s)}=min(T_{FU: k,i}^{(s)}, C_{k,i}^{(s)})$
. The event that an AE is observed is equivalent to
$T_{AE:k,i}^{(s)}$
occurring before the end of follow-up,
$\tilde {T}_{FU: k,i}^{(s)}$
. Thus, the incidence proportion in each individual study is represented as
$p_k^{(s)}=P \left ( T_{AE: k,i}^{(s)} \le \tilde {T}_{FU: k,i}^{(s)} \right )$
. By simple algebraic manipulation, with the assumption
$\tilde {T}_{FU: k,i}^{(s)} \perp T_{AE: k,i}^{(s)}$
, it holds that
$$ \begin{align} p_k^{(s)}&= \int_0^\infty P \left(T^{(s)}_{AE: k,i} \le t \mid \tilde{T}^{(s)}_{FU: k,i}=t \right) d\tilde{F}^{(s)}_{FU: k} (t) \nonumber \\ &=\int_0^\infty P\left(T_{AE: k,i}^{(s)} \le t \right) d\tilde{F}_{FU: k}^{(s)} (t) \nonumber \\ &=\int_0^\infty F_{AE: k}(t) d\left(1-S_{FU: k}^{(s)}(t) G_{k}^{(s)}(t) \right). \end{align} $$
The cumulative distribution function
$F_{AE: k}(t)$
is a key quantity to summarize the occurrence of AEs. Recall that we assume that it is common over the studies. We suppose some parametric model for
$F_{AE: k}(t)$
, such as the exponential, the Weibull, and the log-normal distributions. The cumulative distribution function of the model is denoted by
$F_{AE: k}(t; \theta _k)$
, where
$\theta _k$
is the unknown parameters. We estimate the censoring distribution
$G_{k}^{(s)}(t)$
by assuming some parametric model
$ G_{k}^{(s)} (t; \eta _k)$
. The unknown parameter
$\eta _k$
is estimated by the method given in Section 3.2, and the resulting estimator is denoted by
$\hat {\eta }_k$
. By replacing
$S_{FU: k}^{(s)}(t)$
and
$ G_{k}^{(s)}(t)$
in (1) with
$\hat {S}_{FU: k}^{(s)}(t)$
and
$ G_{k}^{(s)} \left (t; \hat {\eta }_k \right )$
, respectively, an empirical version of
$p_k^{(s)}$
is given by
$$ \begin{align} p_k^{(s)}(\theta_k) &=\int_0^\infty F_{AE: k}(t; \theta_k) d\left(1-\hat{S}_{FU: k}^{(s)}(t) G_{k}^{(s)} \left(t; \hat{\eta}_k\right) \right). \end{align} $$
Here,
$\left (1-\hat {S}_{FU: k}^{(s)}(t) G_{k}^{(s)} \left (t; \hat {\eta }_k\right ) \right )$
is a cumulative distribution function of a discrete distribution with masses at the jump points of
$\hat {S}_{FU: k}^{(s)}(t)$
and the integral on the right-hand side of (2) is defined by the Riemann (Lebesgue)–Stieltjes integral. To be more specific, let the number of distinct time points, at which failure times of some subjects were observed, be denoted by
$H^{(s)}$
, and the ordered observed failure times of the sth study be denoted by
$0=w_0^{(s)} \le w_1^{(s)} \le \dots w_{H^{(s)}}^{(s)}$
. Then, equation (2) is calculated by
$$ \begin{align*} p_k^{(s)}(\theta_k) &=\sum_{k=1}^{H^{(s)}} {F_{AE: k}(w_{k-1}^{(s)}; \theta_k)} \Big\{ \hat{S}_{FU: k}^{(s)}(w_{k-1}^{(s)}) G_{k}^{(s)} \left(w_{k-1}^{(s)}; \hat{\eta}_k\right)- \hat{S}_{FU: k}^{(s)}(w_{k}^{(s)}) G_{k}^{(s)} \left(x_{k}^{(s)}; \hat{\eta}_k\right) \Big\}. \end{align*} $$
Assuming
$X_k^{(s)} \sim $
Binomial(
$n_k^{(s)}$
,
$p^{(s)}_k(\theta _k)$
), we propose to estimate the unknown parameter
$\theta _k$
by applying the maximum likelihood method. The likelihood function is defined by
$$ \begin{align*} &\prod_{s=1}^S {n^{(s)}_k \choose X^{(s)}_k} \{p_k^{(s)}(\theta_k)\}^{X^{(s)}_k}\left\{1- p_k^{(s)}(\theta_k) \right\}^{n^{(s)}_k-X^{(s)}_k}. \end{align*} $$
Using the maximum likelihood estimator
$\widehat {\theta _k}$
, the cumulative distribution function of the occurrence of the AE before time t is estimated by
$ F_{AE: k}(t; \hat {\theta }_k)$
and its
$95\%$
CIs is constructed as
$$ \begin{align*} F_{AE: k}(t; \hat{\theta}_k) \pm1.96 \times \frac{\sqrt {Var \left( F_{AE: k}(t; \hat{\theta}_k)\right) } }{\sqrt {S}}, \end{align*} $$
where the variance formula
$Var \left ( F_{AE: k}(t; \theta _k) \right ) $
is derived by using the delta method with the variance estimates of the maximum likelihood estimator.
3.2 Estimation of the censoring distribution
Recall that
$a_k^{(s)}(t)$
is the number of subjects at-risk for the efficacy endpoint at t, which is available at
$t=t_1^{(s)}, t_2^{(s)}, \ldots , t_{K^{(s)}}^{(s)}$
. The proportion of subjects at-risk of the AE at time t is denoted by
$q_k^{(s)}(t)=a_k^{(s)}(t)/n_k^{(s)}$
. We assume that
$C_{k,i}^{(s)}$
follows the exponential distribution with the parameter
$\eta _k$
and its survival function is denoted by
$G_k^{(s)}(t; \eta _k)$
. Under the assumption of independence between
$T_{k,i}^{(s)}$
and
$C_{k,i}^{(s)}$
,
$q_k^{(s)}(t)=S_k^{(s)}(t) G_{k}^{(s)} \left (t; \eta _k \right )$
holds. It motivates us to estimate the unknown parameter
$\eta _k$
by applying the least square method for each study separately:
$$ \begin{align*} minimize \sum \left\{ q_k^{(s)}(t)- \widehat{S}_{k}^{(s)}(t)G_{k}^{(s)} \left(t; \eta_k \right) \right\}^2, \end{align*} $$
where the summation is taken over
$t=t_1^{(s)}, t_2^{(s)}, \ldots , t_{K^{(s)}}^{(s)}$
.
For studies without
$a_k^{(s)}(t)$
reported, the following approach can be taken. Suppose the median follow-up time
$m^{(s)}$
is reported for the primary endpoint, which is defined by a simple sample median of
$\tilde {T}_{k,i}^{(s)}=min (T_{k,i}^{(s)}, C_{k,i}^{(s)})$
over the subjects of the two treatment groups for the sth study. It is often done as presented in Table 1. Note that the median follow-up time is reported with various definitions.Reference Betensky
12
Suppose the censoring distribution is common over the two treatment groups, and the survival function of the censoring is denoted by
$G^{(s)}(t)$
. We suppose a parametric model
$G^{(s)}(t; \eta )$
for
$G(t)$
. Since we assume that the independence between
$T_{k,i}^{(s)}$
and
$C_{k,i}^{(s)}$
,
$P(\tilde {T}_{k,i}^{(s)}>t)=S^{(s)}(t) G^{(s)}(t; \eta )$
holds, where
$S^{(s)}(t)$
is the marginal survival function of the sth study and can be estimated by
$\hat {S}^{(s)}(t)=\frac {n_0^{(s)}}{n_1^{(s)}+n_0^{(s)}} \hat {S}_0^{(s)}(t)+\frac {n_1^{(s)}}{n_1^{(s)}+n_0^{(s)}} \hat {S}_1^{(s)}(t)$
. Then, since we assume that
$C_{k,i}^{(s)}$
follows the exponential distribution with the hazard parameter
$\eta $
, we can determine
$\eta $
by solving
$0.5=\hat {S}^{(s)}(m^{(s)}) exp(-\eta m^{(s)})$
.
In the above argument, we assumed that
$C_{k,i}^{(s)}$
follows the exponential distribution, which had a single parameter. For the studies with
$a_k^{(s)}(t)$
observed at multiple time points, distributions with multi-dimensional parameters can be handled in a similar way.
3.3 Methods under the dependent follow-up duration
The independent follow-up assumption
$\tilde {T}_{FU: k,i}^{(s)} \perp T_{AE: k,i}^{(s)}$
in Section 3.1 may not hold in practice. This section considers an approach not requiring the independence assumption. By simple algebraic manipulation, the incidence proportion
$\tilde {p}_k^{(s)}=P \left ( T_{AE: k,i}^{(s)} \le \tilde {T}_{FU: k,i}^{(s)} \right )$
is represented as
$$ \begin{align} \tilde{p}_k^{(s)}= \int_0^\infty \frac{P\left(T^{(s)}_{AE: k,i} \le t, \tilde{T}^{(s)}_{FU: k,i}=t \right)}{P \left(\tilde{T}^{(s)}_{FU: k,i}=t \right)} d\left(1-S_{FU: k}^{(s)}(t) G_{k}^{(s)}(t) \right). \end{align} $$
We propose to model the joint probability
$P\left (T_{AE: k,i}^{(s)} \le t_1, \tilde {T}^{(s)}_{FU: k,i} \le t_2 \right )$
in (3) using a two-dimensional parametric copula, such as Frank copula,Reference Frank
13
Clayton–Cook–Johnson copula,Reference Clayton
14
,
Reference Cook and Johnson
15
and Gumbel–Hougaard copula.Reference Gumbel
16
Recall that
$F_{AE: k}{(t)}=P\left (T_{AE: k,i}^{(s)} \le t \right )$
and denote
$\tilde {F}_{FU: k}^{(s)} (t)=P\left (\tilde {T}^{(s)}_{FU: k,i} \le t \right )$
. Thus, we model the joint cumulative distribution function by
$P\left (T^{(s)}_{AE: k,i} \le t_1, \tilde {T}^{(s)}_{FU: k,i} \le t_2 \right )=C_{k,\psi _k}(F_{AE: k}(t_1), \tilde {F}_{FU: k}^{(s)}(t_2))$
, where
$C_{k,\psi _k}(u,v)$
is a parametric copula model for the kth group,
$\psi _k$
is an unknown parameter in the copula;
$C_{k,\psi _k}(u, v)=\psi _k^{-1} log\{1+\frac {(\exp {(-\psi _k u)-1})(\exp {(-\psi _k v)-1})}{\exp {(-\psi _k u)+1}}\}$
with
$\psi _k \ne 0$
for Frank copula,
$C_{\psi _k}(u, v)=(u^{-\psi _k}+v^{-\psi _k})$
with
$-1 \le \psi _k, \psi _k \ne 0$
for Clayton–Cook–Johnson copula, and
$C_{k,\psi _k}(u, v)=exp\{-[(-log(u))^{\psi }+(-log(v))^{\psi _k}]^{1/\psi _k}\}$
with
$-1 \le \psi _k$
for Gumbel–Hougaard copula. By simple algebra, it holds that
$$ \begin{align*} P\left(T^{(s)}_{AE: k,i} \le t_1, \tilde{T}^{(s)}_{FU: k,i} = t_2 \right)= \frac{d}{dv} C_{k,\psi_k}(F_{AE: k}(t_1), v)|_{v=\tilde{F}_{FU: k}^{(s)}(t_2)} \times \frac{d}{dt_2} \tilde{F}_{FU: k}^{(s)}(t_2). \end{align*} $$
Then (3) is expressed as
$$ \begin{align} \tilde{p}_k^{(s)}=\int_0^\infty \frac{\partial}{\partial u} C_{k,\psi_k} \left(F_{AE: k}(t), v \right)|_{v=\left(1-S_{FU:k}^{(s)}(t) G_{k}^{(s)}(t) \right)} d\left(1-S_{FU:k}^{(s)}(t) G_{k}^{(s)}(t) \right). \end{align} $$
As done in the case of the independent follow-up, by replacing
$S_k^{(s)}(t)$
and
$G_{k}^{(s)}(t)$
in (4) with
$\hat {S}_k^{(s)}(t)$
and
$G_{k}^{(s)}(t; \hat {\eta }_k)$
, respectively, an empirical version of
$\tilde {p}_k^{(s)}$
is given by
$$ \begin{align*} \tilde{p}_k^{(s)}(\theta_k, \psi_k)=\int_0^\infty \frac{\partial}{\partial u}C_{k,\psi_k} \left(F_{AE: k}(t; \theta_k), v \right)|_{v=(1-\hat{S}_{FU: k}^{(s)}(t)G_{k}^{(s)}(t; \hat{\eta}_k))} d\left(1-\hat{S}_k^{(s)}(t) G_{k}^{(s)}(t; \hat{\eta}_k) \right). \end{align*} $$
Assuming
$X_k^{(s)} \sim $
Binomial(
$n_k^{(s)}$
,
$\widetilde {p}^{(s)}_k(\theta _k, \psi _k)$
), one can estimate unknown parameters
$\theta _k$
and
$\psi _k$
and conduct statistical inference following the maximum likelihood method.
To obtain the maximum likelihood estimators in Sections 3.1 and 3.3, since
$\tilde {p}_k^{(s)}(\theta _k)$
and
$\tilde {p}_k^{(s)}(\theta _k, \psi _k)$
are a non-linear function with respect to unknown parameters, we need to apply non-linear optimization. We used NLMIXED procedure in SAS (SAS Institute Co. Ltd.). An SAS code to implement the proposed method is available at 10.5281/zenodo.16872253.
4 Example: Re-analysis of the motivating meta-analysis
We illustrate our proposed method by re-analyzing the motivating meta-analysis dataset in Section 2.1. As presented in Table 1, among the 17 studies in the meta-analysis, 12 studies had the PFS or OS as the primary endpoint and reported the Kaplan–Meier estimates of its survival function for the two treatment groups as a graphical display. Among them, the study 11 reported neither the number of at-risk nor the median follow-up time for the primary endpoint. We exclude the study 11 in our illustration. The studies 3 and 5 were designed as a study with the two bevacizumab groups and a control one. In this illustration, we included the two bevacizumab groups of these studies separately in the model and regarded them as independent studies for illustration. Then, there were 13 studies for the both groups. For these studies, to determine the censoring distribution, we suppose the exponential distribution for
$C_{k,i}^{(s)}$
. The unknown parameter
$\eta _k$
was estimated by the least-square method for the number of subjects at-risk in Section 3.2 if it was available. For the studies without the number of subjects at-risk, we utilized the median follow-up time to this purpose by the method in Section 3.2.
At first, we apply the proposed method under the independence follow-up assumption. We considered three candidate models for
$T_{AE: k,i}^{(s)}$
; we assumed it followed the exponential, the Weibull, or the log-normal distributions. To select the most relevant model, we calculated the Akaike information criterion (AIC) for each model;
$AIC$
=-2
$loglikelihood$
+2p, where p is the number of unknown parameters to estimate.Reference Akaike
17
The model with smaller AIC is regarded as fitting better to data. The AICs of the candidate distributions are shown in Table 3. For the bevacizumab group, the log-normal distribution had the smallest AIC, suggesting that the bevacizumab group exhibited a better fit when assuming it, while the control group an exponential distribution. To further evaluate the goodness-of-fit of the selected model, we calculated and contrasted the model-based incidence proportion with the crude one. The model-based incidence proportion is defined as
$$ \begin{align*} p_k^{(s)*} = \int_0^\infty F_{AE: k}(t;\hat{\theta}_k) d\left(1-\widehat{S}_{FU: k}^{(s)}(t) G_{k}^{(s)} \left(t; \hat{\eta}_k \right) \right) \end{align*} $$
and the crude one as
$\hat {p}_k^{(s)}=X_k^{(s)}/n_k^{(s)}$
. In Figure 1, pairs of (
$\hat {p}_k^{(s)}$
,
$p_k^{(s)*}$
) are plotted, for the bevacizumab (left panel) and the control (right panel) groups, respectively. The plots should lie around the diagonal line when the model successfully captures the true structure. To make a formal evaluation of the goodness-of-fit, one may use chi-squared goodness-of-fit test. Let
$$ \begin{align} Z_{GOF, k}^2=\sum_{s=1}^S \frac{(X_k^{(s)}-n_k^{(s)} p_k^{(s)*})^2}{n_k^{(s)} p_k^{(s)*}}. \end{align} $$
Note that
$n_k^{(s)} p_k^{(s)*}$
is the expected number of the occurrence of the AE if the fitted model is correctly specified. Under the null hypothesis that the model is correctly specified,
$Z_{GOF, k}^2$
follows the chi-squared distribution of the degree of freedom
$S-1$
. It attaches a quantitative measure of the goodness-of-fit. Figure 1 suggests that the bevacizumab group generally showed a good fit, while the control group exhibited slightly poorer. On the other hand, the p-value of the goodness-of-fit test for the bevacizumab and the control groups were p = 0.021 and p = 0.012, respectively. Thus, the goodness-of-fit was not satisfactory in the both groups. The observed poor fit in the control group can be attributed to the fact that there were several trials in the control group (9 out of 13 trials) with no AEs reported.
Akaike information criterion (AIC) for comparison of models for the distribution of the time-to-AE and the copula for dependent follow-up in the meta-analysis of the motivating data; NA, not available due to non-convergence of non-linear optimization to maximize the log-likelihood function
Goodness-of-fit of the best model under the independent follow-up duration for the bevacizumab group (S = 13, left panel) and the control group (S = 13, right panel); p-value is for the Pearson chi-squared goodness-of-fit test.
Figure 2 depicts the cumulative incidence proportion of AEs,
$F_{AE: k}(t)$
, (left panel) and the corresponding hazard function (right panel) when fitting the distribution with the best AIC for both groups. The cumulative incidence proportion was higher uniformly over time in the bevacizumab groups. The exponential distribution was selected for the control group, and then the hazard of the occurrence of the AE was constant over time. On the other hand, the hazard of AEs increased in the bevacizumab group (see the right panel of Figure 2).
Estimated cumulative incidence functions
$F_{AE:k}(t; \hat {\theta }_k)$
(solid curve) with point-wise two-sided 95
$\%$
confidence bands (broken curve) of the bevacizumab and the control groups (left panel) and the corresponding hazard functions (right panel) under the independent follow-up duration.
Next, we conducted analyses under the dependent follow-up duration. For the copula function
$C_{k,\psi _k}(u, v)$
, we consider Frank copula, Clayton–Cook–Jonhson copula, and Gumbel–Hougaard copula as candidate models. In Table 3, we report the AICs for the candidate models. For the bevacizumab group, the best fit was observed when assuming the log-normal distribution for
$F_{AE: k}(t)$
and the Gumbel–Hougaard family for the copula function. For the control group, the best fit was found when assuming the exponential distribution and the Gumbel–Hougaard family. As done under the independent follow-up case, the goodness-of-fit was assessed by plotting the model-based incidence proportion and the crude one. The model-based incidence proportion was defined as
$$ \begin{align*} \tilde{p}_k^{(s)*}=\int_0^\infty \frac{\partial}{\partial v}C_{k, \hat{\psi}_k} \left(F_{AE: k}(t;\hat{\theta}_k), v \right)\mid _{v=\left(1-\hat{S}_k^{(s)}(t) G_{k}^{(s)}\left(t; \widehat{\eta}_k \right)\right)} d\left(1-\hat{S}_{FU: k}^{(s)}(t) G_{k}^{(s)}\left(t; \widehat{\eta}_k \right)\right), \end{align*} $$
and the plot is given in Figure 3 for the bevacizumab group (left panel) and the control group (right panel), respectively. A similar goodness-of-fit test can be constructed by replacing
$p_k^{(s)*}$
with
$\tilde {p}_k^{(s)*}$
in (5). For the bevacizumab group, the best model under the dependent follow-up duration had a slightly larger AIC value than the best model under the independent follow-up duration. Correspondingly, Figure 3 suggested a similar goodness-of-fit with a significance p-value of 0.030 for the Pearson chi-squared goodness-of-fit test. For the control group, AIC was also not improved accounting for dependence, and Figure 3 suggested unsatisfactory fit to the data with the p-value 0.013 of the goodness-of-fit test. Although these observations necessitate further improvement of the models, we compare the occurrence of AEs between the two treatment groups using the best-performing models. In Figure 4, we present the cumulative incidence functions of the two groups (left panel) and the corresponding hazard functions (right panel). The confidence interval of the cumulative incidence function was much wider under the dependent follow-up duration. In the control group, some studies did not report any occurrence of the AE. Thus, the model under dependent follow-up duration had extra parameters in the copula model and then estimation was instable with limited information. In Table 4, we present the RR of the bevacizumab group to the control with respect to the cumulative incidence up to selected time points with the best models of each group. As suggested by Figures 2 and 4, the hazard of the AE in the bevacizumab increased over time. Correspondingly, significance differences of the cumulative incidence were observed consistently after t = 30. Recall that the aggregated RR by the Mantel–Haenszel method was 3.28 (95
$\%$
CI: 1.97–5.48). Table 4 indicates that the RR was not homogeneous over time.
Goodness-of-fit of the best model under the dependent follow-up duration for the bevacizumab group (S = 13, left panel) and the control group (S = 13, right panel); p-value is for the Pearson chi-squared goodness-of-fit test.
Estimated cumulative incidence functions
$\ F_{AE:k}(t; \hat {\theta }_k)$
(solid curve) with point-wise two-sided 95
$\%$
confidence bands (broken curve) of the bevacizumab and the control groups (left panel) and the corresponding hazard functions (right panel) under the dependent follow-up duration.
The relative risk of the cumulative occurrence of the AE until selected time points by the model minimizing the AIC under the independent and the dependent follow-up durations
5 Simulation study
5.1 Data generation and applied methods
To evaluate the validity of the method proposed in Section 3, we conducted a simulation study. The setting for the simulation study was designed mimicking the dataset of Zhu et al.Reference Zhu, Chen, Liu, Xiao and Liu
1
We only considered the bevacizumab group (only
$k=1$
) since our proposal was to fit the model separately by the group. We selected 13 studies with the median time of the primary time-to-event endpoint reported. The median survival times for the selected 13 studies are listed Table 5. In generating simulation dataset, although the median survival time was for the OS in some studies, we regarded it as the time for the PFS for simplicity and assume that the PFS was the primary endpoint and the follow-up of the AEs were made until the PFS. We suppose that the Kaplan–Meier estimate for the PFS was available and the at-risk data were available at three time points (6, 12, and 18 months) for all the studies. Two setting on the follow-up duration for AEs were considered: (1) the independent follow-up duration (Dataset 1) and (2) the dependent follow-up duration (Dataset 2). For each setting, we generated 1,000 meta-analyses of 13 studies, whose sample size of each group was the same as the bevacizumab group of each study in Table 1 (see Table 5).
List of parameters of studies used in generating simulation datasets: n shows the number of subjects for efficacy analysis in the bevacizumab group;
$\lambda ^{(s)}$
: the hazard of the exponential distribution for PFS; and
$\lambda _C^{(s)}$
: the hazard of the exponential distribution for censoring
Dataset 1A: Independence
We assumed that
$T_{AE: k,i}^{(s)}$
followed an exponential distribution with a hazard of
$\lambda _{AE}$
. Three settings of
$\lambda _{AE}$
were considered;
$\lambda _{AE}=0.002$
(scenario 1),
$=0.01$
(scenario 2), and
$=0.05$
(scenario 3). The distribution of
$T_{k,i}^{(s)}$
was supposed as an exponential distribution, whose hazard
$\lambda ^{(s)}$
was set based on the median PFS in Table 5. The censoring for the PFS,
$C_{k,i}^{(s)}$
was generated from an exponential distribution, whose hazard
$\lambda _C^{(s)}$
was set a study-specific value listed in Table 5. Then, we set
$\tilde {T}_{k,i}^{(s)}=min(T_{k,i}^{(s)}, C_{k,i}^{(s)})$
and
$\Delta _{k,i}^{(s)}=I(T_{k,i}^{(s)} \le C_{k,i}^{(s)})$
and the occurrence of AE was determined by whether
$T_{AE: k,i}^{(s)} \leq \tilde {T}_{k,i}^{(s)}$
met or not. With these individual-level simulated data, we calculate the Kaplan–Meier estimate for the PFS and the number of AEs observed within the follow-up duration, which is denoted as
$X_k^{(s)}$
. The number of at-risk for the PFS,
$a_k^{(s)}(t)$
was calculated at
$t=6,12,18$
. We regard them as data.
Dataset 2A: Dependence
We assume the same three scenarios for
$\lambda _{AE}$
as Dataset 1. For the marginal distribution of
$T_{k, i}^{(s)}$
and
$C_{k,i}^{(s)}$
, the same exponential parameters
$\lambda ^{(s)}$
and
$\lambda _C^{(s)}$
were considered, respectively, as given in the case of the independent follow-up duration (see Table 5). From the independence assumption between
$T_{k,i}^{(s)}$
and
$C_{k,i}^{(s)}$
, the cumulative distribution function of
$\tilde {T}_{FU: k,i}^{(s)}=min(T_{k,i}^{(s)}, C_{k,i}^{(s)})$
is given by
$\tilde {F}_{FU,k}^{(s)}=1-exp(-t/(\lambda ^{(s)}+\lambda _C^{(s)}))$
. We assumed that the dependence between
$T_{AE: k,i}^{(s)}$
and
$\tilde {T}_{FU: k,i}^{(s)}$
was determined by the copula function. We considered three copulas: Clayton–Cook–Johnson, Gumbel–Hougaard, and Frank copulas. We set the parameter
$\psi _K$
in each copula so that the resulting Kendall’s
$\tau $
values were 0.2 or 0.5. We generated
$T_{AE: k,i}^{(s)}$
and
$\tilde {T}_{FU: k,i}^{(s)}$
from the joint distribution defined by the assumed copula with the supposed marginal distributions of
$T_{AE: k,i}^{(s)}$
and
$\tilde {T}_{FU: k,i}^{(s)}$
according to the method given in Tsukahara.Reference Tsukahara
18
We set
$\tilde {T}_{k,i}^{(s)}=\tilde {T}_{FU: k,i}^{(s)}$
. By setting
$\Delta _{k,i}^{(s)}=1$
with the probability
$P(T_{k,i}^{(s)} \le C_{k,i}^{(s)})$
and otherwise
$\Delta _{k,i}^{(s)}=0$
, we created the observation (
$\tilde {T}_{k,i}^{(s)}$
,
$\Delta _{k,i}^{(s)}$
) for the PFS. The occurrence of the AE was also determined according to
$T_{AE: k,i}^{(s)} \leq \tilde {T}_{k,i}^{(s)}$
. Then, the marginal survival function for
$T_{k,i}^{(s)}$
was estimated by the Kaplan–Meier method and the number of AEs
$X_{k}^{(s)}$
was calculated. The number of at-risk for the PFS,
$a_k^{(s)}(t)$
was calculated at
$t=6, 12,18$
.
We conducted an additional simulation study with S = 30. The number of subjects in each study was generated from the discrete uniform distribution on [50, 500]. The other settings were not altered. We generated datasets under the independent and dependent follow-up durations and the corresponding datasets were called Datasets 1B and 2B, respectively.
In the setting explained above, the at-risk data were available at 6, 12, and 18 months for all the studies. It is called Case 1. To evaluate impacts of availability of the at-risk data, in addition to Case 1, we consider an alternative setting, in which the at-risk data were available at two time points randomly selected on
$\{2, 3, \ldots ,15\}$
for each study. It is called Case 2. We evaluate performance of the method in Section 3.2 in estimating the exponential parameter for the study-specific censoring distribution for Cases 1 2.
5.2 Results
We analyzed both Datasets 1A and 2A using the two approaches: the method in Section 3.1 (independent follow-up duration) and that in Section 3.3 (dependent follow-up duration). In both Datasets 1A and 2A, the marginal distribution for the time until the first occurrence of the AE
$T_{AE: k,i}^{(s)}$
was the exponential distribution with the hazard
$\lambda _{AE}$
. For each of the three scenarios of
$\lambda _{AE}$
(scenario 1,
$\lambda _{AE}=0.002$
; scenario 2,
$\lambda _{AE}=0.01$
; and scenario 3;
$\lambda _{AE}=0.05$
), we fit the exponential distribution and evaluated bias and coverage probability (CP) on estimation of
$\lambda _{AE}$
. Recall that the target estimand of our proposed method is the cumulative distribution function
$F_{AE:k}(t)$
. Since
$T_{AE: k,i}^{(s)}$
was generated from the exponential distribution, it is given by
$F_{AE:k}(t)=1-exp(-\lambda _{AE} t)$
. Then, we evaluate estimates for
$\lambda _{AE}$
. To see how likely the AE occurred in each scenario, the average of the observed incidence proportion was calculated over 1,000 simulated datasets for Dataset 1; they were 0.023, 0.116, and 0.556 for scenarios 1–3, respectively.
Summary of estimates (median (q1, q3)) and empirical coverage probabilities (CP) of the exponential parameter for the time-to-AE and frequencies that each model was selected with AIC for Dataset 1A (S = 13, independence):
$True$
implies the value of
$\lambda _{AE} (\times 100)$
, and
$NoC$
is the number of realizations with successful convergence in iteration for the maximum likelihood estimation
Summary of estimates (median (q1,q3)) and empirical coverage probabilities (CP) of the exponential parameter and frequencies that each model was selected with AIC for the time-to-AE for Dataset 2A (S = 13, dependence): TRUE implies the value of
$\lambda _{AE} (\times 100)$
, dep. shows the applied copula, and
$NoC$
is the number of realizations with successful convergence in iteration for the maximum likelihood estimation
Table 6 presents the summary of the estimates and the coverage probabilities for the hazard parameter with the proposed methods under the independent and dependent follow-up duration assumptions for Dataset 1A. When assuming independence correctly for Dataset 1A, the proposed method had successful convergence for all the simulated dataset (it is given as
$NoC$
in Table 6) in maximizing the likelihood function. We found that the bias on
$\lambda _{AE}$
was negligible and the empirical CP was close to the nominal level in all scenarios. Although modeling dependence with copula functions is unnecessary for independent data (Dataset 1A), predetermining independence is challenging. Thus, we evaluated performance of the method under the dependent follow-up duration applied to the independent data. We observed that the proposed method sometimes could not obtain estimates in particular for scenarios 1 and 2 and for Frank copula. We also observed that the extra modeling of dependence with the Clayton copula led negligible biases for
$\lambda _{AE}$
and provided conservative, but still valid confidence intervals. With the Frank and Gumbel copulas, there were certain biased when the event rate was low (scenario 1) and the empirical coverage probabilities were far from the nominal level except when the event rate was high (scenario 3). Table 6 also shows frequencies that each model was selected with AIC. In approximately 75
$\%$
simulated datasets, the independence model was correctly selected across all the scenarios.
The results for Dataset 2A were summarized in Table 7. When the event rate was moderate (scenario 2) or high (scenario 3), the proposed method had only negligible biases on
$\lambda _{AE}$
with a correctly-specified copula function. The model with correctly-specified copula provided the empirical coverage probabilities close to the nominal level except for some cases such as Frank with
$\psi _k$
=1.87 in scenario 2. On the other hand, when the event rate was low (scenario 1), even if the supposed copula model covered the true dependence, there might be certain biases and empirical coverage probabilities away from the nominal level. Table 7 also indicates that misspecified copula modeling of dependence might lead seriously biased estimates for
$\lambda _{AE}$
without relevant control of the nominal level for the confidence interval. Frequencies that each model was selected with AIC were presented in Table 7. When the event rate was low (scenario 1) or dependence was low with Knedall’s
$\tau $
=0.2, the independence model was likely to be selected although the datasets were generated from dependent copulas. If event rate was high (scenario 3), AIC selected the true copula the most frequently.
Summary of estimates (median (q1,q3)) and empirical coverage probabilities (CP) of the exponential parameter for the time-to-AE and frequencies that each model was selected with AIC for Dataset 1B (S = 30, independence):
$True$
implies the value of
$\lambda _{AE} (\times 100)$
, and
$NoC$
is the number of realizations with successful convergence in iteration for the maximum likelihood estimation
Summary of estimates (median (q1,q3)) and empirical coverage probabilities (CP) of the exponential parameter and frequencies that each model was selected with AIC for the time-to-AE for Dataset 2B (S = 30, dependence): TRUE implies the value of
$\lambda _{AE} (\times 100)$
, dep. shows the applied copula, and
$NoC$
is the number of realizations with successful convergence in iteration for the maximum likelihood estimation
The results under the independent follow-up duration with S = 30 (Dataset 1B) are presented in Table 8. If the independent copula is correctly applied, the proposed method had negligible biases and empirical coverage probabilities close to the nominal level similarly to the case of S = 13 (Table 6). When dependent copulas were applied, increasing the number of studies led substantial reduction of biases. The empirical coverage probabilities might be anti-conservative for a rare event case (scenario 1) even with S = 30; when dependence was supposed with Frank and Gumble copulas to the dataset under independence, empirical coverage probabilities were far from the nominal level and anti-conservative. For scenarios 2 and 3, our methods had empirical coverage probabilities close to the nominal level. The results under the dependent follow-up duration with S = 30 (Dataset 2B) are presented in Table 9. The proposed method had negligible biases when the copula was correctly specified. With a larger number of studies, interquartile ranges of estimates for
$\lambda _E$
became shorter. When the event rate was low (scenario 1), empirical coverage probabilities might be far from the nominal level. Overall, the proposed method had empirical coverage probabilities close to the nominal level with the correctly specified copula.
In Table 10, we gave a summary of estimates of the exponential hazard parameter for the censoring distribution of each study. Recall that the at-risk data were available at 6, 12, and 18 months for all the studies in Case 1 and at randomly selected two time points in Case 2. In both cases, all the exponential parameters were appropriately estimated. It suggests that estimation of the censoring distribution with the method in Section 3.2 is not sensitive to what time points are utilized.
6 Discussion
Since very limited safety assessment can be made in a single clinical trial, meta-analysis would be a very powerful tool to clarify the safety profile of drugs. An FDA Guidance for Industry 7 argues various issues in meta-analysis of AEs in randomized clinical trials. The issue of heterogeneous follow-up duration is repeatedly emphasized (see Section III-D in the guidance). Importance of careful consideration of eligibility criterion in study selection is emphasized; without subject-level data, some trials may need to be excluded in the presence of heterogeneous follow-up duration. However, no specific statistical method is mentioned in the guidance to handle the issue except for recommending use of the IPD if possible. To our best knowledge, as of this moment, no established methodologies exist for literature-based meta-analysis of AEs that accounts for heterogeneous follow-up and our proposed method is the first one to handle this situation.
The idea of our proposed method is to utilize the Kaplan–Meier estimates of the primary time-to-event efficacy endpoints as information of the follow-up duration. Displaying the Kaplan–Meier estimates of the OS or the PFS is routinely used in oncology clinical trials. Then, our method is of particularly practical benefit in oncology clinical trials, in which the heterogeneous follow-up duration often matters in the nature of the study design. In oncology clinical trials, the Kaplan–Meier estimates for both of the OS and the PFS may be reported. Note that our idea is to obtain the distribution of the end of follow-up of AEs with them, then we should choose the endpoint (the OS or the PFS) defining the end of the follow-up in the study protocol. The follow-up duration for the AE is determined by other censoring, such as the treatment discontinuation due to some safety issues or lack of efficacy. In most oncology clinical trials, the subject is off the study once the follow-up of the efficacy endpoint is stopped. Then, such subjects are regarded as censored in evaluating the primary time-to-event endpoint. Since the distribution of the follow-up duration was determined by the distribution of the efficacy endpoint and that of the censoring as seen in Section 3.1, the treatment discontinuation is appropriately handled in the proposed method. Although we focus on oncology in this manuscript, there are some other therapeutic areas of potential application of the proposed method, possessing the features discussed in this paragraph, such as cardiovascular disease.
Summary of the estimates of the study-specific exponential parameter for the censoring distribution
Through numerical studies, we found that our method works well if parametric models are appropriately employed. If misspecified, the estimate may not be obtained due to failure in iteration steps for the maximum likelihood estimation and even if the estimate is successfully obtained, it may be biased. When the independence is satisfied, the extra modeling of the dependence might lead poor performance when the number of studies was small (S = 13). In our simulation study, we observed that the Clayton copula does not lead any invalid inference, whereas with the Frank and Gumbel copulas, the proposed methods might suffer from bias in estimation, poor control of the nominal level of the confidence interval and non-convergence in iteration steps for estimation probably due to difficulty to handle extra parameters in dependence with limited information from summary statistic data. When the number of studies was large (S = 30), the performance of the model improved; the extra copula modeling of dependence did not lead serious deficiencies.
Here, we summarize the steps and the points to be taken care of in applying our proposed method. Recall our proposed method for meta-analysis occurs within treatment groups. To make a valid comparison, the subjects’ characteristics should be similar between the groups and then, careful consideration is needed on inclusion of studies. To apply the proposed method, we need frequencies of the AE of each study and the survival functions of the end of the follow-up for the AE. The latter can be obtained from a graphical plot of the Kaplan–Meier estimate for the time-to-event efficacy primary endpoint (OS or PFS). To estimate the censoring distribution for the efficacy endpoint, information on the number of at-risk, which is often attached below the graphical plot of the Kaplan–Meier estimates, should be obtained (see Table 2). The key quantity to summarize the occurrence of the AE is the cumulative distribution function of the time to occurrence of the AE. Key steps to apply the method are as follows:
1. Specify parametric models for the distribution of the time to occurrence of the AE,
$F_{AE,k}(t; \theta _k)$
.
2. Specify a parametric model for the censoring distribution
$G_k^{(s)}(t; \eta _k)$
and estimate the unknown parameter
$\eta _k$
with the method in Section 3.2.
3. Specify copula families
$C_{k, \psi _k}(u,v)$
for the dependence between the time to occurrence of the AE and the duration of the follow-up for the AE.
4. Estimate the unknown parameter
$\theta _k$
with the method in Section 3.1 under the assumption of the independent follow-up duration. Estimate
$\theta _k$
, as well as the unknown parameters of the copula, with the method in Section 3.3 under the dependent follow-up duration.
5. Select the best combination of the models for the cumulative distribution function for the time to occurrence of the AE and that for the copula attaining the minimum AIC value. The goodness-of-fit can be evaluated by the graphical plots of the observed and model-based incidence proportions with the chi-squared goodness-of-fit test as illustrated in Section 4.
6. Present the cumulative distribution functions
$F_{AE:k}(t; \hat {\theta }_k)$
and their RR at several time points as done in Figures 2 and 4 and Table 4.
One may apply the Poisson regression model with the follow-up durations of each study as the offset term if the follow-up duration is available for all the studies. This was suggested by Cai et al.,Reference Cai, Parast and Ryan 4 although they did not discuss how to determine the follow-up duration in practice. It can be obtained with Kaplan–Meier estimate of the efficacy endpoint and the censoring distribution with the method discussed in this article. The Poisson regression assumes the exponential distribution for the time-to-AE and then our method can handle more flexible class of the distributions for the AE.
The motivation of this research came from difficulty in interpretation of the incidence proportion in the presence of heterogeneous follow-up duration. However, we do not criticize importance of the incidence proportion in practice. It is very useful to capture overall tendency in safety profiles of drugs. On the other hand, when one need to make delicate evaluation of the usefulness of drugs considering balance between efficacy and risk of AEs, it is not appealing to simply compare incidence proportions among the treatments. In oncology clinical trials, the follow-up duration is often restricted by the time-to-event efficacy endpoint and then superiority with respect to the primary endpoint leads longer follow-up duration of AEs. Thus, the drug of superior efficacy may have overly estimated incidence proportions of AEs. To compare the risk and benefit of drugs, proper adjustment of the heterogeneous follow-up durations should be considered and the proposed method would be useful in such situations.
In our empirical data analysis, the precision of fit was compromised due to the inclusion of trials with no AEs in the control group. In the simulation study, we also observed that our proposed method might not work well for AEs of a low event rate. However, our method worked well for AEs of a certainly high event rate. In addressing the risk and benefit of drugs, we usually focus on commonly observed AEs. In such situations, our method worked well and would be useful in practice.
On the other hand, one may need to compare rare AEs among several drugs and the issue of heterogeneous follow-up durations can also have serious impacts on the comparison. Then, it is important to extend our inference procedure so as to provide stable estimates for rare events. To achieve this, it will be necessary to explore extensions, such as Bayesian estimation or penalized inference methods, which we consider as future research topics.
A limitation of our method is it assumes the common cumulative incidence function for the AE over studies. It corresponds to the fixed-effect model in meta-analysis. On the other hand, one may hope to evaluate heterogeneity of occurrence of AEs over studies and then the mixed-effect models would be more appealing. It would warrant to develop inference procedures to incorporate random effects to evaluate heterogeneity. In the proposed method, we fit the model and then estimate the cumulative incidence over time separately within the groups. It allows to fit difference parametric models to each group flexibly. On the other hand, it also implies that outcomes from the same randomized study may be analyzed compromising randomization. This issue is argued in the context of network meta-analysis. In the standard pairwise meta-analysis of randomized trials, a treatment contrast, such as the log-transformed odds ratio between the two treatment groups, is used as the outcome for meta-analysis respecting randomization. This approach is called the contrast-based approach since the outcome is the treatment contrast estimated within a study. In network meta-analysis, studies of various combinations of treatments are included and then background characteristics may not be homogeneous over the treatment groups to compare even if the contrast-based approach is taken. Alternatively, the arm-based approach is of interest and its pros and cons were discussed; White et al.Reference White, Turner, Karahalios and Salanti
19
reported that the arm-based approach may lead to more biased results than the contrast-based one. Our method corresponds to the arm-based method and then studies included in meta-analysis should be carefully determined. An approach is to include only randomized studies with both treatment groups to be compared in meta-analysis or studies with only either group are excluded, respecting randomization. On the other hand, it is desirable to include a wider range of studies in meta-analysis. To this end, it is important to develop methods for adjusting confounders, possibly through regression modeling, for example, for
$\theta _k$
in
$F_{AE: k}(t; \theta _k)$
, that incorporate Bayesian estimation and penalized inference techniques to ensure stable implementation.
Same as for efficacy evaluation, the publication bias may be a serious threat in meta-analysis of AEs.Reference Kicinski, Springate and Kontopantelis 20 – Reference Zorzela, Loke and Ioannidis 22 For the proposed method, we need to consider other type of reporting bias. Our method requires the Kaplan–Meier estimates for the efficacy time-to-event endpoints and they may not be reported; as explained in our example, 5 studies among 17 studies did not report the Kaplan–Meier estimates and then they were excluded in our meta-analysis. It is very valuable to develop methods to address potential biases due to such kind of selective sampling in our proposed method.
Acknowledgements
The authors are grateful to the editors and reviewers for their careful reading of the manuscript and for insightful suggestions and comments.
Author contributions
Conceptualization: S.H.; Methodology: S.K., S.H.; Software: S.K.; Writing—original draft: S.K.; Writing—review and editing: S.H.
Competing interest statement
The authors declare that no competing interests exist.
Data availability statement
The dataset was created using a graph scan software with a dataset available to the public. It is available at 10.5281/zenodo.16872253 with an SAS code to implement the proposed method.
Funding statement
The last author’s research was partly supported by Grant-in-Aid for Scientific Research (25K03086) from the Ministry of Education, Science, Sports and Technology of Japan.
Open access
