A Bayesian approach to jointly estimate centre and treatment by centre heterogeneity in a proportional hazards model

Posterior cervical foraminotomy versus anterior cervical discectomy for Cervical Brachialgia: the FORVAD RCT

Background

Posterior cervical foraminotomy and anterior cervical discectomy are routinely used operations to treat cervical brachialgia, although definitive evidence supporting superiority of either is lacking.

Objective

The primary objective was to investigate whether or not posterior cervical foraminotomy is superior to anterior cervical discectomy in improving clinical outcome.

Design

This was a Phase III, unblinded, prospective, United Kingdom multicentre, parallel-group, individually randomised controlled superiority trial comparing posterior cervical foraminotomy with anterior cervical discectomy. A rapid qualitative study was conducted during the close-down phase, involving remote semistructured interviews with trial participants and health-care professionals.

Setting

National Health Service trusts.

Participants

Patients with symptomatic unilateral cervical brachialgia for at least 6 weeks.

Interventions

Participants were randomised to receive posterior cervical foraminotomy or anterior cervical discectomy. Allocation was not blinded to participants, medical staff or trial staff. Health-care use from providing the initial surgical intervention to hospital discharge was measured and valued using national cost data.

Main outcome measures

The primary outcome measure was clinical outcome, as measured by patient-reported Neck Disability Index score 52 weeks post operation. Secondary outcome measures included complications, reoperations and restricted American Spinal Injury Association score over 6 weeks post operation, and patient-reported Eating Assessment Tool-10 items, Glasgow-Edinburgh Throat Scale, Voice Handicap Index-10 items, PainDETECT and Numerical Rating Scales for neck and upper-limb pain over 52 weeks post operation.

Results

The target recruitment was 252 participants. Owing to slow accrual, the trial closed after randomising 23 participants from 11 hospitals. The qualitative substudy found that there was support and enthusiasm for the posterior cervical FORaminotomy Versus Anterior cervical Discectomy in the treatment of cervical brachialgia trial and randomised clinical trials in this area. However, clinical equipoise appears to have been an issue for sites and individual surgeons. Randomisation on the day of surgery and processes for screening and approaching participants were also crucial factors in some centres. The median Neck Disability Index scores at baseline (pre surgery) and at 52 weeks was 44.0 (interquartile range 36.0-62.0 weeks) and 25.3 weeks (interquartile range 20.0-42.0 weeks), respectively, in the posterior cervical foraminotomy group (n = 14), and 35.6 weeks (interquartile range 34.0-44.0 weeks) and 45.0 weeks (interquartile range 20.0-57.0 weeks), respectively, in the anterior cervical discectomy group (n = 9). Scores appeared to reduce (i.e. improve) in the posterior cervical foraminotomy group, but not in the anterior cervical discectomy group. The median Eating Assessment Tool-10 items score for swallowing was higher (worse) after anterior cervical discectomy (13.5) than after posterior cervical foraminotomy (0) on day 1, but not at other time points, whereas the median Glasgow-Edinburgh Throat Scale score for globus was higher (worse) after anterior cervical discectomy (15, 7, 6, 6, 2, 2.5) than after posterior cervical foraminotomy (3, 0, 0, 0.5, 0, 0) at all postoperative time points. Five postoperative complications occurred within 6 weeks of surgery, all after anterior cervical discectomy. Neck pain was more severe on day 1 following posterior cervical foraminotomy (Numerical Rating Scale - Neck Pain score 8.5) than at the same time point after anterior cervical discectomy (Numerical Rating Scale - Neck Pain score 7.0). The median health-care costs of providing initial surgical intervention were £2610 for posterior cervical foraminotomy and £4411 for anterior cervical discectomy.

Conclusions

The data suggest that posterior cervical foraminotomy is associated with better outcomes, fewer complications and lower costs, but the trial recruited slowly and closed early. Consequently, the trial is underpowered and definitive conclusions cannot be drawn. Recruitment was impaired by lack of individual equipoise and by concern about randomising on the day of surgery. A large prospective multicentre trial comparing anterior cervical discectomy and posterior cervical foraminotomy in the treatment of cervical brachialgia is still required.

Trial registration

This trial is registered as ISRCTN10133661.

Funding

This project was funded by the National Institute for Health and Care Research (NIHR) Health Technology Assessment programme and will be published in full in Health Technology Assessment; Vol. 27, No. 21. See the NIHR Journals Library website for further project information.

