MAIN PAPER (wileyonlinelibrary.com) DOI: 10.1002/pst.1624
Published online 16 June 2014 in Wiley Online Library
Missing data sensitivity analysis for recurrent event data using controlled imputation Oliver N. Keene,a * James H. Roger,b Benjamin F. Hartley,a and Michael G. Kenwardb Statistical analyses of recurrent event data have typically been based on the missing at random assumption. One implication of this is that, if data are collected only when patients are on their randomized treatment, the resulting de jure estimator of treatment effect corresponds to the situation in which the patients adhere to this regime throughout the study. For confirmatory analysis of clinical trials, sensitivity analyses are required to investigate alternative de facto estimands that depart from this assumption. Recent publications have described the use of multiple imputation methods based on pattern mixture models for continuous outcomes, where imputation for the missing data for one treatment arm (e.g. the active arm) is based on the statistical behaviour of outcomes in another arm (e.g. the placebo arm). This has been referred to as controlled imputation or reference-based imputation. In this paper, we use the negative multinomial distribution to apply this approach to analyses of recurrent events and other similar outcomes. The methods are illustrated by a trial in severe asthma where the primary endpoint was rate of exacerbations and the primary analysis was based on the negative binomial model. Copyright © 2014 John Wiley & Sons, Ltd. Keywords: missing; sensitivity; recurrent event; exacerbation; multiple imputation; MNAR
1. INTRODUCTION
258
The potential for bias in study conclusions introduced by missing data has become an important aspect in the interpretation of clinical trials designed to demonstrate the effectiveness of new medicines. The National Research Council of the National Academy of Science report on the prevention and treatment of missing data in clinical trials [1] discussed many of these issues and in particular recommended that ‘examining sensitivity to the assumptions about the missing data mechanism should be a mandatory component of reporting’. In a follow-up paper in the New England Journal of Medicine [2], the authors stated that analysts of clinical trials need to ‘decide on a primary set of assumptions about the missing-data mechanism...In some cases, the primary assumptions can be that data are missing at random’ and to ‘assess the robustness of inferences about treatment effects to various missing-data assumptions by conducting a sensitivity analysis that relates inferences to one or more parameters that capture departures from the primary missing data assumption’. The ‘missing at random’ (MAR) assumption requires that conditional on the data observed for each patient including their covariates in the model (together making their ‘history’), their unavailable data are randomly missing. Or equivalently, for withdrawals, the conditional future statistical behaviour of outcomes given the history is the same for those who remain in the trial and those who withdraw, provided their treatment is the same throughout and they share the same history. This implies that, under MAR, future statistical behaviour for those who withdraw can be modelled from those who remain. It follows that, under MAR, if data are collected only when patients are on their randomised treatment, the resulting estimator of treatment effect
Pharmaceut. Statist. 2014, 13 258–264
corresponds to the situation in which the patients adhere to this regime throughout the study. This has been referred to as the de jure estimand of treatment effect [3] and elsewhere as an efficacy estimand [4], although others view the ‘efficacy/effectiveness terminology’ as confusing [3,5]. A MAR analysis is often labelled as an intention-to-treat (ITT) analysis, but it is actually inconsistent with a strict interpretation of the ITT principle [6]. A full ITT analysis is one that includes all data for a patient regardless of what treatment they take and, if the patient withdraws, must reflect the actual treatment taken following withdrawal [7,8]. Such an analysis including data collected after cessation of randomised treatment would provide a de facto estimand of treatment effect [3] and has been labelled as an effectiveness estimand because it would compare the outcomes of treatment strategies rather than the actual pharmacological effects of treatment itself. An alternative approach to derive such a de facto estimand compares the outcome of treatment strategies by making assumptions about the possible outcome when subjects switch to control, without using the actual off-treatment data that may include the confounding influence of rescue medication received after withdrawal [3,4]. The issue of missing data is a focus of attention for regulatory review of confirmatory trials, and the requirement for sensitivity analyses is clearly stated in the Committee for Medicinal Proda
GlaxoSmithKline Research and Development, Middlesex, UK
b
Medical Statistics Department, London School of Hygiene & Tropical Medicine, London, UK *Correspondence to: Oliver N. Keene, GlaxoSmithKline Research and Development, 1-3 Iron Bridge Road, Uxbridge, Middlesex UB11 1BT, UK. E-mail: [email protected]
Copyright © 2014 John Wiley & Sons, Ltd.
O. N. Keene et al. ucts for Human Use guideline [9], which states, ‘as the choice of primary analysis will be based on assumptions that cannot be verified it will almost always be necessary to investigate the robustness of trial results through appropriate sensitivity analyses that make different assumptions’. Most developments in this area of sensitivity analyses have been for longitudinal continuous outcomes. Little and Yau [10]; Carpenter, Roger and Kenward [3]; Mallinckrodt [11] and O’Kelly and Ratitch [12] describe the use of multiple imputation (MI) methods based on pattern mixture models to investigate the impact of different assumptions about the statistical behaviour of postwithdrawal outcomes. For trials that compare a test and reference product, this approach allows modelling for the postwithdrawal outcome for the test arm to reflect the behaviour of later outcomes on the reference arm. Similar approaches where the postwithdrawal outcomes on the active arm has been chosen to reflect the corresponding observed outcomes on the placebo arm have been called ‘placebo MI’ [4]. The term ‘controlled imputation’ has been used to describe all of these approaches [5]. Recurrent event data occur in clinical trials when multiple events are observed for each patient. Examples of recurrent events include exacerbations in asthma and chronic obstructive pulmonary disease and bone fractures in osteoporosis. Analysis that takes account of the repeated nature of these events is to be preferred to a simple time-to-first-event analysis [13]. A valuable approach that allows simultaneous modelling of the event rate and the variability of this rate between patients is the use of the negative binomial model with log link [14]. This paper describes an extension of these controlled imputation methods for the analysis of recurrent event data. The methods are illustrated by an example in severe asthma.
2. EXAMPLE TRIAL The DREAM trial was a study of mepolizumab in severe asthma [15]. The primary endpoint was the frequency of clinically significant asthma exacerbations, collected over a 1-year treatment
period. The trial randomised subjects equally to four treatment groups: placebo and three doses of mepolizumab. Of the 616 subjects recruited to the trial, 520 (84%) completed treatment. Withdrawal rates were 18% in the placebo group, 16% for mepolizumab 75 mg, 14% for mepolizumab 250 mg and 15% for mepolizumab 750 mg. The time to withdrawal is shown in Figure 1. Reasons for withdrawal of the 96 subjects who withdrew were adverse event (28 subjects), withdrawal of consent (28 subjects), lack of efficacy (22 subjects), lost to follow-up (6 subjects) and other reasons (12 subjects). There was no obvious imbalance in numbers of patients reporting each reason across treatment groups. Each dose of mepolizumab gave a similar reduction in the frequency of exacerbations and a similar level of withdrawal so for simplicity these dose groups have been combined to illustrate the impact of the different methods of analysis. The primary analysis was performed using a negative binomial generalised linear model with a log link function. For this model, as for Poisson regression, the response variable is the total number of exacerbations and time at risk of an exacerbation is fitted as an offset variable. Because the response variable is the total number of events, this model and Poisson regression essentially assume that events occur at a constant rate within an individual and that there is no serial correlation. In the negative binomial model, variation in event rate across patients is accounted for through use of an explicit dispersion parameter in the model. In contrast, for Poisson regression, the use of an overdispersion correction is an alternative way of inflating the variance to allow for this variation from patient to patient. The merits of using the negative binomial to analyse exacerbation rates have been discussed elsewhere [14]. The covariates used in the model were treatment group and other important baseline predictors of outcome. The model included covariates for treatment group, use of maintenance oral corticosteroids, region, number of exacerbations in the previous year and baseline percentage of predicted FEV1 . The log of the
Pharmaceut. Statist. 2014, 13 258–264
259
Figure 1. Time to withdrawal in the DREAM trial.
Copyright © 2014 John Wiley & Sons, Ltd.
O. N. Keene et al. length of time observed in the trial for exacerbations was used as an offset variable. This likelihood-based analysis is valid under the MAR assumption, assuming that the model is correct and estimates the treatment effect that corresponds to what would be observed if those who withdrew had continued to take their randomised treatment throughout the study. That is, it forms a de jure analysis. It is of interest to investigate the sensitivity of the conclusions of this trial to different assumptions regarding the implied behaviour of those that withdrew, that is, to specific departures from the MAR assumption (de facto analysis).
3. TECHNICAL DETAILS OF SENSITIVITY ANALYSIS The negative binomial distribution is conventionally defined as the probability distribution of the number of successes y before k failures are seen in a series of independent Bernoulli trials with probability p of success and .1 p/ of failure. Equation (1) is the resulting probability of y successes.
PŒY D y D
.y C k/ .1 p/k p y .y C 1/.k/
(1)
The distribution has mean D kp=.1 p/ and variance kp=.1 p/2 . The inverse of the parameter k is often referred to as the dispersion parameter as the variance can be reexpressed as . C 2 =k/, noting that the variance for the Poisson distribution is simply . Confusingly, the inverse of the parameter k is sometimes also denoted by k (e.g. in documentation for the SAS PROC GENMOD procedure). For recurrent events, the negative binomial distribution can be derived as a mixture of Poisson distributions where the mixing distribution acts multiplicatively on the intensity œ of the Poisson process to give an effective intensity œ” conditional upon ”, where ” has a gamma distribution with mean k. The negative binomial is the marginal distribution of the number of occurrences after integrating out gamma. The link between œ and p is that p D œ=.1 C œ/ and in reverse œ D p=.1 p/. This recurrent event definition allows for noninteger values of k. In this context, it is often called a gamma-Poisson distribution. The distribution can be extended to cover a series of distinct outcomes. The negative multinomial distribution is the probability distribution of the number of occurrences .y1 , y2 , : : : , ym / of m types of success before we see y0 failures in series of independent trials where probability of outcome j is pj for j D 0, 1, : : : , m and P failure is outcome zero with probability p0 D 1 jD1,:::,m pj . Here, y0 takes the role of the parameter k. The joint probability is given in equation (2). y0 PŒY D . y1 , y2 , : : : , ym / D Qm
P
m jD0 yj
jD0 . yj C 1/
Ym jD0
yj
pj
(2)
260
As with the negative binomial distribution, the negative multinomial distribution can be derived as a series of m independent Poisson values Yj sampled with rates œj ” using the same value
Copyright © 2014 John Wiley & Sons, Ltd.
” from the mixing distribution, with œj D p0 =pj . Then, p0 D 1 P 1 C j 1=œj and pj D p0 =œj . The negative binomial distribution is a special case where m D 1. In this paper, we only use the extension of negative binomial to the negative trinomial as we only consider two periods, before and after withdrawal. The underlying negative multinomial theory can be used to extend this model to a series of predefined discrete periods where the event rate within an individual is assumed constant within each period. For each subject, the period in which withdrawal occurs would then be split into two as we do here. For recurrent events, the mixing parameter ” is a property of the individual patient, and the Yij are the number of events in a series of distinct periods indexed by j for the ith patient. The assumption is that, conditional on the patient’s underlying propensity ”, the numbers of events seen in each period are independent. In particular, there is no serial correlation. The impact of a patient’s propensity works proportionally on the rate predicted by the other components of the model, such as treatment, baseline covariates and length of time at risk. Under this model, the underlying rate of events within a trial might change from period to period and also be modified by covariates. For a patient i with covariates Xi the numbers of events Y1 , : : : , Ym in periods j D 1, : : : , m have independent Poisson distribution with rate: i exp.j C Xi ˇ/. It follows from the gamma-Poisson derivation that the marginal distribution of any subsetP S of the Yj is also nega tive multinomial, with pj D pj = s2S ps C p0 and p0 D P p0 = s2S ps C p0 . So each individual Yi has a marginal distribution which is negative binomial with parameter p D pj =.pj C p0 / and the same k. Similarly, any weighted sum of the Yj also has a negative binomial distribution with parameter p D P P P j2S wj pj = j2S wj pj C p0 j2S wj where k remains the same. The gamma-Poisson model is a special case of a wider class of model. Consider the classic generalised linear mixed model with Poisson distribution and log link. In this model, repeated observations Yit are made on the ith subject in interval t: Yit jsi Poisson.it / and ln.it / D Xit ˇ C si where si is a subject effect with some distribution on the real line. The design matrix Xit describes both the baseline covariates, as well as the imposed treatment structure. Set ui D exp.si /, and it D exp.Xit ˇ/, then: Yit jui Poisson.ui
it /
The margin distribution PŒYit D Ry for Yit is obtained by integrating out the random effect ui , as PŒYit D yjui f .ui /dui The negative binomial is a special case where the distribution of ui is Gamma, while a very similar distribution is generated when si is normally distributed and ui has a log-normal distribution. Clearly, any weighted sum across intervals t for subject i has a P similar distribution within the family with mean t .ui it /. The distribution of Yij conditional upon Yik for the same subject i is
Pharmaceut. Statist. 2014, 13 258–264
O. N. Keene et al. presented in equation (3) where §1 D §ik and §2 D §ij are the event rates for Yik and Yij .
PŒYij D y2 jYik D y1 D
.y1 Cy2 /Š y1 Šy2 Š. 1 C
y1 1 2
y2 2 y / 1Cy2
PŒYij CYik D y2 Cy1 PŒYik D y1
(3) This conditional probability is obtained by noticing that the joint distribution of Yij and Yik can be factored into PŒYij D y2 , Yik D y1 D PŒYij D y2 , Yik D y1 jYij C Yik D y2 C y1 PŒYij C Yik D y2 C y1 (4) where the first term on the right-hand side of (4) describes the allocation of .y1Cy2/ events between the two intervals as a binomial distribution with probabilities §2 =.§1 C §2 / and §1 =.§1 C §2 /. This result relies on the nature of the log link function and is irrespective of the mixing distribution The importance of this result for missing data is that the expression (3) is free of any requirement to integrate out the actual random effect u. That is, the binomial term relies only on the relative sizes of §1 and §2 , while the remaining part involves simple marginal distributions for Yik and .Yij CYik /. This contrasts strongly with the analogous multivariate normal setting, where such a random patient effect is not eliminated through conditioning on the observed part of the sequence. Here, the impact of the distribution of the random effect u on the conditional distribution is summarised through the ratio of the probability distribution of the sum and of the observed first observation. That is, we only require the marginal impact of the distribution of the random effect u on each of these to obtain the marginal impact on the conditional distribution. This means that MAR and MAR-related methods can be implemented simply through these marginal distributions, ignoring the underlying random effects that have been conditioned out. The importance is however limited by the practicality of sampling from the conditional distribution (3) whenever an explicit expression is not available for PŒYij C Yik D y2 C y1 as a function of y2 given y1 . For this reason, we limit consideration to the following special case and ignore the otherwise interesting case where ui has a log-normal distribution. This later might be more easily solved by a data augmentation approach using an MCMC approach. When the mixing distribution is gamma, then both Yij and Yik have negative binomial distributions with probabilities given by equation (1) with the probability p taking values §1 =.1 C §1 / and §2 =.1 C §2 / giving a general form for probability of y events as expressed in terms of § in equation (5): PŒY D y D
.k C y/ y .y C 1/.k/.1 C /kCy
(5)
The conditional distribution of Yij given as Yik D y1 given in equation (6) that can be derived directly from equations (3) and (5) as both Yik and .Yij C Yik / also have negative binomial distributions with § taking values §1 and .§1 C §2 /.
PŒYij D y2 jYik D y1 D
Pharmaceut. Statist. 2014, 13 258–264
(6)
PŒYij D y2 jYik D y1 D
..kCy1 /Cy2 /. .y2 C1/.kCy1 /.1C
y 2 =.1C 1 //2 .kCy1 /Cy2 2 =.1C 1 //
(7) The equivalent result for the negative multinomial has y1 and §1 replaced by the sum of their equivalents in the conditional statement. This property is extremely useful for handling missing data, because the conditional distribution of the missing values for partially observed negative multinomial data is also negative multinomial. Missing data can be imputed using the appropriate negative multinomial models. In its simplest form for any specific subject, conditional upon their own withdrawal time t, the number of events before and after withdrawal can be modelled under a MAR assumption as having a joint negative trinomial distribution, with observed Y1 before withdrawal and unobserved Y2 after withdrawal This is equivalent to the classic primary analysis where Y is assumed to have a negative binomial distribution before withdrawal. We simply extend the model to include the unobserved number of occurrences between withdrawal and study end. Ordinarily, an analysis under a MAR assumption is implemented using direct likelihood, by ignoring the unobserved Y2 and introducing an offset into the log-linear model to allow for the length of the observation period. Using the aforementioned results, an alternative MI approach can be used in which values for Y2 are imputed using the distribution of Y2 conditional on Y1 D y1 . This requires fitting a Bayesian log-linear negative binomial model with offset to the observed data and drawing independent samples from the posterior distribution for parameters k and “ created by multiplying a noninformative prior with the likelihood function. It is assumed for the subsequent MI steps that the priors are sufficiently vague and the sample size large enough, for the data to dominate the posterior. This stage can be implemented simply in SAS using the BAYES statement in the GENMOD procedure. Alternatively, MCMC can be used programming the likelihood directly in either Winbugs or its derivatives or in the MCMC procedure in SAS. For each imputed data set, one sampled set of values k and “ is used to generate a single sample of values y2 for each missing value Y2 conditional upon their associated observed Y1 D y1 . This stage can be programmed using a function to sample from the negative binomial distribution such as the SAS RAND function. The imputed and observed values for each patient are summed to form Y D y1 C y2, and these are analyzed as though the data were complete. For instance, least-square mean treatment differences can be extracted using the GENMOD procedure for each imputation using the BY statement. The resulting estimates and their standard errors are combined across imputations using Rubin’s rules [16]. This last stage can be performed using the MIANALYZE procedure in SAS. This procedure will be valid under MAR and gives results that are very close to those obtained by direct likelihood It is important to note that this approach is not introduced simply as an alternative to using the standard direct likelihood
Copyright © 2014 John Wiley & Sons, Ltd.
261
y2 kCy1 2 .kCy1 Cy2 /.1C 1 / .y2 C1/.kCy1 /.1C§1 C§2 /kCy1Cy2
This conditional distribution is also negative binomial with the k parameter k D k C y1 , while § D §2 =.1 C §1 /, and the mean is §2 .k C y1 /=.1 C §1 / This can be seen by rearranging equation (6) in the form (7) and identifying its relationship to equation (5).
O. N. Keene et al. approach to a MAR analysis. Rather, it provides a convenient framework for sensitivity analysis that can explore the implications of particular MNAR departures from MAR. Sensitivity analyses can be based on the pattern-mixture approach, through modifications of the postwithdrawal imputation models, either of the following: (i) by using a reference-based approach and borrowing parameters from the reference arm or (ii) by using a ‘Delta method’ approach and explicitly altering the rate from the prewithdrawal period [17] In the reference-based approach (i), the rate in the reference arm is used to impute the missing data for both arms. The model for the postwithdrawal outcomes is modified from what would be used under MAR but, most importantly, remains calculable from some fixed linear function of the parameters “ in the original model. Approach (ii) in the current setting would be implemented by multiplying the event rate by a specific ratio, the delta, after withdrawal, equivalent to adding a fixed constant to the linear predictor in the model for §2 . In the sensitivity analysis, the total number of events (observed prewithdrawal + imputed postwithdrawal) is analysed over the full treatment period using the original primary analysis model. For the reference-based approach, possible options include the jump to reference (J2R) and copy reference (CR) approaches analogous to those used for longitudinal normal data (3). The J2R approach represents the situation where the patient’s expected event rate is shifted to that of the reference arm. Such a change may be seen as extreme and might be used as a worst-case scenario in terms of reducing any treatment effect. In this approach, it is assumed that withdrawn patients on the active treatment will quickly lose the positive effect of their treatment when it is stopped. The impact of sampling from the conditional distribution is that if their event rate prior to withdrawal is worse than would be expected (positive residual) on the test treatment; their imputed event rate after withdrawal will be worse than the expected event rate on the reference treatment. Postwithdrawal events in the reference arm are imputed under randomised arm just as under MAR.
In the CR approach, a patient’s expected event rate both before and after withdrawal is assumed to be the same as the reference group. This mimics the case where those withdrawing are in effect nonresponders. This has a less extreme impact than the aforementioned J2R. When a patient on test treatment has more events before withdrawal than expected for them if they had been in the reference arm, then this ‘positive residual’ will feed through into a higher than expected event rate in the postwithdrawal period. This is because their earlier observed higher than normal event rate suggests a personal high propensity to events. Postwithdrawal data in the reference arm are again imputed under randomised arm MAR. An alternative approach based simply on results from the reference arm is unconditional reference (UR). The basis of this approach is that withdrawal from test treatment represents a new episode for the patient, and the previous history of events is not used in the imputation model for events postwithdrawal. Instead, events postwithdrawal for the test arm are imputed using the overall mean for the reference arm, conditional only on baseline covariates. Postwithdrawal data in the reference arm are again imputed under randomised arm MAR. Alternatively, the reference arm could be imputed unconditionally, on the basis of an argument of sharing imputation approaches across arms. This approach may work well when there are other informative covariates other than the previous event history, on which to condition, and make the implicit MAR assumption likely to hold. One of the great advantages of these approaches to sensitivity analysis is that different scenarios can be used for different patients, depending on the reasons for withdrawal. For example, it may be reasonable to use a MAR assumption for patients who withdraw because of moving location and a J2R approach for patients who withdraw for a reason potentially related to treatment such as an adverse event or lack of efficacy. Whichever imputation model is chosen, the subsequent procedure follows that of standard MI. For each imputation, a single set of data is created for each subject’s missing data conditional upon their observed data. After summing observed and imputed
262
Figure 2. Event rate ratios from sensitivity analyses.
Copyright © 2014 John Wiley & Sons, Ltd.
Pharmaceut. Statist. 2014, 13 258–264
O. N. Keene et al. values for each subject, each imputed data set is analysed using the model used for the primary analysis. At least 1000 imputations are recommended and more should be considered, to provide sufficiently stable and reproducible results. We note that when deviations from MAR are included in the imputation model, then the imputation model is what is termed ‘uncongenial’ [18]. Details of how to implement the approach in the SAS system are given in the appendix (available online as Supporting Information).
4. RESULTS OF SENSITIVITY ANALYSIS These sensitivity analyses have been implemented for the example DREAM trial described earlier. Figure 2 shows a Forest plot of the results. The MAR imputation model provides essentially the same results as that for the direct likelihood model. As expected CR, J2R and UR approaches reduce the size of the treatment difference, but the trial remains overwhelmingly positive in terms of size of reduction in exacerbations. For the direct likelihood analysis, the ratio of mepolizumab to placebo is 0.53 (95% CI: 0.43–0.67), while the UR approach gives an estimated ratio of 0.58 (95% CI: 0.46–0.73). The figure also includes an analysis based on reason for withdrawal (labelled on the plot as MAR and J2R imputation). For this analysis, imputation was performed as MAR for patients with a reason viewed as plausibly unrelated to treatment (e.g. lost to follow-up) and J2R for reasons potentially more plausibly related to treatment (e.g. adverse event). This approach provides a result intermediate between the MAR and J2R approaches. The plot also shows the result of one delta method imputation. In this case, the predicted postwithdrawal rate was doubled for active patients compared with active patients remaining in the study, while no delta was applied to placebo patients. This approach provides a result similar to the UR and J2R approaches.
5. DISCUSSION
Pharmaceut. Statist. 2014, 13 258–264
Copyright © 2014 John Wiley & Sons, Ltd.
263
The use of controlled imputation methods for the analysis of recurrent event data allows sensitivity analyses to be performed, which investigate the consequences of departures from the MAR assumption. These departures correspond to different assumptions about postwithdrawal treatments and so map to different scenarios underlying de facto estimators. Alternative approaches to sensitivity analyses for nonnormally distributed data have often relied on filling in the missing data using an ad hoc rule. For example, the Committee for Medicinal Products for Human Use guideline [9] states that for recurrent events, ‘baseline count imputation could also be a plausible approach’. Such an approach is analogous to the baseline observation carried forward (BOCF) approach for longitudinal normally distributed data. These methods are commonly described as ‘single imputation’ and are problematic for several reasons. They are valid only under specific and arguably unrealistic, assumptions [19], and they do not appropriately reflect the variability inherent in the process of imputation. They can mask treatment effects when these exist and create treatment effects when they do not [20]. Molnar, Hutton and Ferguson [21] provide a critique of the use of last observation carried forward (LOCF)-type methods in the neurosciences. Mallinckrodt et al. [22] have stated: ‘Regarding the primary analysis for confirmatory longitudinal clinical trials, conclusive evidence has demonstrated the need to abandon the simple, ad hoc methods such as LOCF and BOCF’. For recurrent
event data, the MI approaches described here are similarly clearly preferable to ad hoc methods. Basing sensitivity analysis on data collected on the reference group (e.g. placebo) has clear attractions as this can reasonably be thought of as what might happen to a patient who discontinues treatment on the test product. However, options remain in how to implement imputations based on the reference group, depending on the nature of the treatment and condition under consideration. For example, does the treatment cause permanent changes or does the treatment have only short-term benefits, as with an analgesic? This paper describes three possibilities: J2R, CR and UR. The J2R imputation method can be thought of as a ‘worst case’ approach; patients who discontinue active medication and who are worse than expected on the test arm are then imputed as worse than expected on the reference arm. More reasonable may be the CR or UR approaches where imputed values for patients on the test treatment are similar to those of patients on the reference arm with similar observed outcomes. The UR arm views withdrawal as a new episode for the patient and does not use their previous exacerbation history. This may be criticised for the reason that data showing poorer outcome on active during treatment may be predictive of poorer outcome than placebo after treatment. In this case, the delta approach to sensitivity analysis imposes proportionally some fixed additional event rate after withdrawal, compared with that before withdrawal. However, here, choice of a specific delta is arbitrary and may be difficult to justify. In any discussion of missing data, it is important to emphasise the need for careful trial design and conduct to limit the impact of missing data. Steps that can be taken to minimise the amount of missing data have been described elsewhere [11,23]. Another important issue is the collection of data after termination of randomised treatment. Such an approach can provide valuable data for assessing the outcome relative to randomised treatment. However, when effective medications are available to such patients, we agree with Mallinckrodt et al. [4] who stated that ‘the confounding influence of rescue medications can render follow-up data of little use in understanding the causal effects of the randomized interventions’. For outcomes other than mortality, it is difficult to obtain complete information, and therefore, there will still be missing data off treatment. The sensitivity analyses described here are therefore still appropriate to assess sensitivity to the remaining missing data in such retrieved dropout analyses. It is important to understand the limitations of the approach outlined in this manuscript. First, it is assumed that the underlying event rate is constant from randomisation until the planned end of trial, although that rate may vary from subject to subject. Second, the implicit variation in rate between subjects is estimated solely from the overdispersion parameter in the negative binomial model. This contrasts with the hierarchical model proposed by Aregay et al. [24], where repeated observations allow the subject to subject variance component to be estimated separate from an overdispersion parameter at the individual count level. This more complex approach could be used by splitting the trial into a series of observation periods. However, we would suggest that in many circumstances the overdispersion does represent subject frailty, making it appropriate to use the simpler model described in this paper.
O. N. Keene et al.
6. CONCLUSION In conclusion, this paper has provided methods for calculating de facto estimands for recurrent event data. There is a role for both de facto and de jure estimands in the analysis of clinical trials. The de facto estimand answers questions such as ‘what would be the effect seen in practice if this treatment were assigned to the target population of eligible patients?’ [3]. The de jure estimand answers a question that an individual patient may have concerning the efficacy of treatment should they take the medication as directed. This estimand is also potentially more useful for determining differences among subgroups and whether a treatment effect persists or reduces over time.
Acknowledgements We are grateful to the reviewers whose constructive comments helped improve the clarity of the paper. SAS is a registered trademark of SAS Institute Inc. in the USA and other countries.
REFERENCES [1] National Research Council. The Prevention and Treatment of Missing Data in Clinical Trials. National Academies Press: Washington, DC, 2010. Available at: http://www.nap.edu/search/?term=12955 (accessed 7 June 2014). [2] Little RJ, D’Agostino R, Cohen ML, Dickersin K, Emerson SS, Farrar JT, Frangakis C, Hogan J, Molenberghs G, Murphy SA, Neaton JD, Rotnitzky A, Scharfstein D, Shih WJ, Siegel JP, Stern H. The prevention and treatment of missing data in clinical trials. New England Journal of Medicine 2012; 367(14):1355–1360. [3] Carpenter JR, Roger JH, Kenward MG. Analysis of longitudinal trials with protocol deviation: a framework for relevant, accessible assumptions, and inference via multiple Imputation. Journal of Biopharmaceutical Statistics 2013; 23:1352–1371. [4] Mallinckrodt CH, Lin Q, Lipkovich I, Molenberghs G. A structured approach to choosing estimands and estimators in longitudinal clinical trials. Pharmaceutical Statistics 2012; 11:456–461. [5] Mallinckrodt C, Roger J, Chuang-Stein C, Molenberghs G, Lane PW, O’Kelly M, Ratitch B, Xu L, Gilbert S, Mehrotra D, Wolfinger R, Thijs H. Missing data: turning guidance into action. Statistics in Biopharmaceutical Research 2013; 5(4):369–382. [6] Keene ON. Intent-to-treat analysis in the presence of off-treatment or missing data. Pharmaceutical Statistics 2011; 10:191–195. [7] Moher D, Hopewell S, Schulz KF, Montori V, Gøtzsche PC, Devereaux PJ, Elbourne D, Egger M, Altman DG. For the CONSORT Group. CONSORT 2010 Explanation and Elaboration: updated guidelines for reporting parallel group randomised trial. BMJ 2010; 340:c869.
[8] Altman DG. Missing outcomes in randomised trials: addressing the dilemma. Open Med 2009; 3(2):51–53. [9] Committee for Medical Products for Human Use. Guideline on missing data in confirmatory clinical trials. European Medicines Evaluation Agency: London, September 2010. Available at: http://www. ema.europa.eu/docs/en_GB/document_library/Scientific_guideline/ 2010/09/WC500096793.pdf (accessed 7 June 2014). [10] Little RJA, Yau L. Intent-to-treat analysis for longitudinal studies with drop-outs. Biometrics 1996; 52:471–483. [11] Mallinckrodt CH. Preventing and Treating Missing Data in Longitudinal Clinical Trials. Cambridge University Press: Cambridge, UK, 2013. [12] O’Kelly M, Ratitch B. Clinical Trials with Missing Data: A Guide for Practitioners. Wiley: Chichester, UK, 2014. [13] Keene ON, Calverley PMA, Jones PW, Vestbo J, Anderson JA. Statistical analysis of exacerbation rates in COPD: TRISTAN and ISOLDE revisited. European Respiratory Journal 2008; 32:17–24. [14] Keene ON, Jones MRK, Lane PW, Anderson J. Analysis of exacerbation rates in asthma and chronic obstructive pulmonary disease: example from the TRISTAN study. Pharmaceutical Statistics 2007; 6:89–97. [15] Pavord ID, Korn S, Howarth P, Bleecker ER, Buhl R, Keene ON, Ortega H, Chanez P. Mepolizumab for severe eosinophilic asthma (DREAM): a multicentre, double-blind, placebo-controlled trial. The Lancet 2012; 380(9842):651–659. [16] Barnard J, Rubin D. Small-sample degrees of freedom with multiple imputation. Biometrika 1999; 86:948–955. [17] Ratitch B, O’Kelly M, Tosiello R. Missing data in clinical trials: from clinical assumptions to statistical analysis using pattern mixture models. Pharmaceutical Statistics 2013; 12:337–347. [18] Meng X-L. Multiple-imputation inferences with uncongenial sources of input. Statistical Science 1994; 9:538–558. [19] Kenward MG, Molenberghs G. Last observation carried forward: a crystal ball? Journal of Biopharmaceutical Statistics 2008; 19:872–888. [20] Carpenter JC, Kenward MG, Evans S, White IR. Last observation carry-forward and last observation analysis. Letter to the Editor. Statistics in Medicine 2004; 23:3241–3244. [21] Molnar FJ, Hutton B, Fergusson G. Does analysis using ‘last observation carried forward’ introduce bias in dementia research? Canadian Medical Association Journal 2008; 179:751–753. [22] Mallinckrodt CH, Lane PW, Schnell D, Peng Y, Mancuso JP. Recommendations for the primary analysis of continuous endpoints in longitudinal clinical trials. Drug Information Journal 2008; 42:303–319. [23] Hughes S, Harris J, Flack N, Cuffe RL. The statistician’s role in the prevention of missing data. Pharmaceutical Statistics 2012; 11:410–416. [24] Aregay M, Shkedy Z, Molenberghs G. A hierarchical Bayesian approach for the analysis of longitudinal count data with overdispersion: a simulation study. Computational Statistics & Data Analysis 2013; 57:233–245.
SUPPORTING INFORMATION Additional supporting information may be found in the online version of this article at the publisher’s web site.
264 Copyright © 2014 John Wiley & Sons, Ltd.
Pharmaceut. Statist. 2014, 13 258–264
Published online 16 June 2014 in Wiley Online Library
Missing data sensitivity analysis for recurrent event data using controlled imputation Oliver N. Keene,a * James H. Roger,b Benjamin F. Hartley,a and Michael G. Kenwardb Statistical analyses of recurrent event data have typically been based on the missing at random assumption. One implication of this is that, if data are collected only when patients are on their randomized treatment, the resulting de jure estimator of treatment effect corresponds to the situation in which the patients adhere to this regime throughout the study. For confirmatory analysis of clinical trials, sensitivity analyses are required to investigate alternative de facto estimands that depart from this assumption. Recent publications have described the use of multiple imputation methods based on pattern mixture models for continuous outcomes, where imputation for the missing data for one treatment arm (e.g. the active arm) is based on the statistical behaviour of outcomes in another arm (e.g. the placebo arm). This has been referred to as controlled imputation or reference-based imputation. In this paper, we use the negative multinomial distribution to apply this approach to analyses of recurrent events and other similar outcomes. The methods are illustrated by a trial in severe asthma where the primary endpoint was rate of exacerbations and the primary analysis was based on the negative binomial model. Copyright © 2014 John Wiley & Sons, Ltd. Keywords: missing; sensitivity; recurrent event; exacerbation; multiple imputation; MNAR
1. INTRODUCTION
258
The potential for bias in study conclusions introduced by missing data has become an important aspect in the interpretation of clinical trials designed to demonstrate the effectiveness of new medicines. The National Research Council of the National Academy of Science report on the prevention and treatment of missing data in clinical trials [1] discussed many of these issues and in particular recommended that ‘examining sensitivity to the assumptions about the missing data mechanism should be a mandatory component of reporting’. In a follow-up paper in the New England Journal of Medicine [2], the authors stated that analysts of clinical trials need to ‘decide on a primary set of assumptions about the missing-data mechanism...In some cases, the primary assumptions can be that data are missing at random’ and to ‘assess the robustness of inferences about treatment effects to various missing-data assumptions by conducting a sensitivity analysis that relates inferences to one or more parameters that capture departures from the primary missing data assumption’. The ‘missing at random’ (MAR) assumption requires that conditional on the data observed for each patient including their covariates in the model (together making their ‘history’), their unavailable data are randomly missing. Or equivalently, for withdrawals, the conditional future statistical behaviour of outcomes given the history is the same for those who remain in the trial and those who withdraw, provided their treatment is the same throughout and they share the same history. This implies that, under MAR, future statistical behaviour for those who withdraw can be modelled from those who remain. It follows that, under MAR, if data are collected only when patients are on their randomised treatment, the resulting estimator of treatment effect
Pharmaceut. Statist. 2014, 13 258–264
corresponds to the situation in which the patients adhere to this regime throughout the study. This has been referred to as the de jure estimand of treatment effect [3] and elsewhere as an efficacy estimand [4], although others view the ‘efficacy/effectiveness terminology’ as confusing [3,5]. A MAR analysis is often labelled as an intention-to-treat (ITT) analysis, but it is actually inconsistent with a strict interpretation of the ITT principle [6]. A full ITT analysis is one that includes all data for a patient regardless of what treatment they take and, if the patient withdraws, must reflect the actual treatment taken following withdrawal [7,8]. Such an analysis including data collected after cessation of randomised treatment would provide a de facto estimand of treatment effect [3] and has been labelled as an effectiveness estimand because it would compare the outcomes of treatment strategies rather than the actual pharmacological effects of treatment itself. An alternative approach to derive such a de facto estimand compares the outcome of treatment strategies by making assumptions about the possible outcome when subjects switch to control, without using the actual off-treatment data that may include the confounding influence of rescue medication received after withdrawal [3,4]. The issue of missing data is a focus of attention for regulatory review of confirmatory trials, and the requirement for sensitivity analyses is clearly stated in the Committee for Medicinal Proda
GlaxoSmithKline Research and Development, Middlesex, UK
b
Medical Statistics Department, London School of Hygiene & Tropical Medicine, London, UK *Correspondence to: Oliver N. Keene, GlaxoSmithKline Research and Development, 1-3 Iron Bridge Road, Uxbridge, Middlesex UB11 1BT, UK. E-mail: [email protected]
Copyright © 2014 John Wiley & Sons, Ltd.
O. N. Keene et al. ucts for Human Use guideline [9], which states, ‘as the choice of primary analysis will be based on assumptions that cannot be verified it will almost always be necessary to investigate the robustness of trial results through appropriate sensitivity analyses that make different assumptions’. Most developments in this area of sensitivity analyses have been for longitudinal continuous outcomes. Little and Yau [10]; Carpenter, Roger and Kenward [3]; Mallinckrodt [11] and O’Kelly and Ratitch [12] describe the use of multiple imputation (MI) methods based on pattern mixture models to investigate the impact of different assumptions about the statistical behaviour of postwithdrawal outcomes. For trials that compare a test and reference product, this approach allows modelling for the postwithdrawal outcome for the test arm to reflect the behaviour of later outcomes on the reference arm. Similar approaches where the postwithdrawal outcomes on the active arm has been chosen to reflect the corresponding observed outcomes on the placebo arm have been called ‘placebo MI’ [4]. The term ‘controlled imputation’ has been used to describe all of these approaches [5]. Recurrent event data occur in clinical trials when multiple events are observed for each patient. Examples of recurrent events include exacerbations in asthma and chronic obstructive pulmonary disease and bone fractures in osteoporosis. Analysis that takes account of the repeated nature of these events is to be preferred to a simple time-to-first-event analysis [13]. A valuable approach that allows simultaneous modelling of the event rate and the variability of this rate between patients is the use of the negative binomial model with log link [14]. This paper describes an extension of these controlled imputation methods for the analysis of recurrent event data. The methods are illustrated by an example in severe asthma.
2. EXAMPLE TRIAL The DREAM trial was a study of mepolizumab in severe asthma [15]. The primary endpoint was the frequency of clinically significant asthma exacerbations, collected over a 1-year treatment
period. The trial randomised subjects equally to four treatment groups: placebo and three doses of mepolizumab. Of the 616 subjects recruited to the trial, 520 (84%) completed treatment. Withdrawal rates were 18% in the placebo group, 16% for mepolizumab 75 mg, 14% for mepolizumab 250 mg and 15% for mepolizumab 750 mg. The time to withdrawal is shown in Figure 1. Reasons for withdrawal of the 96 subjects who withdrew were adverse event (28 subjects), withdrawal of consent (28 subjects), lack of efficacy (22 subjects), lost to follow-up (6 subjects) and other reasons (12 subjects). There was no obvious imbalance in numbers of patients reporting each reason across treatment groups. Each dose of mepolizumab gave a similar reduction in the frequency of exacerbations and a similar level of withdrawal so for simplicity these dose groups have been combined to illustrate the impact of the different methods of analysis. The primary analysis was performed using a negative binomial generalised linear model with a log link function. For this model, as for Poisson regression, the response variable is the total number of exacerbations and time at risk of an exacerbation is fitted as an offset variable. Because the response variable is the total number of events, this model and Poisson regression essentially assume that events occur at a constant rate within an individual and that there is no serial correlation. In the negative binomial model, variation in event rate across patients is accounted for through use of an explicit dispersion parameter in the model. In contrast, for Poisson regression, the use of an overdispersion correction is an alternative way of inflating the variance to allow for this variation from patient to patient. The merits of using the negative binomial to analyse exacerbation rates have been discussed elsewhere [14]. The covariates used in the model were treatment group and other important baseline predictors of outcome. The model included covariates for treatment group, use of maintenance oral corticosteroids, region, number of exacerbations in the previous year and baseline percentage of predicted FEV1 . The log of the
Pharmaceut. Statist. 2014, 13 258–264
259
Figure 1. Time to withdrawal in the DREAM trial.
Copyright © 2014 John Wiley & Sons, Ltd.
O. N. Keene et al. length of time observed in the trial for exacerbations was used as an offset variable. This likelihood-based analysis is valid under the MAR assumption, assuming that the model is correct and estimates the treatment effect that corresponds to what would be observed if those who withdrew had continued to take their randomised treatment throughout the study. That is, it forms a de jure analysis. It is of interest to investigate the sensitivity of the conclusions of this trial to different assumptions regarding the implied behaviour of those that withdrew, that is, to specific departures from the MAR assumption (de facto analysis).
3. TECHNICAL DETAILS OF SENSITIVITY ANALYSIS The negative binomial distribution is conventionally defined as the probability distribution of the number of successes y before k failures are seen in a series of independent Bernoulli trials with probability p of success and .1 p/ of failure. Equation (1) is the resulting probability of y successes.
PŒY D y D
.y C k/ .1 p/k p y .y C 1/.k/
(1)
The distribution has mean D kp=.1 p/ and variance kp=.1 p/2 . The inverse of the parameter k is often referred to as the dispersion parameter as the variance can be reexpressed as . C 2 =k/, noting that the variance for the Poisson distribution is simply . Confusingly, the inverse of the parameter k is sometimes also denoted by k (e.g. in documentation for the SAS PROC GENMOD procedure). For recurrent events, the negative binomial distribution can be derived as a mixture of Poisson distributions where the mixing distribution acts multiplicatively on the intensity œ of the Poisson process to give an effective intensity œ” conditional upon ”, where ” has a gamma distribution with mean k. The negative binomial is the marginal distribution of the number of occurrences after integrating out gamma. The link between œ and p is that p D œ=.1 C œ/ and in reverse œ D p=.1 p/. This recurrent event definition allows for noninteger values of k. In this context, it is often called a gamma-Poisson distribution. The distribution can be extended to cover a series of distinct outcomes. The negative multinomial distribution is the probability distribution of the number of occurrences .y1 , y2 , : : : , ym / of m types of success before we see y0 failures in series of independent trials where probability of outcome j is pj for j D 0, 1, : : : , m and P failure is outcome zero with probability p0 D 1 jD1,:::,m pj . Here, y0 takes the role of the parameter k. The joint probability is given in equation (2). y0 PŒY D . y1 , y2 , : : : , ym / D Qm
P
m jD0 yj
jD0 . yj C 1/
Ym jD0
yj
pj
(2)
260
As with the negative binomial distribution, the negative multinomial distribution can be derived as a series of m independent Poisson values Yj sampled with rates œj ” using the same value
Copyright © 2014 John Wiley & Sons, Ltd.
” from the mixing distribution, with œj D p0 =pj . Then, p0 D 1 P 1 C j 1=œj and pj D p0 =œj . The negative binomial distribution is a special case where m D 1. In this paper, we only use the extension of negative binomial to the negative trinomial as we only consider two periods, before and after withdrawal. The underlying negative multinomial theory can be used to extend this model to a series of predefined discrete periods where the event rate within an individual is assumed constant within each period. For each subject, the period in which withdrawal occurs would then be split into two as we do here. For recurrent events, the mixing parameter ” is a property of the individual patient, and the Yij are the number of events in a series of distinct periods indexed by j for the ith patient. The assumption is that, conditional on the patient’s underlying propensity ”, the numbers of events seen in each period are independent. In particular, there is no serial correlation. The impact of a patient’s propensity works proportionally on the rate predicted by the other components of the model, such as treatment, baseline covariates and length of time at risk. Under this model, the underlying rate of events within a trial might change from period to period and also be modified by covariates. For a patient i with covariates Xi the numbers of events Y1 , : : : , Ym in periods j D 1, : : : , m have independent Poisson distribution with rate: i exp.j C Xi ˇ/. It follows from the gamma-Poisson derivation that the marginal distribution of any subsetP S of the Yj is also nega tive multinomial, with pj D pj = s2S ps C p0 and p0 D P p0 = s2S ps C p0 . So each individual Yi has a marginal distribution which is negative binomial with parameter p D pj =.pj C p0 / and the same k. Similarly, any weighted sum of the Yj also has a negative binomial distribution with parameter p D P P P j2S wj pj = j2S wj pj C p0 j2S wj where k remains the same. The gamma-Poisson model is a special case of a wider class of model. Consider the classic generalised linear mixed model with Poisson distribution and log link. In this model, repeated observations Yit are made on the ith subject in interval t: Yit jsi Poisson.it / and ln.it / D Xit ˇ C si where si is a subject effect with some distribution on the real line. The design matrix Xit describes both the baseline covariates, as well as the imposed treatment structure. Set ui D exp.si /, and it D exp.Xit ˇ/, then: Yit jui Poisson.ui
it /
The margin distribution PŒYit D Ry for Yit is obtained by integrating out the random effect ui , as PŒYit D yjui f .ui /dui The negative binomial is a special case where the distribution of ui is Gamma, while a very similar distribution is generated when si is normally distributed and ui has a log-normal distribution. Clearly, any weighted sum across intervals t for subject i has a P similar distribution within the family with mean t .ui it /. The distribution of Yij conditional upon Yik for the same subject i is
Pharmaceut. Statist. 2014, 13 258–264
O. N. Keene et al. presented in equation (3) where §1 D §ik and §2 D §ij are the event rates for Yik and Yij .
PŒYij D y2 jYik D y1 D
.y1 Cy2 /Š y1 Šy2 Š. 1 C
y1 1 2
y2 2 y / 1Cy2
PŒYij CYik D y2 Cy1 PŒYik D y1
(3) This conditional probability is obtained by noticing that the joint distribution of Yij and Yik can be factored into PŒYij D y2 , Yik D y1 D PŒYij D y2 , Yik D y1 jYij C Yik D y2 C y1 PŒYij C Yik D y2 C y1 (4) where the first term on the right-hand side of (4) describes the allocation of .y1Cy2/ events between the two intervals as a binomial distribution with probabilities §2 =.§1 C §2 / and §1 =.§1 C §2 /. This result relies on the nature of the log link function and is irrespective of the mixing distribution The importance of this result for missing data is that the expression (3) is free of any requirement to integrate out the actual random effect u. That is, the binomial term relies only on the relative sizes of §1 and §2 , while the remaining part involves simple marginal distributions for Yik and .Yij CYik /. This contrasts strongly with the analogous multivariate normal setting, where such a random patient effect is not eliminated through conditioning on the observed part of the sequence. Here, the impact of the distribution of the random effect u on the conditional distribution is summarised through the ratio of the probability distribution of the sum and of the observed first observation. That is, we only require the marginal impact of the distribution of the random effect u on each of these to obtain the marginal impact on the conditional distribution. This means that MAR and MAR-related methods can be implemented simply through these marginal distributions, ignoring the underlying random effects that have been conditioned out. The importance is however limited by the practicality of sampling from the conditional distribution (3) whenever an explicit expression is not available for PŒYij C Yik D y2 C y1 as a function of y2 given y1 . For this reason, we limit consideration to the following special case and ignore the otherwise interesting case where ui has a log-normal distribution. This later might be more easily solved by a data augmentation approach using an MCMC approach. When the mixing distribution is gamma, then both Yij and Yik have negative binomial distributions with probabilities given by equation (1) with the probability p taking values §1 =.1 C §1 / and §2 =.1 C §2 / giving a general form for probability of y events as expressed in terms of § in equation (5): PŒY D y D
.k C y/ y .y C 1/.k/.1 C /kCy
(5)
The conditional distribution of Yij given as Yik D y1 given in equation (6) that can be derived directly from equations (3) and (5) as both Yik and .Yij C Yik / also have negative binomial distributions with § taking values §1 and .§1 C §2 /.
PŒYij D y2 jYik D y1 D
Pharmaceut. Statist. 2014, 13 258–264
(6)
PŒYij D y2 jYik D y1 D
..kCy1 /Cy2 /. .y2 C1/.kCy1 /.1C
y 2 =.1C 1 //2 .kCy1 /Cy2 2 =.1C 1 //
(7) The equivalent result for the negative multinomial has y1 and §1 replaced by the sum of their equivalents in the conditional statement. This property is extremely useful for handling missing data, because the conditional distribution of the missing values for partially observed negative multinomial data is also negative multinomial. Missing data can be imputed using the appropriate negative multinomial models. In its simplest form for any specific subject, conditional upon their own withdrawal time t, the number of events before and after withdrawal can be modelled under a MAR assumption as having a joint negative trinomial distribution, with observed Y1 before withdrawal and unobserved Y2 after withdrawal This is equivalent to the classic primary analysis where Y is assumed to have a negative binomial distribution before withdrawal. We simply extend the model to include the unobserved number of occurrences between withdrawal and study end. Ordinarily, an analysis under a MAR assumption is implemented using direct likelihood, by ignoring the unobserved Y2 and introducing an offset into the log-linear model to allow for the length of the observation period. Using the aforementioned results, an alternative MI approach can be used in which values for Y2 are imputed using the distribution of Y2 conditional on Y1 D y1 . This requires fitting a Bayesian log-linear negative binomial model with offset to the observed data and drawing independent samples from the posterior distribution for parameters k and “ created by multiplying a noninformative prior with the likelihood function. It is assumed for the subsequent MI steps that the priors are sufficiently vague and the sample size large enough, for the data to dominate the posterior. This stage can be implemented simply in SAS using the BAYES statement in the GENMOD procedure. Alternatively, MCMC can be used programming the likelihood directly in either Winbugs or its derivatives or in the MCMC procedure in SAS. For each imputed data set, one sampled set of values k and “ is used to generate a single sample of values y2 for each missing value Y2 conditional upon their associated observed Y1 D y1 . This stage can be programmed using a function to sample from the negative binomial distribution such as the SAS RAND function. The imputed and observed values for each patient are summed to form Y D y1 C y2, and these are analyzed as though the data were complete. For instance, least-square mean treatment differences can be extracted using the GENMOD procedure for each imputation using the BY statement. The resulting estimates and their standard errors are combined across imputations using Rubin’s rules [16]. This last stage can be performed using the MIANALYZE procedure in SAS. This procedure will be valid under MAR and gives results that are very close to those obtained by direct likelihood It is important to note that this approach is not introduced simply as an alternative to using the standard direct likelihood
Copyright © 2014 John Wiley & Sons, Ltd.
261
y2 kCy1 2 .kCy1 Cy2 /.1C 1 / .y2 C1/.kCy1 /.1C§1 C§2 /kCy1Cy2
This conditional distribution is also negative binomial with the k parameter k D k C y1 , while § D §2 =.1 C §1 /, and the mean is §2 .k C y1 /=.1 C §1 / This can be seen by rearranging equation (6) in the form (7) and identifying its relationship to equation (5).
O. N. Keene et al. approach to a MAR analysis. Rather, it provides a convenient framework for sensitivity analysis that can explore the implications of particular MNAR departures from MAR. Sensitivity analyses can be based on the pattern-mixture approach, through modifications of the postwithdrawal imputation models, either of the following: (i) by using a reference-based approach and borrowing parameters from the reference arm or (ii) by using a ‘Delta method’ approach and explicitly altering the rate from the prewithdrawal period [17] In the reference-based approach (i), the rate in the reference arm is used to impute the missing data for both arms. The model for the postwithdrawal outcomes is modified from what would be used under MAR but, most importantly, remains calculable from some fixed linear function of the parameters “ in the original model. Approach (ii) in the current setting would be implemented by multiplying the event rate by a specific ratio, the delta, after withdrawal, equivalent to adding a fixed constant to the linear predictor in the model for §2 . In the sensitivity analysis, the total number of events (observed prewithdrawal + imputed postwithdrawal) is analysed over the full treatment period using the original primary analysis model. For the reference-based approach, possible options include the jump to reference (J2R) and copy reference (CR) approaches analogous to those used for longitudinal normal data (3). The J2R approach represents the situation where the patient’s expected event rate is shifted to that of the reference arm. Such a change may be seen as extreme and might be used as a worst-case scenario in terms of reducing any treatment effect. In this approach, it is assumed that withdrawn patients on the active treatment will quickly lose the positive effect of their treatment when it is stopped. The impact of sampling from the conditional distribution is that if their event rate prior to withdrawal is worse than would be expected (positive residual) on the test treatment; their imputed event rate after withdrawal will be worse than the expected event rate on the reference treatment. Postwithdrawal events in the reference arm are imputed under randomised arm just as under MAR.
In the CR approach, a patient’s expected event rate both before and after withdrawal is assumed to be the same as the reference group. This mimics the case where those withdrawing are in effect nonresponders. This has a less extreme impact than the aforementioned J2R. When a patient on test treatment has more events before withdrawal than expected for them if they had been in the reference arm, then this ‘positive residual’ will feed through into a higher than expected event rate in the postwithdrawal period. This is because their earlier observed higher than normal event rate suggests a personal high propensity to events. Postwithdrawal data in the reference arm are again imputed under randomised arm MAR. An alternative approach based simply on results from the reference arm is unconditional reference (UR). The basis of this approach is that withdrawal from test treatment represents a new episode for the patient, and the previous history of events is not used in the imputation model for events postwithdrawal. Instead, events postwithdrawal for the test arm are imputed using the overall mean for the reference arm, conditional only on baseline covariates. Postwithdrawal data in the reference arm are again imputed under randomised arm MAR. Alternatively, the reference arm could be imputed unconditionally, on the basis of an argument of sharing imputation approaches across arms. This approach may work well when there are other informative covariates other than the previous event history, on which to condition, and make the implicit MAR assumption likely to hold. One of the great advantages of these approaches to sensitivity analysis is that different scenarios can be used for different patients, depending on the reasons for withdrawal. For example, it may be reasonable to use a MAR assumption for patients who withdraw because of moving location and a J2R approach for patients who withdraw for a reason potentially related to treatment such as an adverse event or lack of efficacy. Whichever imputation model is chosen, the subsequent procedure follows that of standard MI. For each imputation, a single set of data is created for each subject’s missing data conditional upon their observed data. After summing observed and imputed
262
Figure 2. Event rate ratios from sensitivity analyses.
Copyright © 2014 John Wiley & Sons, Ltd.
Pharmaceut. Statist. 2014, 13 258–264
O. N. Keene et al. values for each subject, each imputed data set is analysed using the model used for the primary analysis. At least 1000 imputations are recommended and more should be considered, to provide sufficiently stable and reproducible results. We note that when deviations from MAR are included in the imputation model, then the imputation model is what is termed ‘uncongenial’ [18]. Details of how to implement the approach in the SAS system are given in the appendix (available online as Supporting Information).
4. RESULTS OF SENSITIVITY ANALYSIS These sensitivity analyses have been implemented for the example DREAM trial described earlier. Figure 2 shows a Forest plot of the results. The MAR imputation model provides essentially the same results as that for the direct likelihood model. As expected CR, J2R and UR approaches reduce the size of the treatment difference, but the trial remains overwhelmingly positive in terms of size of reduction in exacerbations. For the direct likelihood analysis, the ratio of mepolizumab to placebo is 0.53 (95% CI: 0.43–0.67), while the UR approach gives an estimated ratio of 0.58 (95% CI: 0.46–0.73). The figure also includes an analysis based on reason for withdrawal (labelled on the plot as MAR and J2R imputation). For this analysis, imputation was performed as MAR for patients with a reason viewed as plausibly unrelated to treatment (e.g. lost to follow-up) and J2R for reasons potentially more plausibly related to treatment (e.g. adverse event). This approach provides a result intermediate between the MAR and J2R approaches. The plot also shows the result of one delta method imputation. In this case, the predicted postwithdrawal rate was doubled for active patients compared with active patients remaining in the study, while no delta was applied to placebo patients. This approach provides a result similar to the UR and J2R approaches.
5. DISCUSSION
Pharmaceut. Statist. 2014, 13 258–264
Copyright © 2014 John Wiley & Sons, Ltd.
263
The use of controlled imputation methods for the analysis of recurrent event data allows sensitivity analyses to be performed, which investigate the consequences of departures from the MAR assumption. These departures correspond to different assumptions about postwithdrawal treatments and so map to different scenarios underlying de facto estimators. Alternative approaches to sensitivity analyses for nonnormally distributed data have often relied on filling in the missing data using an ad hoc rule. For example, the Committee for Medicinal Products for Human Use guideline [9] states that for recurrent events, ‘baseline count imputation could also be a plausible approach’. Such an approach is analogous to the baseline observation carried forward (BOCF) approach for longitudinal normally distributed data. These methods are commonly described as ‘single imputation’ and are problematic for several reasons. They are valid only under specific and arguably unrealistic, assumptions [19], and they do not appropriately reflect the variability inherent in the process of imputation. They can mask treatment effects when these exist and create treatment effects when they do not [20]. Molnar, Hutton and Ferguson [21] provide a critique of the use of last observation carried forward (LOCF)-type methods in the neurosciences. Mallinckrodt et al. [22] have stated: ‘Regarding the primary analysis for confirmatory longitudinal clinical trials, conclusive evidence has demonstrated the need to abandon the simple, ad hoc methods such as LOCF and BOCF’. For recurrent
event data, the MI approaches described here are similarly clearly preferable to ad hoc methods. Basing sensitivity analysis on data collected on the reference group (e.g. placebo) has clear attractions as this can reasonably be thought of as what might happen to a patient who discontinues treatment on the test product. However, options remain in how to implement imputations based on the reference group, depending on the nature of the treatment and condition under consideration. For example, does the treatment cause permanent changes or does the treatment have only short-term benefits, as with an analgesic? This paper describes three possibilities: J2R, CR and UR. The J2R imputation method can be thought of as a ‘worst case’ approach; patients who discontinue active medication and who are worse than expected on the test arm are then imputed as worse than expected on the reference arm. More reasonable may be the CR or UR approaches where imputed values for patients on the test treatment are similar to those of patients on the reference arm with similar observed outcomes. The UR arm views withdrawal as a new episode for the patient and does not use their previous exacerbation history. This may be criticised for the reason that data showing poorer outcome on active during treatment may be predictive of poorer outcome than placebo after treatment. In this case, the delta approach to sensitivity analysis imposes proportionally some fixed additional event rate after withdrawal, compared with that before withdrawal. However, here, choice of a specific delta is arbitrary and may be difficult to justify. In any discussion of missing data, it is important to emphasise the need for careful trial design and conduct to limit the impact of missing data. Steps that can be taken to minimise the amount of missing data have been described elsewhere [11,23]. Another important issue is the collection of data after termination of randomised treatment. Such an approach can provide valuable data for assessing the outcome relative to randomised treatment. However, when effective medications are available to such patients, we agree with Mallinckrodt et al. [4] who stated that ‘the confounding influence of rescue medications can render follow-up data of little use in understanding the causal effects of the randomized interventions’. For outcomes other than mortality, it is difficult to obtain complete information, and therefore, there will still be missing data off treatment. The sensitivity analyses described here are therefore still appropriate to assess sensitivity to the remaining missing data in such retrieved dropout analyses. It is important to understand the limitations of the approach outlined in this manuscript. First, it is assumed that the underlying event rate is constant from randomisation until the planned end of trial, although that rate may vary from subject to subject. Second, the implicit variation in rate between subjects is estimated solely from the overdispersion parameter in the negative binomial model. This contrasts with the hierarchical model proposed by Aregay et al. [24], where repeated observations allow the subject to subject variance component to be estimated separate from an overdispersion parameter at the individual count level. This more complex approach could be used by splitting the trial into a series of observation periods. However, we would suggest that in many circumstances the overdispersion does represent subject frailty, making it appropriate to use the simpler model described in this paper.
O. N. Keene et al.
6. CONCLUSION In conclusion, this paper has provided methods for calculating de facto estimands for recurrent event data. There is a role for both de facto and de jure estimands in the analysis of clinical trials. The de facto estimand answers questions such as ‘what would be the effect seen in practice if this treatment were assigned to the target population of eligible patients?’ [3]. The de jure estimand answers a question that an individual patient may have concerning the efficacy of treatment should they take the medication as directed. This estimand is also potentially more useful for determining differences among subgroups and whether a treatment effect persists or reduces over time.
Acknowledgements We are grateful to the reviewers whose constructive comments helped improve the clarity of the paper. SAS is a registered trademark of SAS Institute Inc. in the USA and other countries.
REFERENCES [1] National Research Council. The Prevention and Treatment of Missing Data in Clinical Trials. National Academies Press: Washington, DC, 2010. Available at: http://www.nap.edu/search/?term=12955 (accessed 7 June 2014). [2] Little RJ, D’Agostino R, Cohen ML, Dickersin K, Emerson SS, Farrar JT, Frangakis C, Hogan J, Molenberghs G, Murphy SA, Neaton JD, Rotnitzky A, Scharfstein D, Shih WJ, Siegel JP, Stern H. The prevention and treatment of missing data in clinical trials. New England Journal of Medicine 2012; 367(14):1355–1360. [3] Carpenter JR, Roger JH, Kenward MG. Analysis of longitudinal trials with protocol deviation: a framework for relevant, accessible assumptions, and inference via multiple Imputation. Journal of Biopharmaceutical Statistics 2013; 23:1352–1371. [4] Mallinckrodt CH, Lin Q, Lipkovich I, Molenberghs G. A structured approach to choosing estimands and estimators in longitudinal clinical trials. Pharmaceutical Statistics 2012; 11:456–461. [5] Mallinckrodt C, Roger J, Chuang-Stein C, Molenberghs G, Lane PW, O’Kelly M, Ratitch B, Xu L, Gilbert S, Mehrotra D, Wolfinger R, Thijs H. Missing data: turning guidance into action. Statistics in Biopharmaceutical Research 2013; 5(4):369–382. [6] Keene ON. Intent-to-treat analysis in the presence of off-treatment or missing data. Pharmaceutical Statistics 2011; 10:191–195. [7] Moher D, Hopewell S, Schulz KF, Montori V, Gøtzsche PC, Devereaux PJ, Elbourne D, Egger M, Altman DG. For the CONSORT Group. CONSORT 2010 Explanation and Elaboration: updated guidelines for reporting parallel group randomised trial. BMJ 2010; 340:c869.
[8] Altman DG. Missing outcomes in randomised trials: addressing the dilemma. Open Med 2009; 3(2):51–53. [9] Committee for Medical Products for Human Use. Guideline on missing data in confirmatory clinical trials. European Medicines Evaluation Agency: London, September 2010. Available at: http://www. ema.europa.eu/docs/en_GB/document_library/Scientific_guideline/ 2010/09/WC500096793.pdf (accessed 7 June 2014). [10] Little RJA, Yau L. Intent-to-treat analysis for longitudinal studies with drop-outs. Biometrics 1996; 52:471–483. [11] Mallinckrodt CH. Preventing and Treating Missing Data in Longitudinal Clinical Trials. Cambridge University Press: Cambridge, UK, 2013. [12] O’Kelly M, Ratitch B. Clinical Trials with Missing Data: A Guide for Practitioners. Wiley: Chichester, UK, 2014. [13] Keene ON, Calverley PMA, Jones PW, Vestbo J, Anderson JA. Statistical analysis of exacerbation rates in COPD: TRISTAN and ISOLDE revisited. European Respiratory Journal 2008; 32:17–24. [14] Keene ON, Jones MRK, Lane PW, Anderson J. Analysis of exacerbation rates in asthma and chronic obstructive pulmonary disease: example from the TRISTAN study. Pharmaceutical Statistics 2007; 6:89–97. [15] Pavord ID, Korn S, Howarth P, Bleecker ER, Buhl R, Keene ON, Ortega H, Chanez P. Mepolizumab for severe eosinophilic asthma (DREAM): a multicentre, double-blind, placebo-controlled trial. The Lancet 2012; 380(9842):651–659. [16] Barnard J, Rubin D. Small-sample degrees of freedom with multiple imputation. Biometrika 1999; 86:948–955. [17] Ratitch B, O’Kelly M, Tosiello R. Missing data in clinical trials: from clinical assumptions to statistical analysis using pattern mixture models. Pharmaceutical Statistics 2013; 12:337–347. [18] Meng X-L. Multiple-imputation inferences with uncongenial sources of input. Statistical Science 1994; 9:538–558. [19] Kenward MG, Molenberghs G. Last observation carried forward: a crystal ball? Journal of Biopharmaceutical Statistics 2008; 19:872–888. [20] Carpenter JC, Kenward MG, Evans S, White IR. Last observation carry-forward and last observation analysis. Letter to the Editor. Statistics in Medicine 2004; 23:3241–3244. [21] Molnar FJ, Hutton B, Fergusson G. Does analysis using ‘last observation carried forward’ introduce bias in dementia research? Canadian Medical Association Journal 2008; 179:751–753. [22] Mallinckrodt CH, Lane PW, Schnell D, Peng Y, Mancuso JP. Recommendations for the primary analysis of continuous endpoints in longitudinal clinical trials. Drug Information Journal 2008; 42:303–319. [23] Hughes S, Harris J, Flack N, Cuffe RL. The statistician’s role in the prevention of missing data. Pharmaceutical Statistics 2012; 11:410–416. [24] Aregay M, Shkedy Z, Molenberghs G. A hierarchical Bayesian approach for the analysis of longitudinal count data with overdispersion: a simulation study. Computational Statistics & Data Analysis 2013; 57:233–245.
SUPPORTING INFORMATION Additional supporting information may be found in the online version of this article at the publisher’s web site.
264 Copyright © 2014 John Wiley & Sons, Ltd.
Pharmaceut. Statist. 2014, 13 258–264