Goodness-of-fit tests in proportional hazards models with random effects

This paper deals with testing the functional form of the covariate effects in a Cox proportional hazards model with random effects. We assume that the responses are clustered and incomplete due to right censoring. The estimation of the model under the null (parametric covariate effect) and the alternative (nonparametric effect) is performed using the full marginal likelihood. Under the alternative, the nonparametric covariate effects are estimated using orthogonal expansions. The test statistic is the likelihood ratio statistic, and its distribution is approximated using a bootstrap method. The performance of the proposed testing procedure is studied through simulations. The method is also applied on two real data sets one from biomedical research and one from veterinary medicine.

Individual participant data meta-analysis of intervention studies with time-to-event outcomes: A review of the methodology and an applied example

Many randomized trials evaluate an intervention effect on time-to-event outcomes. Individual participant data (IPD) from such trials can be obtained and combined in a so-called IPD meta-analysis (IPD-MA), to summarize the overall intervention effect. We performed a narrative literature review to provide an overview of methods for conducting an IPD-MA of randomized intervention studies with a time-to-event outcome. We focused on identifying good methodological practice for modeling frailty of trial participants across trials, modeling heterogeneity of intervention effects, choosing appropriate association measures, dealing with (trial differences in) censoring and follow-up times, and addressing time-varying intervention effects and effect modification (interactions).We discuss how to achieve this using parametric and semi-parametric methods, and describe how to implement these in a one-stage or two-stage IPD-MA framework. We recommend exploring heterogeneity of the effect(s) through interaction and non-linear effects. Random effects should be applied to account for residual heterogeneity of the intervention effect. We provide further recommendations, many of which specific to IPD-MA of time-to-event data from randomized trials examining an intervention effect.We illustrate several key methods in a real IPD-MA, where IPD of 1225 participants from 5 randomized clinical trials were combined to compare the effects of Carbamazepine and Valproate on the incidence of epileptic seizures.

Multilevel model with random effects for clustered survival data with multiple failure outcomes

We present a multilevel frailty model for handling serial dependence and simultaneous heterogeneity in survival data with a multilevel structure attributed to clustering of subjects and the presence of multiple failure outcomes. One commonly observes such data, for example, in multi-institutional, randomized placebo-controlled trials in which patients suffer repeated episodes (eg, recurrent migraines) of the disease outcome being measured. The model extends the proportional hazards model by incorporating a random covariate and unobservable random institution effect to respectively account for treatment-by-institution interaction and institutional variation in the baseline risk. Moreover, a random effect term with correlation structure driven by a first-order autoregressive process is attached to the model to facilitate estimation of between patient heterogeneity and serial dependence. By means of the generalized linear mixed model methodology, the random effects distribution is assumed normal and the residual maximum likelihood and the maximum likelihood methods are extended for estimation of model parameters. Simulation studies are carried out to evaluate the performance of the residual maximum likelihood and the maximum likelihood estimators and to assess the impact of misspecifying random effects distribution on the proposed inference. We demonstrate the practical feasibility of the modeling methodology by analyzing real data from a double-blind randomized multi-institutional clinical trial, designed to examine the effect of rhDNase on the occurrence of respiratory exacerbations among patients with cystic fibrosis.

Bias and precision of methods for estimating the difference in restricted mean survival time from an individual patient data meta-analysis

Background

The difference in restricted mean survival time (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ rmstD\left({t}^{\ast}\right) $$\end{document}rmstDt∗), the area between two survival curves up to time horizon \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {t}^{\ast } $$\end{document}t∗, is often used in cost-effectiveness analyses to estimate the treatment effect in randomized controlled trials. A challenge in individual patient data (IPD) meta-analyses is to account for the trial effect. We aimed at comparing different methods to estimate the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ rmstD\left({t}^{\ast}\right) $$\end{document}rmstDt∗ from an IPD meta-analysis.

Methods

We compared four methods: the area between Kaplan-Meier curves (experimental vs. control arm) ignoring the trial effect (Naïve Kaplan-Meier); the area between Peto curves computed at quintiles of event times (Peto-quintile); the weighted average of the areas between either trial-specific Kaplan-Meier curves (Pooled Kaplan-Meier) or trial-specific exponential curves (Pooled Exponential). In a simulation study, we varied the between-trial heterogeneity for the baseline hazard and for the treatment effect (possibly correlated), the overall treatment effect, the time horizon \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {t}^{\ast } $$\end{document}t∗, the number of trials and of patients, the use of fixed or DerSimonian-Laird random effects model, and the proportionality of hazards. We compared the methods in terms of bias, empirical and average standard errors. We used IPD from the Meta-Analysis of Chemotherapy in Nasopharynx Carcinoma (MAC-NPC) and its updated version MAC-NPC2 for illustration that included respectively 1,975 and 5,028 patients in 11 and 23 comparisons.

Results

The Naïve Kaplan-Meier method was unbiased, whereas the Pooled Exponential and, to a much lesser extent, the Pooled Kaplan-Meier methods showed a bias with non-proportional hazards. The Peto-quintile method underestimated the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ rmstD\left({t}^{\ast}\right) $$\end{document}rmstDt∗, except with non-proportional hazards at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {t}^{\ast } $$\end{document}t∗= 5 years. In the presence of treatment effect heterogeneity, all methods except the Pooled Kaplan-Meier and the Pooled Exponential with DerSimonian-Laird random effects underestimated the standard error of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ rmstD\left({t}^{\ast}\right) $$\end{document}rmstDt∗. Overall, the Pooled Kaplan-Meier method with DerSimonian-Laird random effects formed the best compromise in terms of bias and variance. The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ rmstD\left({t}^{\ast },=,10,\kern0.5em ,\mathrm{years}\right) $$\end{document}rmstDt∗=10years estimated with the Pooled Kaplan-Meier method was 0.49 years (95 % CI: [−0.06;1.03], p = 0.08) when comparing radiotherapy plus chemotherapy vs. radiotherapy alone in the MAC-NPC and 0.59 years (95 % CI: [0.34;0.84], p < 0.0001) in the MAC-NPC2.

Conclusions

We recommend the Pooled Kaplan-Meier method with DerSimonian-Laird random effects to estimate the difference in restricted mean survival time from an individual-patient data meta-analysis.

Electronic supplementary material

The online version of this article (doi:10.1186/s12874-016-0137-z) contains supplementary material, which is available to authorized users.

Semi-parametric proportional hazards models with crossed random effects for psychometric response times

The semi-parametric proportional hazards model with crossed random effects has two important characteristics: it avoids explicit specification of the response time distribution by using semi-parametric models, and it captures heterogeneity that is due to subjects and items. The proposed model has a proportionality parameter for the speed of each test taker, for the time intensity of each item, and for subject or item characteristics of interest. It is shown how all these parameters can be estimated by Markov chain Monte Carlo methods (Gibbs sampling). The performance of the estimation procedure is assessed with simulations and the model is further illustrated with the analysis of response times from a visual recognition task.

Multilevel mixed effects parametric survival models using adaptive Gauss-Hermite quadrature with application to recurrent events and individual participant data meta-analysis

Multilevel mixed effects survival models are used in the analysis of clustered survival data, such as repeated events, multicenter clinical trials, and individual participant data (IPD) meta-analyses, to investigate heterogeneity in baseline risk and covariate effects. In this paper, we extend parametric frailty models including the exponential, Weibull and Gompertz proportional hazards (PH) models and the log logistic, log normal, and generalized gamma accelerated failure time models to allow any number of normally distributed random effects. Furthermore, we extend the flexible parametric survival model of Royston and Parmar, modeled on the log-cumulative hazard scale using restricted cubic splines, to include random effects while also allowing for non-PH (time-dependent effects). Maximum likelihood is used to estimate the models utilizing adaptive or nonadaptive Gauss-Hermite quadrature. The methods are evaluated through simulation studies representing clinically plausible scenarios of a multicenter trial and IPD meta-analysis, showing good performance of the estimation method. The flexible parametric mixed effects model is illustrated using a dataset of patients with kidney disease and repeated times to infection and an IPD meta-analysis of prognostic factor studies in patients with breast cancer. User-friendly Stata software is provided to implement the methods.

The Survival Kit: software to analyze survival data including possibly correlated random effects

The Survival Kit is a Fortran 90 Software intended for survival analysis using proportional hazards models and their extension to frailty models with a single response time. The hazard function is described as the product of a baseline hazard function and a positive (exponential) function of possibly time-dependent fixed and random covariates. Stratified Cox, grouped data and Weibull models can be used. Random effects can be either log-gamma or normally distributed and can account for a pedigree structure. Variance parameters are estimated in a Bayesian context. It is possible to account for the correlated nature of two random effects either by specifying a known correlation coefficient or estimating it from the data. An R interface of the Survival Kit provides a user friendly way to run the software.

Interval estimation of random effects in proportional hazards models with frailties

Semi-parametric frailty models are widely used to analyze clustered survival data. In this article, we propose the use of the hierarchical likelihood interval for individual frailties. We study the relationship between hierarchical likelihood, empirical Bayesian, and fully Bayesian intervals for frailties. We show that our proposed interval can be interpreted as a frequentist confidence interval and Bayesian credible interval under a uniform prior. We also propose an adjustment of the proposed interval to avoid null intervals. Simulation studies show that the proposed interval preserves the nominal confidence level. The procedure is illustrated using data from a multicenter lung cancer clinical trial.

Individual patient data meta-analysis of survival data using Poisson regression models

Background

An Individual Patient Data (IPD) meta-analysis is often considered the gold-standard for synthesising survival data from clinical trials. An IPD meta-analysis can be achieved by either a two-stage or a one-stage approach, depending on whether the trials are analysed separately or simultaneously. A range of one-stage hierarchical Cox models have been previously proposed, but these are known to be computationally intensive and are not currently available in all standard statistical software. We describe an alternative approach using Poisson based Generalised Linear Models (GLMs).

Methods

We illustrate, through application and simulation, the Poisson approach both classically and in a Bayesian framework, in two-stage and one-stage approaches. We outline the benefits of our one-stage approach through extension to modelling treatment-covariate interactions and non-proportional hazards. Ten trials of hypertension treatment, with all-cause death the outcome of interest, are used to apply and assess the approach.

Results

We show that the Poisson approach obtains almost identical estimates to the Cox model, is additionally computationally efficient and directly estimates the baseline hazard. Some downward bias is observed in classical estimates of the heterogeneity in the treatment effect, with improved performance from the Bayesian approach.

Conclusion

Our approach provides a highly flexible and computationally efficient framework, available in all standard statistical software, to the investigation of not only heterogeneity, but the presence of non-proportional hazards and treatment effect modifiers.

Understanding the outcomes of multi-centre clinical trials: a qualitative study of health professional experiences and views

All trials use protocols to standardize practice within and between trial centres and to enable replication of an experiment across space and time. However, while 'centre effects' have been noted in the literature, the processes and mechanisms by which trial staff convert a protocol into practice, and create 'evidence', is a relatively understudied phenomenon. We undertook a qualitative investigation of a multi-centre, UK-based, insulin trial, where differences were found between participating centres in their attainment of the trial's primary clinical endpoint (HbA(1c)), a measure of patients' average blood glucose control. In-depth interviews were conducted with 12 research nurses and nine clinicians recruited from 11 centres in 2009, and explored their views about trial participation and experiences of trial delivery from inception to closeout. Staff accounts highlighted mixed agendas and/or ambivalent views about involvement in pharmaceutically funded trials, and discursive and temporal strategies by which they attempted to separate research from clinical practice and to convert commercially funded work into better patient care. Staff in different centres also reported divergent practices by which they recruited patients into the trial and 'enacted' the protocol to enhance trial outcomes and/or to individualise and improve patient care. By exploring, and comparing, the experiences of staff who worked on the same trial but in different centres, this study highlights the importance of understanding, and exploring, the enactment of protocols in ways which situate individual practices within both local (institutional) and global contexts.

Frailty modelling for survival data from multi-centre clinical trials

Despite the use of standardized protocols in, multi-centre, randomized clinical trials, outcome may vary between centres. Such heterogeneity may alter the interpretation and reporting of the treatment effect. Below, we propose a general frailty modelling approach for investigating, inter alia, putative treatment-by-centre interactions in time-to-event data in multi-centre clinical trials. A correlated random effects model is used to model the baseline risk and the treatment effect across centres. It may be based on shared, individual or correlated random effects. For inference we develop the hierarchical-likelihood (or h-likelihood) approach which facilitates computation of prediction intervals for the random effects with proper precision. We illustrate our methods using disease-free time-to-event data on bladder cancer patients participating in an European Organization for Research and Treatment of Cancer trial, and a simulation study. We also demonstrate model selection using h-likelihood criteria.

Comparing methods to estimate treatment effects on a continuous outcome in multicentre randomized controlled trials: a simulation study

BACKGROUND

Multicentre randomized controlled trials (RCTs) routinely use randomization and analysis stratified by centre to control for differences between centres and to improve precision. No consensus has been reached on how to best analyze correlated continuous outcomes in such settings. Our objective was to investigate the properties of commonly used statistical models at various levels of clustering in the context of multicentre RCTs.

METHODS

Assuming no treatment by centre interaction, we compared six methods (ignoring centre effects, including centres as fixed effects, including centres as random effects, generalized estimating equation (GEE), and fixed- and random-effects centre-level analysis) to analyze continuous outcomes in multicentre RCTs using simulations over a wide spectrum of intraclass correlation (ICC) values, and varying numbers of centres and centre size. The performance of models was evaluated in terms of bias, precision, mean squared error of the point estimator of treatment effect, empirical coverage of the 95% confidence interval, and statistical power of the procedure.

RESULTS

While all methods yielded unbiased estimates of treatment effect, ignoring centres led to inflation of standard error and loss of statistical power when within centre correlation was present. Mixed-effects model was most efficient and attained nominal coverage of 95% and 90% power in almost all scenarios. Fixed-effects model was less precise when the number of centres was large and treatment allocation was subject to chance imbalance within centre. GEE approach underestimated standard error of the treatment effect when the number of centres was small. The two centre-level models led to more variable point estimates and relatively low interval coverage or statistical power depending on whether or not heterogeneity of treatment contrasts was considered in the analysis.

CONCLUSIONS

All six models produced unbiased estimates of treatment effect in the context of multicentre trials. Adjusting for centre as a random intercept led to the most efficient treatment effect estimation across all simulations under the normality assumption, when there was no treatment by centre interaction.

Validation of prognostic indices using the frailty model

A major issue when proposing a new prognostic index is its generalisibility to daily clinical practice. Validation is therefore required. Most validation techniques assess whether "on average" the results obtained by the prognostic index in classifying patients in a new sample of patients are similar to the results obtained in the construction set. We introduce a new important aspect of the generalisibility of a prognostic index: the heterogeneity of the prognostic index risk group hazard ratios over different centers. If substantial variability between centers exists, the prognostic index may have no discriminatory capability in some of the centers. To model such heterogeneity, we use a frailty model including a random center effect and a random prognostic index by center interaction. Statistical inference is based on a Bayesian approach using a Laplacian approximation for the marginal posterior distribution of the variances of the random effects. We investigate different ways to summarize the information available from this marginal posterior distribution. Our approach is applied to a real bladder cancer database for which we demonstrate how to investigate and interpret heterogeneity in prognostic index effect over centers.

Investigating trial and treatment heterogeneity in an individual patient data meta-analysis of survival data by means of the penalized maximum likelihood approach

In a meta-analysis combining survival data from different clinical trials, an important issue is the possible heterogeneity between trials. Such intertrial variation can not only be explained by heterogeneity of treatment effects across trials but also by heterogeneity of their baseline risk. In addition, one might examine the relationship between magnitude of the treatment effect and the underlying risk of the patients in the different trials. Such a scenario can be accounted for by using additive random effects in the Cox model, with a random trial effect and a random treatment-by-trial interaction. We propose to use this kind of model with a general correlation structure for the random effects and to estimate parameters and hazard function using a semi-parametric penalized marginal likelihood method (maximum penalized likelihood estimators). This approach gives smoothed estimates of the hazard function, which represents incidence in epidemiology. The idea for the approach in this paper comes from the study of heterogeneity in a large meta-analysis of randomized trials in patients with head and neck cancers (meta-analysis of chemotherapy in head and neck cancers) and the effect of adding chemotherapy to locoregional treatment. The simulation study and the application demonstrate that the proposed approach yields satisfactory results and they illustrate the need to use a flexible variance-covariance structure for the random effects.

Bayesian statistics in medicine: a 25 year review

This review examines the state of Bayesian thinking as Statistics in Medicine was launched in 1982, reflecting particularly on its applicability and uses in medical research. It then looks at each subsequent five-year epoch, with a focus on papers appearing in Statistics in Medicine, putting these in the context of major developments in Bayesian thinking and computation with reference to important books, landmark meetings and seminal papers. It charts the growth of Bayesian statistics as it is applied to medicine and makes predictions for the future. From sparse beginnings, where Bayesian statistics was barely mentioned, Bayesian statistics has now permeated all the major areas of medical statistics, including clinical trials, epidemiology, meta-analyses and evidence synthesis, spatial modelling, longitudinal modelling, survival modelling, molecular genetics and decision-making in respect of new technologies.