,
Valentijn de Jong [aut] | Title: | Meta-Analysis of Diagnosis and Prognosis Research Studies |
|---|---|
| Description: | Facilitate frequentist and Bayesian meta-analysis of diagnosis and prognosis research studies. It includes functions to summarize multiple estimates of prediction model discrimination and calibration performance (Debray et al., 2019) <doi:10.1177/0962280218785504>. It also includes functions to evaluate funnel plot asymmetry (Debray et al., 2018) <doi:10.1002/jrsm.1266>. Finally, the package provides functions for developing multivariable prediction models from datasets with clustering (de Jong et al., 2021) <doi:10.1002/sim.8981>. |
| Authors: | Thomas Debray [aut, cre] ,
Valentijn de Jong [aut] |
| Maintainer: | Thomas Debray <[email protected]> |
| License: | GPL-3 |
| Version: | 0.4.0.9000 |
| Built: | 2024-11-19 04:57:09 UTC |
| Source: | https://github.com/smartdata-analysis-and-statistics/metamisc |
Facilitate frequentist and Bayesian meta-analysis of diagnosis and prognosis research studies. It includes functions to summarize multiple estimates of prediction model discrimination and calibration performance (Debray et al., 2019) <doi:10.1177/0962280218785504>. It also includes functions to evaluate funnel plot asymmetry (Debray et al., 2018) <doi:10.1002/jrsm.1266>. Finally, the package provides functions for developing multivariable prediction models from datasets with clustering (de Jong et al., 2021) <doi:10.1002/sim.8981>.
The following functionality is currently implemented: univariate meta-analysis of summary data (uvmeta), bivariate meta-analysis of correlated outcomes (riley), meta-analysis of prediction model performance (valmeta), evaluation of funnel plot asymmetry (fat).
The metamisc package also provides a comprehensive framework for developing prediction models when patient-level data from multiple studies or settings are available (metapred).
Thomas Debray <[email protected]>, Valentijn de Jong <[email protected]>
de Jong VMT, Moons KGM, Eijkemans MJC, Riley RD, Debray TPA. Developing more generalizable prediction models from pooled studies and large clustered data sets. Stat Med. 2021;40(15):3533–59.
Debray TPA, Moons KGM, Ahmed I, Koffijberg H, Riley RD. A framework for developing, implementing, and evaluating clinical prediction models in an individual participant data meta-analysis. Stat Med. 2013;32(18):3158–80.
Debray TPA, Damen JAAG, Riley R, Snell KIE, Reitsma JB, Hooft L, et al. A framework for meta-analysis of prediction model studies with binary and time-to-event outcomes. Stat Methods Med Res. 2019 Sep;28(9):2768–86.
Debray TPA, Damen JAAG, Snell KIE, Ensor J, Hooft L, Reitsma JB, et al. A guide to systematic review and meta-analysis of prediction model performance. BMJ. 2017;356:i6460.
Debray TPA, Moons KGM, Riley RD. Detecting small-study effects and funnel plot asymmetry in meta-analysis of survival data: a comparison of new and existing tests. Res Syn Meth. 2018;9(1):41–50.
Riley RD, Moons K, Snell KIE, Ensor J, Hooft L, Altman D, et al. A guide to systematic review and meta-analysis of prognostic factor studies. BMJ. 2019;364:k4597.
Riley RD, Tierney JF, Stewart LA. Individual participant data meta-analysis: a handbook for healthcare research. Hoboken, NJ: Wiley; 2021. ISBN: 978-1-119-33372-2.
Schmid CH, Stijnen T, White IR. Handbook of meta-analysis. First edition. Boca Raton: Taylor and Francis; 2020. ISBN: 978-1-315-11940-3.
Steyerberg EW, Nieboer D, Debray TPA, Van Houwelingen JC. Assessment of heterogeneity in an individual participant data meta-analysis of prediction models: An overview and illustration. Stat Med. 2019;38(22):4290–309.
fat, metapred, riley, uvmeta, valmeta
Function to display autocorrelation of a fitted Bayesian meta-analysis model.
acplot(...)acplot(...)
... |
Additional arguments, which are currently ignored. |
This is a generic function.
A ggplot object.
Thomas Debray <[email protected]>
Function to display autocorrelation of a fitted Bayesian meta-analysis model.
## S3 method for class 'mcmc.list' acplot(x, P, nLags = 50, greek = FALSE, ...)## S3 method for class 'mcmc.list' acplot(x, P, nLags = 50, greek = FALSE, ...)
x |
An object of class |
P |
Optional dataframe describing the parameters to plot and their respective names |
nLags |
Integer indicating the number of lags of the autocorrelation plot. |
greek |
Logical value indicating whether parameter labels have to be parsed to get Greek letters. Defaults to FALSE. |
... |
Additional arguments which are passed to ggs_autocorrelation |
A ggplot object.
Thomas Debray <[email protected]>
Function to display autocorrelation of a fitted Bayesian meta-analysis model.
## S3 method for class 'uvmeta' acplot(x, ...)## S3 method for class 'uvmeta' acplot(x, ...)
x |
An object of class |
... |
Additional arguments which are currently not used |
Results are displayed for the estimated mean (mu) and standard-deviation (tau) of the meta-analysis model.
A ggplot object.
An object of class ggplot
Thomas Debray <[email protected]>
## Not run: data(Roberts) fit <- with(Roberts, uvmeta(r=SDM, r.se=SE, labels=rownames(Roberts), method="BAYES")) acplot(fit) ## End(Not run)## Not run: data(Roberts) fit <- with(Roberts, uvmeta(r=SDM, r.se=SE, labels=rownames(Roberts), method="BAYES")) acplot(fit) ## End(Not run)
Function to display autocorrelation of a fitted Bayesian meta-analysis model.
## S3 method for class 'valmeta' acplot(x, ...)## S3 method for class 'valmeta' acplot(x, ...)
x |
An object of class |
... |
Additional arguments which are currently not used |
Results are displayed for the estimated mean (mu) and standard-deviation (tau) of the meta-analysis model.
A ggplot object.
An object of class ggplot
Thomas Debray <[email protected]>
## Not run: data(EuroSCORE) fit <- valmeta(cstat=c.index, cstat.se=se.c.index, cstat.cilb=c.index.95CIl, cstat.ciub=c.index.95CIu, N=n, O=n.events, data=EuroSCORE, method="BAYES", slab=Study) acplot(fit) ## End(Not run)## Not run: data(EuroSCORE) fit <- valmeta(cstat=c.index, cstat.se=se.c.index, cstat.cilb=c.index.95CIl, cstat.ciub=c.index.95CIu, N=n, O=n.events, data=EuroSCORE, method="BAYES", slab=Study) acplot(fit) ## End(Not run)
The function calculates (transformed versions of) the concordance (c-) statistic with the corresponding sampling variance.
ccalc( cstat, cstat.se, cstat.cilb, cstat.ciub, cstat.cilv, sd.LP, N, O, Po, data, slab, subset, g = NULL, level = 0.95, approx.se.method = 4, ... )ccalc( cstat, cstat.se, cstat.cilb, cstat.ciub, cstat.cilv, sd.LP, N, O, Po, data, slab, subset, g = NULL, level = 0.95, approx.se.method = 4, ... )
cstat |
vector to specify the estimated c-statistics. |
cstat.se |
Optional vector to specify the corresponding standard errors. |
cstat.cilb |
Optional vector to specify the lower limits of the confidence interval. |
cstat.ciub |
Optional vector to specify the upper limits of the confidence interval. |
cstat.cilv |
Optional vector to specify the levels of aformentioned confidence interval limits. (default: 0.95, which corresponds to the 95% confidence interval). |
sd.LP |
Optional vector to specify the standard deviations of the linear predictor (prognostic index). |
N |
Optional vector to specify the sample/group sizes. |
O |
Optional vector to specify the total number of observed events. |
Po |
Optional vector to specify the observed event probabilities. |
data |
Optional data frame containing the variables given to the arguments above. |
slab |
Optional vector with labels for the studies. |
subset |
Optional vector indicating the subset of studies that should be used. This can be a logical vector or a numeric vector indicating the indices of the studies to include. |
g |
a quoted string that is the function to transform estimates of the c-statistic; see the details below. |
level |
Optional numeric to specify the level for the confidence interval, default |
approx.se.method |
integer specifying which method should be used for estimating the standard error of the
c-statistic (Newcombe, 2006). So far, only method |
... |
Additional arguments. |
The c-statistic is a measure of discrimination, and indicates the ability of a prediction model to distinguish between patients developing and not developing the outcome. The c-statistic typically ranges from 0.5 (no discriminative ability) to 1 (perfect discriminative ability).
By default, the function ccalc will derive the c-statistic of each study, together with
the corresponding standard error and 95% confidence interval. However, it is also possible to calculate transformed
versions of the c-statistic. Appropriate standard errors are then derived using the Delta method.
For instance, the logit transformation can be applied by specifying g="log(cstat/(1-cstat))".
For studies where the c-statistic is missing, it is estimated from the standard deviation of the linear predictor
(theta.source="std.dev(LP)"). The corresponding method is described by White et al. (2015).
When missing, the standard error of the c-statistic can be estimated from the confidence interval. Alternatively,
the standard error can be approximated from a combination of the reported c-statistic, the total sample size and
the total number of events (Newcombe, 2006). This can be achieved by adopting (a modification of) the method
proposed by Hanley and McNeil, as specified in approx.se.method.
An object of class c("mm_perf","data.frame") with the following columns:
The (transformed) c-statistics.
Standard errors of the (transformed) c-statistics.
Lower confidence interval of the (transformed) c-statistics. The level is specified in
level. Intervals are calculated on the same scale as theta by assuming a Normal distribution.
Upper confidence interval of the (transformed) c-statistics. The level is specified in
level. Intervals are calculated on the same scale as theta by assuming a Normal distribution.
Method used for calculating the (transformed) c-statistic.
Method used for calculating the standard error of the (transformed) c-statistic.
Thomas Debray <[email protected]>
Debray TPA, Damen JAAG, Snell KIE, Ensor J, Hooft L, Reitsma JB, et al. A guide to systematic review and meta-analysis of prediction model performance. BMJ. 2017;356:i6460.
Debray TPA, Damen JAAG, Riley R, Snell KIE, Reitsma JB, Hooft L, et al. A framework for meta-analysis of prediction model studies with binary and time-to-event outcomes. Stat Methods Med Res. 2018; In press.
Hanley JA, McNeil BJ. The meaning and use of the area under a receiver operating characteristic (ROC) curve. Radiology. 1982; 143(1):29–36.
Newcombe RG. Confidence intervals for an effect size measure based on the Mann-Whitney statistic. Part 2: asymptotic methods and evaluation. Stat Med. 2006; 25(4):559–73.
Snell KI, Ensor J, Debray TP, Moons KG, Riley RD. Meta-analysis of prediction model performance across multiple studies: Which scale helps ensure between-study normality for the C -statistic and calibration measures? Statistical Methods in Medical Research. 2017.
White IR, Rapsomaniki E, the Emerging Risk Factors Collaboration. Covariate-adjusted measures of discrimination for survival data. Biom J. 2015;57(4):592–613.
######### Validation of prediction models with a binary outcome ######### data(EuroSCORE) # Calculate the c-statistic and its standard error est1 <- ccalc(cstat = c.index, cstat.se = se.c.index, cstat.cilb = c.index.95CIl, cstat.ciub = c.index.95CIu, N = n, O = n.events, data = EuroSCORE, slab = Study) est1 # Calculate the logit c-statistic and its standard error est2 <- ccalc(cstat = c.index, cstat.se = se.c.index, cstat.cilb = c.index.95CIl, cstat.ciub = c.index.95CIu, N = n, O = n.events, data = EuroSCORE, slab = Study, g = "log(cstat/(1-cstat))") est2 # Display the results of all studies in a forest plot plot(est1)######### Validation of prediction models with a binary outcome ######### data(EuroSCORE) # Calculate the c-statistic and its standard error est1 <- ccalc(cstat = c.index, cstat.se = se.c.index, cstat.cilb = c.index.95CIl, cstat.ciub = c.index.95CIu, N = n, O = n.events, data = EuroSCORE, slab = Study) est1 # Calculate the logit c-statistic and its standard error est2 <- ccalc(cstat = c.index, cstat.se = se.c.index, cstat.cilb = c.index.95CIl, cstat.ciub = c.index.95CIu, N = n, O = n.events, data = EuroSCORE, slab = Study, g = "log(cstat/(1-cstat))") est2 # Display the results of all studies in a forest plot plot(est1)
A meta-analysis of nine clinical trials investigating the effect of taking diuretics during pregnancy on the risk of pre-eclampsia.
data(Collins)data(Collins)
A data frame with 9 observations on the following 2 variables.
logORa numeric vector with treatment effect sizes (log odds ratio)
SEa numeric vector with the standard error of the treatment effect sizes
Collins, R., Yusuf, S., Peto, R. Overview of randomised trials of diuretics in pregnancy. British Medical Journal 1985, 290, 17–23.
Hardy, R.J. Thompson, S.G. A likelihood approach to meta-analysis with random effects. Statistics in Medicine 1996; 15:619–629.
data(Collins)data(Collins)
Convert a correlation matrix into a covariance matrix
cor2cov(sigma, cormat)cor2cov(sigma, cormat)
sigma |
vector of standard deviations. The order of standard deviations should correspond to the column order in 'cormat'. |
cormat |
a symmetric numeric correlation matrix |
The covariance matrix
Thomas Debray <[email protected]>
Data frame with treatment differences in CD4 cell count.
data("Daniels")data("Daniels")
A data frame with 15 observations on the following 2 variables.
Y1Treatment differences for the log hazard ratio for the development of AIDS or death over 2 years.
vars1Error variances of Y1.
Y2Difference in mean change in CD4 cell count between baseline and 6 month for studies of the AIDS Clinical Trial Group
vars2Error variances of Y2.
The Daniels data comprises 15 phase II/III randomized clinical trials of the HIV Disease Section of the Adult AIDS Clinical Trials Group of the National Institutes of Health, which had data available as of May 1996, which had at least six months of follow-up on some patients and in which at least one patient developed AIDS or died. The data were previously used
by Daniels and Hughes (1997) to assess whether the change in CD4 cell count is a surrogate for time to either development of AIDS or death in drug trials of patients with HIV.
Daniels MJ, Hughes MD. Meta-analysis for the evaluation of potential surrogate markers. Statistics in Medicine 1997; 16: 1965–1982.
Generate a plot of the posterior distribution
dplot(...)dplot(...)
... |
Additional arguments, which are currently ignored. |
This is a generic function.
A ggplot object.
Thomas Debray <[email protected]>
Generate a plot of the posterior distribution
## S3 method for class 'mcmc.list' dplot(x, P, plot_type = "dens", ...)## S3 method for class 'mcmc.list' dplot(x, P, plot_type = "dens", ...)
x |
An object of class |
P |
Optional dataframe describing the parameters to plot and their respective names |
plot_type |
Optional character string to specify whether a density plot ( |
... |
Additional arguments which are currently not used |
A ggplot object.
Thomas Debray <[email protected]>
Function to generate plots for the prior and posterior distribution of a Bayesian meta-analysis.
## S3 method for class 'uvmeta' dplot(x, par, distr_type, plot_type = "dens", ...)## S3 method for class 'uvmeta' dplot(x, par, distr_type, plot_type = "dens", ...)
x |
An object of class |
par |
Character string to specify for which parameter a plot should be generated. Options are |
distr_type |
Character string to specify whether the prior distribution ( |
plot_type |
Character string to specify whether a density plot ( |
... |
Additional arguments which are currently not used |
An object of class ggplot
Thomas Debray <[email protected]>
## Not run: data(Roberts) fit <- with(Roberts, uvmeta(r=SDM, r.se=SE, method="BAYES")) dplot(fit) dplot(fit, distr_type = "posterior") dplot(fit, par = "tau", distr_type = "prior") dplot(fit, plot_type = "hist") ## End(Not run)## Not run: data(Roberts) fit <- with(Roberts, uvmeta(r=SDM, r.se=SE, method="BAYES")) dplot(fit) dplot(fit, distr_type = "posterior") dplot(fit, par = "tau", distr_type = "prior") dplot(fit, plot_type = "hist") ## End(Not run)
Function to generate plots for the prior and posterior distribution of a Bayesian meta-analysis.
## S3 method for class 'valmeta' dplot(x, par, distr_type, plot_type = "dens", ...)## S3 method for class 'valmeta' dplot(x, par, distr_type, plot_type = "dens", ...)
x |
An object of class |
par |
Character string to specify for which parameter a plot should be generated. Options are |
distr_type |
Character string to specify whether the prior distribution ( |
plot_type |
Character string to specify whether a density plot ( |
... |
Additional arguments which are currently not used |
A ggplot object.
An object of class ggplot
Thomas Debray <[email protected]>
## Not run: data(EuroSCORE) # Meta-analysis of the concordance statistic fit <- valmeta(cstat=c.index, cstat.se=se.c.index, cstat.cilb=c.index.95CIl, cstat.ciub=c.index.95CIu, N=n, O=n.events, data=EuroSCORE, method="BAYES", slab=Study) dplot(fit) dplot(fit, distr_type = "posterior") dplot(fit, par = "tau", distr_type = "prior") # Meta-analysis of the O:E ratio EuroSCORE.new <- EuroSCORE EuroSCORE.new$n[c(1, 2, 5, 10, 20)] <- NA pars <- list(hp.tau.dist="dhalft", # Prior for the between-study standard deviation hp.tau.sigma=1.5, # Standard deviation for 'hp.tau.dist' hp.tau.df=3, # Degrees of freedom for 'hp.tau.dist' hp.tau.max=10) # Maximum value for the between-study standard deviation fit2 <- valmeta(measure="OE", O=n.events, E=e.events, N=n, data=EuroSCORE.new, method="BAYES", slab=Study, pars=pars) dplot(fit2, plot_type = "hist") ## End(Not run)## Not run: data(EuroSCORE) # Meta-analysis of the concordance statistic fit <- valmeta(cstat=c.index, cstat.se=se.c.index, cstat.cilb=c.index.95CIl, cstat.ciub=c.index.95CIu, N=n, O=n.events, data=EuroSCORE, method="BAYES", slab=Study) dplot(fit) dplot(fit, distr_type = "posterior") dplot(fit, par = "tau", distr_type = "prior") # Meta-analysis of the O:E ratio EuroSCORE.new <- EuroSCORE EuroSCORE.new$n[c(1, 2, 5, 10, 20)] <- NA pars <- list(hp.tau.dist="dhalft", # Prior for the between-study standard deviation hp.tau.sigma=1.5, # Standard deviation for 'hp.tau.dist' hp.tau.df=3, # Degrees of freedom for 'hp.tau.dist' hp.tau.max=10) # Maximum value for the between-study standard deviation fit2 <- valmeta(measure="OE", O=n.events, E=e.events, N=n, data=EuroSCORE.new, method="BAYES", slab=Study, pars=pars) dplot(fit2, plot_type = "hist") ## End(Not run)
A hypothetical dataset with 500 subjects suspected of having deep vein thrombosis (DVT).
data(DVTipd)data(DVTipd)
A data frame with 500 observations of 16 variables.
sexgender (0=female, 1=male)
malignactive malignancy (0=no active malignancy, 1=active malignancy)
parparesis (0=no paresis, 1=paresis)
surgrecent surgery or bedridden
tendtenderness venous system
oachstoral contraceptives or hst
legentire leg swollen
notraumabsence of leg trauma
calfdif3calf difference >= 3 cm
pitpitting edema
veinvein distension
altdiagnalternative diagnosis present
histdvthistory of previous DVT
ddimdichdichotimized D-dimer value
dvtfinal diagnosis of DVT
studystudy indicator
Hypothetical dataset derived from the Individual Participant Data Meta-Analysis from Geersing et al (2014). The dataset consists of consecutive outpatients with suspected deep vein thrombosis, with documented information on the presence or absence of proximal deep vein thrombosis (dvt) by an acceptable reference test. Acceptable such tests were either compression ultrasonography or venography at initial presentation, or, if venous imaging was not performed, an uneventful follow-up for at least three months.
Geersing GJ, Zuithoff NPA, Kearon C, Anderson DR, Ten Cate-Hoek AJ, Elf JL, et al. Exclusion of deep vein thrombosis using the Wells rule in clinically important subgroups: individual patient data meta-analysis. BMJ. 2014;348:g1340.
data(DVTipd) str(DVTipd) summary(apply(DVTipd,2,as.factor)) ## Develop a prediction model to predict presence of DVT model.dvt <- glm("dvt~sex+oachst+malign+surg+notraum+vein+calfdif3+ddimdich", family=binomial, data=DVTipd) summary(model.dvt)data(DVTipd) str(DVTipd) summary(apply(DVTipd,2,as.factor)) ## Develop a prediction model to predict presence of DVT model.dvt <- glm("dvt~sex+oachst+malign+surg+notraum+vein+calfdif3+ddimdich", family=binomial, data=DVTipd) summary(model.dvt)
Previously published prediction models for predicting the presence of DVT.
data(DVTmodels)data(DVTmodels)
An object of the class litmodels with the following information for each literature model: the study-level descriptives ("descriptives"), the regression coefficient or weight for each predictor ("weights") and the error variance for each regression coefficient or weight ("weights.var").
Previously, several models (Gagne, Oudega) and score charts (Wells, modified Wells, and Hamilton) have been published for evaluating the presence of DVT in suspected patients. These models combine information on mulitple predictors into a weighted sum, that can subsequently be used to obtain estimates of absolute risk. See DVTipd for more information on the predictors.
Wells PS, Anderson DR, Bormanis J, Guy F, Mitchell M, Gray L, Clement C, Robinson KS, Lewandowski B. Value of assessment of pretest probability of deep-vein thrombosis in clinical management. Lancet 1997; 350(9094):1795–1798. DOI: 10.1016/S0140-6736(97)08140-3.
Wells PS, Anderson DR, Rodger M, Forgie M, Kearon C, Dreyer J, Kovacs G, Mitchell M, Lewandowski B, Kovacs MJ. Evaluation of D-dimer in the diagnosis of suspected deep-vein thrombosis. The New England Journal of Medicine 2003; 349(13):1227–1235. DOI: 10.1056/NEJMoa023153.
Gagne P, Simon L, Le Pape F, Bressollette L, Mottier D, Le Gal G. Clinical prediction rule for diagnosing deep vein thrombosis in primary care. La Presse Medicale 2009; 38(4):525–533. DOI: 10.1016/j.lpm.2008.09.022.
Subramaniam RM, Snyder B, Heath R, Tawse F, Sleigh J. Diagnosis of lower limb deep venous thrombosis in emergency department patients: performance of Hamilton and modified Wells scores. Annals of Emergency Medicine 2006; 48(6):678–685. DOI: 10.1016/j.annemergmed.2006.04.010.
Oudega R, Moons KGM, Hoes AW. Ruling out deep venous thrombosis in primary care. a simple diagnostic algorithm including D-dimer testing. Thrombosis and Haemostasis 2005; 94(1):200–205. DOI: 10.1160/TH04-12-0829.
Debray TPA, Koffijberg H, Nieboer D, Vergouwe Y, Steyerberg EW, Moons KGM. Meta-analysis and aggregation of multiple published prediction models. Stat Med. 2014 Jun 30;33(14):2341–62.
Debray TPA, Koffijberg H, Vergouwe Y, Moons KGM, Steyerberg EW. Aggregating published prediction models with individual participant data: a comparison of different approaches. Stat Med. 2012 Oct 15;31(23):2697–712.
data(DVTmodels)data(DVTmodels)
This data set contains estimates on the predictive performance of the European system for cardiac operative risk evaluation (EuroSCORE II) in patients undergoing cardiac surgery. Results are based on the original development study and 22 validations identified by Guida et al.
data("EuroSCORE")data("EuroSCORE")
A data frame with 23 observations on the following 13 variables.
Studya vector with the first author of each validation study
na numeric vector with the total number of patients on which performance estimates are based
n.eventsa numeric vector with the total number of observed events
c.indexa numeric vector with the estimated concordance statistic of each validation
se.c.indexa numeric vector with the standard error of the concordance statistics
c.index.95CIla numeric vector with the lower bound of the 95% confidence interval of the estimated concordance statistics
c.index.95CIua numeric vector with the upper bound of the 95% confidence interval of the estimated concordance statistics
Poa numeric vector with the overall observed event probability of each validation
Pea numeric vector with the overall expected event probability of each validation
SD.Pea numeric vector with the standard error of Pe
e.eventsa numeric vector with the total number of expected events in each validation
multicentrea logical vector describing whether the study was a multicentre study
mean.agea numeric vector describing the mean age of the patients
sd.agea numeric vector with the spread of the age of the patients
pts.before.2010a logical vector describing whether studies included patients before 2010 (i.e., before EuroSCORE II was developed)
Published in 2012, EuroSCORE II was developed using logistic regression in a dataset comprising 16,828 adult patients undergoing major cardiac surgery from 154 hospitals in 43 countries over a 12-week period (May-July) in 2010. EuroSCORE II was developed to predict in-hospital mortality for patients undergoing any type of cardiac surgery. In 2014, a systematic review of published evidence on the performance value of the euroSCORE II was undertaken by Guida et al. Twenty-two validations, including more 145,592 patients from 21 external validation articles (one study included two validations) and a split-sample validation contained within original development article were included in the review; 23 validation studies in total.
Guida P, Mastro F, Scrascia G, Whitlock R, Paparella D. Performance of the European System for Cardiac Operative Risk Evaluation II: a meta-analysis of 22 studies involving 145,592 cardiac surgery procedures. J Thorac Cardiovasc Surg. 2014; 148(6):3049–3057.e1.
Nashef SAM, Roques F, Sharples LD, Nilsson J, Smith C, Goldstone AR, et al. EuroSCORE II. Eur J Cardiothorac Surg. 2012; 41(4):734-744; discussion 744-745.
data(EuroSCORE)data(EuroSCORE)
The presence of small-study effects is a common threat to systematic reviews and meta-analyses, especially when
it is due to publication bias, which occurs when small primary studies are more likely to be reported (published)
if their findings were positive. The presence of small-study effects can be verified by visual inspection of
the funnel plot, where for each included study of the meta-analysis, the estimate of the reported effect size is
depicted against a measure of precision or sample size.
The premise is that the scatter of plots should reflect a funnel shape, if small-study
effects do not exist. However, when small studies are predominately in one direction (usually the
direction of larger effect sizes), asymmetry will ensue.
The fat function implements several tests for detecting funnel plot asymmetry,
which can be used when the presence of between-study heterogeneity in treatment effect is relatively low.
fat(b, b.se, n.total, d.total, d1, d2, method = "E-FIV")fat(b, b.se, n.total, d.total, d1, d2, method = "E-FIV")
b |
Vector with the effect size of each study. Examples are log odds ratio, log hazards ratio, log relative risk. |
b.se |
Optional vector with the standard error of the effect size of each study |
n.total |
Optional vector with the total sample size of each study |
d.total |
Optional vector with the total number of observed events for each study |
d1 |
Optional vector with the total number of observed events in the exposed groups |
d2 |
Optional vector with the total number of observed events in the unexposed groups |
method |
Method for testing funnel plot asymmetry, defaults to |
A common approach to test the presence of small-study effects is to
estimate a regression model where the standardized effect estimate
(effect/SE) is regressed on a measure of precision (1/SE),
(method="E-UW", Egger 1997).
It is possible to allow for between-study heterogeneity by adopting a
multiplicative overdispersion parameter by which the variance in each
study is multiplied (method="E-FIV", Sterne 2000).
Unfortunately, it has been demonstrated that the aforementioned two tests are biased because: (i) the independent variable is subject to sampling variability; (ii) the standardized treatment effect is correlated with its estimated precision; and (iii) for binary data, the independent regression variable is a biased estimate of the true precision, with larger bias for smaller sample sizes (Macaskill et al. 2001).
The standard approach estimates a regression model with the effect size as
a function of the study size (method="M-FIV", Macaskill et al.
2001). Each study is weighted by the precision of the treatment effect
estimate to allow for possible heteroscedasticity.
An alternative approach is to weight each study by a pooled' estimate of the
outcome proportion (method="M-FPV")
For studies with zero events, a continuity correction is applied by adding 0.5 to all cell counts.
This approach (method="P-FPV") estimates a regression model with the
treatment effect as a function of the inverse of the total sample size
(Peters et al. 2006).
For studies with zero events, a continuity correction is applied by adding 0.5 to all cell counts.
This approach was proposed for survival data, and uses the total
number of events as independent variable in the weighted regression model
(Debray et al. 2017). The study weights are based on the inverse variance
(method="D-FIV") or on an approximation thereof
(method="D-FAV").
a list containing the following entries:
A two-sided P-value indicating statistical significance of the funnel plot asymettry test. Values below the significance level (usually defined as 10%) support the presence of funnel plot asymmetry, and thus small-study effects.
A fitted glm object, representing the estimated regression model used for testing funnel
plot asymmetry.
Thomas Debray <[email protected]>
Debray TPA, Moons KGM, Riley RD. Detecting small-study effects and funnel plot asymmetry in meta-analysis of
survival data: a comparison of new and existing tests. Res Syn Meth. 2018;9(1):41–50.
Egger M, Davey Smith G, Schneider M, Minder C. Bias in meta-analysis detected by a simple, graphical test.
BMJ. 1997;315(7109):629–34.
Macaskill P, Walter SD, Irwig L. A comparison of methods to detect publication bias in meta-analysis.
Stat Med. 2001;20(4):641–54.
Peters JL, Sutton AJ, Jones DR, Abrams KR, Rushton L. Comparison of two methods to detect publication bias
in meta-analysis. JAMA. 2006 Feb 8;295(6):676–80.
Sterne JA, Gavaghan D, Egger M. Publication and related bias in meta-analysis: power of statistical tests
and prevalence in the literature. J Clin Epidemiol. 2000;53(11):1119–29.
data(Fibrinogen) b <- log(Fibrinogen$HR) b.se <- ((log(Fibrinogen$HR.975) - log(Fibrinogen$HR.025))/(2*qnorm(0.975))) n.total <- Fibrinogen$N.total d.total <- Fibrinogen$N.events fat(b=b, b.se=b.se) fat(b=b, b.se=b.se, d.total=d.total, method="D-FIV") # Note that many tests are also available via metafor require(metafor) fat(b=b, b.se=b.se, n.total=n.total, method="M-FIV") regtest(x=b, sei=b.se, ni=n.total, model="lm", predictor="ni")data(Fibrinogen) b <- log(Fibrinogen$HR) b.se <- ((log(Fibrinogen$HR.975) - log(Fibrinogen$HR.025))/(2*qnorm(0.975))) n.total <- Fibrinogen$N.total d.total <- Fibrinogen$N.events fat(b=b, b.se=b.se) fat(b=b, b.se=b.se, d.total=d.total, method="D-FIV") # Note that many tests are also available via metafor require(metafor) fat(b=b, b.se=b.se, n.total=n.total, method="M-FIV") regtest(x=b, sei=b.se, ni=n.total, model="lm", predictor="ni")
The Fibrinogen data set is a meta-analysis of 31 studies in which the association between plasma fibrinogen concentration and the risk of coronary heath disease (CHD) was estimated.
data("Fibrinogen")data("Fibrinogen")
A data frame with 5 variables:
N.totala numeric vector describing the total number of patients for each study
N.eventsa numeric vector describing the number of observed events within each study
HRa numeric vector describing the estimated hazard ratio of each study
HR.025a numeric vector describing the lower boundary of the 95% confidence interval of HR
HR.975a numeric vector describing the upper boundary of the 95% confidence interval of HR
Fibrinogen Studies Collaboration. Collaborative meta-analysis of prospective studies of plasma fibrinogen and cardiovascular disease. Eur J Cardiovasc Prev Rehabil. 2004 Feb;11(1):9-17.
Thompson S, Kaptoge S, White I, Wood A, Perry P, Danesh J, et al. Statistical methods for the time-to-event analysis of individual participant data from multiple epidemiological studies. Int J Epidemiol. 2010 Oct;39(5):1345-59.
data(Fibrinogen) ## maybe str(Fibrinogen) ; plot(Fibrinogen) ...data(Fibrinogen) ## maybe str(Fibrinogen) ; plot(Fibrinogen) ...
Extract the fitted values of a metapred object. By default returns fitted values of the model in the
cross-validation procedure.
## S3 method for class 'metapred' fitted( object, select = "cv", step = NULL, model = NULL, as.stratified = TRUE, type = "response", ... )## S3 method for class 'metapred' fitted( object, select = "cv", step = NULL, model = NULL, as.stratified = TRUE, type = "response", ... )
object |
object of class metapred |
select |
character. Select fitted values from "cv" (default) or from "global" model. |
step |
character or numeric. Name or number of step to select if |
model |
character or numeric. Name or number of model to select if |
as.stratified |
logical. |
type |
character. Type of fitted value. |
... |
For compatibility only. |
Function still under development, use with caution.
Only returns type = "response".
Valentijn de Jong
Generate a forest plot by specifying the various effect sizes, confidence intervals and summary estimate.
forest(...)forest(...)
... |
Additional arguments, which are currently ignored. |
This is a generic function. See forest.default for making forest plots of summary statistics,
forest.metapred for plotting metapred objects, and forest.mp.cv.val for plotting
mp.cv.val objects.
Thomas Debray <[email protected]>
Valentijn de Jong <[email protected]>
Generate a forest plot by specifying the various effect sizes, confidence intervals and summary estimate.
## Default S3 method: forest( theta, theta.ci.lb, theta.ci.ub, theta.slab, theta.summary, theta.summary.ci.lb, theta.summary.ci.ub, theta.summary.pi.lb, theta.summary.pi.ub, title, sort = "asc", theme = theme_bw(), predint.linetype = 1, xlim, xlab = "", refline = 0, label.summary = "Summary Estimate", label.predint = "Prediction Interval", nrows.before.summary = 1, study.digits = 2, study.shape = 15, col.diamond = "white", col.predint = "black", size.study = 0.5, size.predint = 1, lty.ref = "dotted", ... )## Default S3 method: forest( theta, theta.ci.lb, theta.ci.ub, theta.slab, theta.summary, theta.summary.ci.lb, theta.summary.ci.ub, theta.summary.pi.lb, theta.summary.pi.ub, title, sort = "asc", theme = theme_bw(), predint.linetype = 1, xlim, xlab = "", refline = 0, label.summary = "Summary Estimate", label.predint = "Prediction Interval", nrows.before.summary = 1, study.digits = 2, study.shape = 15, col.diamond = "white", col.predint = "black", size.study = 0.5, size.predint = 1, lty.ref = "dotted", ... )
theta |
Numeric vector with effect size for each study |
theta.ci.lb |
Numeric vector specifying the lower bound of the confidence interval of the effect sizes |
theta.ci.ub |
Numeric vector specifying the upper bound of the confidence interval of the effect sizes |
theta.slab |
Character vector specifying the study labels |
theta.summary |
Meta-analysis summary estimate of the effect sizes |
theta.summary.ci.lb |
Lower bound of the confidence (or credibility) interval of the summary estimate |
theta.summary.ci.ub |
Upper bound of the confidence (or credibility) interval of the summary estimate |
theta.summary.pi.lb |
Lower bound of the (approximate) prediction interval of the summary estimate. |
theta.summary.pi.ub |
Upper bound of the (approximate) prediction interval of the summary estimate. |
title |
Title of the forest plot |
sort |
By default, studies are sorted by ascending effect size ( |
theme |
Theme to generate the forest plot. By default, the classic dark-on-light ggplot2 theme is used. See ggtheme for more information. |
predint.linetype |
The linetype of the prediction interval |
xlim |
The |
xlab |
Optional character string specifying the X label |
refline |
Optional numeric specifying a reference line |
label.summary |
Optional character string specifying the label for the summary estimate |
label.predint |
Optional character string specifying the label for the (approximate) prediction interval |
nrows.before.summary |
How many empty rows should be introduced between the study results and the summary estimates |
study.digits |
How many significant digits should be used to print the stuy results |
study.shape |
Plotting symbol to use for the study results. By default, a filled square is used. |
col.diamond |
The filling color for the diamond representing the summary estimate. E.g. "red", "blue", or hex color code ("#2e8aff") |
col.predint |
Line color for the prediction interval. E.g. "red", "blue", or hex color code ("#2e8aff") |
size.study |
Line width for the study results in mm |
size.predint |
Line width for the prediction interval in mm |
lty.ref |
Line type for the reference line |
... |
Additional arguments, which are currently ignored. |
An object of class ggplot
Thomas Debray <[email protected]>
Draw a forest plot of the performance of an internally-externally cross-validated model. By default the final model is shown.
## S3 method for class 'metapred' forest(x, perfFUN = 1, step = NULL, method = "REML", model = NULL, ...)## S3 method for class 'metapred' forest(x, perfFUN = 1, step = NULL, method = "REML", model = NULL, ...)
x |
A |
perfFUN |
Numeric or character. Which performance statistic should be plotted? Defaults to the first. |
step |
Which step should be plotted? Defaults to the best step. numeric is converted to name of the step: 0 for an unchanged model, 1 for the first change... |
method |
character string specifying whether a fixed- or a random-effects model should be used to summarize the prediction model performance. A fixed-effects model is fitted when using method="FE". Random-effects models are fitted by setting method equal to one of the following: "DL", "HE", "SJ", "ML", "REML", "EB", "HS", or "GENQ". Default is "REML". |
model |
Which model change should be plotted? NULL (default, best change) or character name of variable or (integer) index of model change. |
... |
Other arguments passed to plotting internals. E.g. |
Valentijn de Jong <[email protected]>
data(DVTipd) # Internal-external cross-validation of a pre-specified model 'f' f <- dvt ~ histdvt + ddimdich + sex + notraum fit <- metapred(DVTipd, strata = "study", formula = f, scope = f, family = binomial) # Display the model's external performance (expressed as mean squared error by default) # for each study forest(fit)data(DVTipd) # Internal-external cross-validation of a pre-specified model 'f' f <- dvt ~ histdvt + ddimdich + sex + notraum fit <- metapred(DVTipd, strata = "study", formula = f, scope = f, family = binomial) # Display the model's external performance (expressed as mean squared error by default) # for each study forest(fit)
Draw a forest plot of the performance of an internally-externally cross-validated model.
## S3 method for class 'mp.cv.val' forest(x, perfFUN = 1, method = "REML", xlab = NULL, ...)## S3 method for class 'mp.cv.val' forest(x, perfFUN = 1, method = "REML", xlab = NULL, ...)
x |
An |
perfFUN |
Numeric or character. Which performance statistic should be plotted? Defaults to the first. |
method |
character string specifying whether a fixed- or a random-effects model should be used to summarize the prediction model performance. A fixed-effects model is fitted when using method="FE". Random-effects models are fitted by setting method equal to one of the following: "DL", "HE", "SJ", "ML", "REML", "EB", "HS", or "GENQ". Default is "REML". |
xlab |
Label on x-axis. Defaults to the name of the performance function. |
... |
Other arguments passed to plotting internals. E.g. |
Valentijn de Jong <[email protected]>
This data set contains estimates on the performance of the Framingham model for predicting coronary heart disease in male populations (Wilson 1998). Results are based on the original development study and 20 validations identified by Damen et al (BMC Med, 2017).
data("Framingham")data("Framingham")
A data frame with 24 observations on the following 19 variables.
AuthorYeara vector describing the study authors
na numeric vector with the total number of patients on which performance estimates are based
n.eventsa numeric vector with the total number of observed events
c.indexa numeric vector with the estimated concordance statistic of each validation
se.c.indexa numeric vector with the standard error of the concordance statistics
c.index.95CIla numeric vector with the lower bound of the 95% confidence interval of the estimated concordance statistics
c.index.95CIua numeric vector with the upper bound of the 95% confidence interval of the estimated concordance statistics
Poa numeric vector with the overall observed event probability of each validation
Pea numeric vector with the overall expected event probability of each validation
t.vala numeric vector describing the time period in which predictive performance was assessed for each validation
mean_agea numeric vector describing the mean age of the patients
sd_agea numeric vector with the spread of the age of the patients
mean_SBPa numeric vector with the mean systolic blood pressure in the validation studies (mm Hg)
sd_SBPa numeric vector with the spread of systolic blood pressure in the validation studies
mean_total_cholesterola numeric vector with the mean total cholesterol in the validation studies (mg/dL)
sd_total_cholesterola numeric vector with the spread of total cholesterol in the validation studies
mean_hdl_cholesterola numeric vector with the mean high-density lipoprotein cholesterol in the validation studies (mg/dL)
sd_hdl_cholesterola numeric vector with the spread of high-density lipoprotein cholesterol in the validation studies
pct_smokera numeric vector with the percentage smokers in the validation studies
The Framingham Risk Score allows physicians to predict 10-year coronary heart disease (CHD) risk in patients without overt CHD. It was developed in 1998 from a middle-aged white population sample, and has subsequently been validated across different populations. The current dataset contains the original (internal validation) results, as well as 23 external validations which were identified through a systematic review. In this review, studies were eligible for inclusion if they described the validation of the original Framingham model and assessed its performance for fatal or nonfatal CHD in males from a general population setting.
Damen JAAG, Hooft L, Schuit E, Debray TPA, Collins GS, Tzoulaki I, et al. Prediction models for cardiovascular disease risk in the general population: systematic review. BMJ. 2016;i2416.
Damen JAAG, Pajouheshnia R, Heus P, Moons KGM, Reitsma JB, Scholten RJPM, et al. Performance of the Framingham risk models and Pooled Cohort Equations: a systematic review and meta-analysis. BMC Med. 2017;17(1):109.
Wilson PW, D'Agostino RB, Levy D, Belanger AM, Silbershatz H, Kannel WB. Prediction of coronary heart disease using risk factor categories. Circulation. 1998; 97(18):1837–47.
data(Framingham)data(Framingham)
This plot shows the evolution of Gelman and Rubin's shrink factor as the number of iterations increases. The code is adapted from the R package coda.
gelmanplot(...)gelmanplot(...)
... |
Additional arguments which are currently not used |
A ggplot object.
An object of class ggplot
Thomas Debray <[email protected]>
This plot shows the evolution of Gelman and Rubin's shrink factor as the number of iterations increases. The code is adapted from the R package coda.
## S3 method for class 'mcmc.list' gelmanplot( x, P, confidence = 0.95, max.bins = 50, autoburnin = TRUE, greek = FALSE, ... )## S3 method for class 'mcmc.list' gelmanplot( x, P, confidence = 0.95, max.bins = 50, autoburnin = TRUE, greek = FALSE, ... )
x |
An mcmc object |
P |
Optional dataframe describing the parameters to plot and their respective names |
confidence |
The coverage probability of the confidence interval for the potential scale reduction factor |
max.bins |
Maximum number of bins, excluding the last one. |
autoburnin |
Logical flag indicating whether only the second half of the series should be used in the computation. If set to TRUE (default) and start(x) is less than end(x)/2 then start of series will be adjusted so that only second half of series is used. |
greek |
Logical value indicating whether parameter labels have to be parsed to get Greek letters. Defaults to false. |
... |
Additional arguments which are currently not used |
A ggplot object.
An object of class ggplot
Thomas Debray <[email protected]>
This plot shows the evolution of Gelman and Rubin's shrink factor as the number of iterations increases. The code is adapted from the R package coda.
## S3 method for class 'uvmeta' gelmanplot(x, confidence = 0.95, ...)## S3 method for class 'uvmeta' gelmanplot(x, confidence = 0.95, ...)
x |
An mcmc object |
confidence |
The coverage probability of the confidence interval for the potential scale reduction factor |
... |
Additional arguments which are currently not used |
A ggplot object.
An object of class ggplot
Thomas Debray <[email protected]>
This plot shows the evolution of Gelman and Rubin's shrink factor as the number of iterations increases. The code is adapted from the R package coda.
## S3 method for class 'valmeta' gelmanplot(x, confidence = 0.95, ...)## S3 method for class 'valmeta' gelmanplot(x, confidence = 0.95, ...)
x |
An mcmc object |
confidence |
The coverage probability of the confidence interval for the potential scale reduction factor |
... |
Additional arguments which are currently not used |
A ggplot object.
An object of class ggplot
Thomas Debray <[email protected]>
## Not run: data(EuroSCORE) # Meta-analysis of the concordance statistic fit <- valmeta(cstat=c.index, cstat.se=se.c.index, cstat.cilb=c.index.95CIl, cstat.ciub=c.index.95CIu, N=n, O=n.events, data=EuroSCORE, method="BAYES", slab=Study) gelmanplot(fit) ## End(Not run)## Not run: data(EuroSCORE) # Meta-analysis of the concordance statistic fit <- valmeta(cstat=c.index, cstat.se=se.c.index, cstat.cilb=c.index.95CIl, cstat.ciub=c.index.95CIu, N=n, O=n.events, data=EuroSCORE, method="BAYES", slab=Study) gelmanplot(fit) ## End(Not run)
Obtain generalizability estimates from a model fit.
gen(object, ...) generalizability(object, ...)gen(object, ...) generalizability(object, ...)
object |
A model fit object, either a metapred or |
... |
By default, the final model is selected. This parameter allows other arguments to be passed to subset.metapred such that the generalizability estimates of other steps/models may be returned.. |
With named values or indices for parameter genFUN one or more estimates
of generalizability can be selected. Use genFUN = 0 to select all.
Valentijn de Jong
The IMPACT dataset comprises 15 studies of patients suffering from traumatic brain injury, including individual patient data from 11 randomized controlled trials and four observational studies.
data("impact")data("impact")
A data frame with 11022 observations on the following 11 variables.
nameName of the study
typeType of study, RCT: randomized controlled trial,OBS: observational cohort
ageAge of the patient
motor_scoreGlasgow Coma Scale motor score
pupilPupillary reactivity
ctMarshall Computerized Tomography classification
hypoxHypoxia (0=no, 1=yes)
hypotsHypotension (0=no, 1=yes)
tsahTraumatic subarachnoid hemorrhage (0=no, 1=yes)
edhEpidural hematoma (0=no, 1=yes)
mort6-month mortality (0=alive, 1=dead)
The included studies were part of the IMPACT project, where a total of 25 prognostic factors were considered for prediction of 6-month mortality. Missing values were imputed using the study as a fixed effect in the imputation model (Steyerberg et al, 2008).
Steyerberg EW, Nieboer D, Debray TPA, Van Houwelingen JC. Assessment of heterogeneity in an individual participant data meta-analysis of prediction models: An overview and illustration. Stat Med. 2019;38(22):4290–309.
Murray GD, Butcher I, McHugh GS, et al. Multivariable prognostic analysis in traumatic brain injury: results from the IMPACT study. J Neurotrauma. 2007;24(2):329–337.
Steyerberg EW, Mushkudiani N, Perel P, et al. Predicting outcome after traumatic brain injury: development and international validation of prognostic scores based on admission characteristics. PLOS Med. 2008;5(8):e165.
data(impact) by(impact, impact$name, summary) # Plot the distribution of age by study library(ggplot2) e <- ggplot(impact, aes(x = name, y = age)) e + geom_violin(aes(fill = type), trim = FALSE) + theme(axis.text.x = element_text(angle = 45)) + xlab("Study")data(impact) by(impact, impact$name, summary) # Plot the distribution of age by study library(ggplot2) e <- ggplot(impact, aes(x = name, y = age)) e + geom_violin(aes(fill = type), trim = FALSE) + theme(axis.text.x = element_text(angle = 45)) + xlab("Study")
This function imputes missing values by their conditional mean
impute_conditional_mean(x, mu, Sigma)impute_conditional_mean(x, mu, Sigma)
x |
A vector with observations, some of which may be missing (indicated by NA) |
mu |
A vector with the population means for 'x'. No missing values are allowed here. |
Sigma |
A matrix describing the population covariance of 'x' |
A vector where missing values for 'x' have been replaced by their conditional mean
Thomas Debray <[email protected]>
# Define the population means mu <- c(0, 1, 2) # Define the covariance of the population Sigma <- diag(1,3) Sigma[1,2] <- Sigma[2,1] <- 0.3 Sigma[2,3] <- Sigma[3,2] <- 0.1 Sigma[1,3] <- Sigma[3,1] <- -0.2 # Generate a 'random' sample from the population that is partially observed x <- c(NA, 2, 4) # Impute the missing values impute_conditional_mean (x=x, mu=mu, Sigma=Sigma)# Define the population means mu <- c(0, 1, 2) # Define the covariance of the population Sigma <- diag(1,3) Sigma[1,2] <- Sigma[2,1] <- 0.3 Sigma[2,3] <- Sigma[3,2] <- 0.1 Sigma[1,3] <- Sigma[3,1] <- -0.2 # Generate a 'random' sample from the population that is partially observed x <- c(NA, 2, 4) # Impute the missing values impute_conditional_mean (x=x, mu=mu, Sigma=Sigma)
Transforms a linear predictor into a probability.
inv.logit(x)inv.logit(x)
x |
A vector of numerics (between -Inf and Inf) |
A vector of numerics between 0 and 1.
Thomas Debray <[email protected]>
Data frame with diagnostic accuracy data from exercise electrocardiography.
data("Kertai")data("Kertai")
One data frame with 4 variables.
integer. number of true positives
integer. number of false negatives
integer. number of false positives
integer. number of true negatives
The Kertai data set is a meta-analysis of prognostic test studies and comprises 7 studies where the diagnostic test accuracy of exercise electrocardiography for predicting cardiac events in patients undergoing major vascular surgery was measured.
Kertai MD, Boersma E, Bax JJ, Heijenbrok-Kal MH, Hunink MGM, L'talien GJ, Roelandt JRTC, van Urk H, Poldermans D. A meta-analysis comparing the prognostic accuracy of six diagnostic tests for predicting perioperative cardiac risk in patients undergoing major vascular surgery. Heart 2003; 89: 1327–1334.
Jackson D, Riley RD, & White IW. Multivariate meta-analysis: Potential and promise. Statistics in Medicine 2010; 30: 2481–2498.
Transforms values between 0 and 1 to values between -Inf and Inf.
logit(x)logit(x)
x |
A vector of numerics (between 0 and 1) |
A vector of numerics (between -Inf and Inf).
Thomas Debray <[email protected]>
This function provides the (restricted) log-likelihood of a fitted model.
## S3 method for class 'riley' logLik(object, ...)## S3 method for class 'riley' logLik(object, ...)
object |
A |
... |
Additional arguments to be passed on to other functions, currently ignored. |
Returns an object of class logLik. This is the (restricted) log-likelihood of the model represented
by object evaluated at the estimated coefficients. It contains at least one attribute,
"df" (degrees of freedom), giving the number of (estimated) parameters in the model.
Thomas Debray <[email protected]>
Riley RD, Thompson JR, Abrams KR. An alternative model for bivariate random-effects meta-analysis when the within-study correlations are unknown. Biostatistics 2008; 9: 172–186.
data(Daniels) fit <- riley(Daniels,control=list(maxit=10000)) logLik(fit)data(Daniels) fit <- riley(Daniels,control=list(maxit=10000)) logLik(fit)
Meta-analysis of the performance or coefficients of a metapred object. Caution: it is still under development.
ma(object, method, ...)ma(object, method, ...)
object |
A model fit object, such as metapred object. |
method |
Character, method for meta-analysis passed to valmeta and uvmeta. Defaults to "REML". |
... |
Produces different object types depending on input.
Valentijn de Jong
Generalized stepwise regression for obtaining a prediction model that is validated with (stepwise) internal-external cross-validation, in or to obtain adequate performance across data sets. Requires data from individuals in multiple studies.
metapred( data, strata, formula, estFUN = "glm", scope = NULL, retest = FALSE, max.steps = 1000, center = FALSE, recal.int = FALSE, cvFUN = NULL, cv.k = NULL, metaFUN = NULL, meta.method = NULL, predFUN = NULL, perfFUN = NULL, genFUN = NULL, selFUN = "which.min", gen.of.perf = "first", ... )metapred( data, strata, formula, estFUN = "glm", scope = NULL, retest = FALSE, max.steps = 1000, center = FALSE, recal.int = FALSE, cvFUN = NULL, cv.k = NULL, metaFUN = NULL, meta.method = NULL, predFUN = NULL, perfFUN = NULL, genFUN = NULL, selFUN = "which.min", gen.of.perf = "first", ... )
data |
data.frame containing the data. Note that |
strata |
Character to specify the name of the strata (e.g. studies or clusters) variable |
formula |
|
estFUN |
Function for estimating the model in the first stage. Currently "lm", "glm" and "logistfirth" are supported. |
scope |
|
retest |
Logical. Should added (removed) terms be retested for removal (addition)? |
max.steps |
Integer. Maximum number of steps (additions or removals of terms) to take. Defaults to 1000, which is essentially as many as it takes. 0 implies no stepwise selection. |
center |
logical. Should numeric predictors be centered around the cluster mean? |
recal.int |
Logical. Should the intercept be recalibrated in each validation? |
cvFUN |
Cross-validation method, on the study (i.e. cluster or stratum) level. "l1o" for leave-one-out cross-validation (default). "bootstrap" for bootstrap. Or "fixed", for one or more data sets which are only used for validation. A user written function may be supplied as well. |
cv.k |
Parameter for cvFUN. For |
metaFUN |
Function for computing the meta-analytic coefficient estimates in two-stage MA.
By default, rma.uni, from the metafor package is used. Default settings are univariate random effects,
estimated with "REML". Method can be passed trough the |
meta.method |
Name of method for meta-analysis. Default is "REML". For more options see rma.uni. |
predFUN |
Function for predicting new values. Defaults to the predicted probability of the outcome, using the link
function of |
perfFUN |
Function for computing the performance of the prediction models. Default: mean squared error
( |
genFUN |
Function or |
selFUN |
Function for selecting the best method. Default: lowest value for |
gen.of.perf |
For which performance measures should generalizability measures be computed? |
... |
To pass arguments to estFUN (e.g. family = "binomial"), or to other FUNctions. |
Use subset.metapred to obtain an individual prediction model from a metapred object.
Note that formula.changes is currently unordered; it does not represent the order of changes in the stepwise
procedure.
metapred is still under development, use with care.
A list of class metapred, containing the final model in global.model, and the stepwise
tree of estimates of the coefficients, performance measures, generalizability measures in stepwise.
Valentijn de Jong <[email protected]>
Debray TPA, Moons KGM, Ahmed I, Koffijberg H, Riley RD. A framework for developing, implementing, and evaluating clinical prediction models in an individual participant data meta-analysis. Stat Med. 2013;32(18):3158-80.
de Jong VMT, Moons KGM, Eijkemans MJC, Riley RD, Debray TPA. Developing more generalizable prediction models from pooled studies and large clustered data sets. Stat Med. 2021;40(15):3533–59.
Riley RD, Tierney JF, Stewart LA. Individual participant data meta-analysis: a handbook for healthcare research. Hoboken, NJ: Wiley; 2021. ISBN: 978-1-119-33372-2.
Schmid CH, Stijnen T, White IR. Handbook of meta-analysis. First edition. Boca Raton: Taylor and Francis; 2020. ISBN: 978-1-315-11940-3.
forest.metapred for generating a forest plot of prediction model performance
data(DVTipd) ## Not run: # Explore heterogeneity in intercept and assocation of 'ddimdich' glmer(dvt ~ 0 + cluster + (ddimdich|study), family = binomial(), data = DVTipd) ## End(Not run) # Scope f <- dvt ~ histdvt + ddimdich + sex + notraum # Internal-external cross-validation of a pre-specified model 'f' fit <- metapred(DVTipd, strata = "study", formula = f, scope = f, family = binomial) fit # Let's try to simplify model 'f' in order to improve its external validity metapred(DVTipd, strata = "study", formula = f, family = binomial) # We can also try to build a generalizable model from scratch ## Not run: # Some additional examples: metapred(DVTipd, strata = "study", formula = dvt ~ 1, scope = f, family = binomial) # Forwards metapred(DVTipd, strata = "study", formula = f, scope = f, family = binomial) # no selection metapred(DVTipd, strata = "study", formula = f, max.steps = 0, family = binomial) # no selection metapred(DVTipd, strata = "study", formula = f, recal.int = TRUE, family = binomial) metapred(DVTipd, strata = "study", formula = f, meta.method = "REML", family = binomial) ## End(Not run) # By default, metapred assumes the first column is the outcome. newdat <- data.frame(dvt=0, histdvt=0, ddimdich=0, sex=1, notraum=0) fitted <- predict(fit, newdata = newdat)data(DVTipd) ## Not run: # Explore heterogeneity in intercept and assocation of 'ddimdich' glmer(dvt ~ 0 + cluster + (ddimdich|study), family = binomial(), data = DVTipd) ## End(Not run) # Scope f <- dvt ~ histdvt + ddimdich + sex + notraum # Internal-external cross-validation of a pre-specified model 'f' fit <- metapred(DVTipd, strata = "study", formula = f, scope = f, family = binomial) fit # Let's try to simplify model 'f' in order to improve its external validity metapred(DVTipd, strata = "study", formula = f, family = binomial) # We can also try to build a generalizable model from scratch ## Not run: # Some additional examples: metapred(DVTipd, strata = "study", formula = dvt ~ 1, scope = f, family = binomial) # Forwards metapred(DVTipd, strata = "study", formula = f, scope = f, family = binomial) # no selection metapred(DVTipd, strata = "study", formula = f, max.steps = 0, family = binomial) # no selection metapred(DVTipd, strata = "study", formula = f, recal.int = TRUE, family = binomial) metapred(DVTipd, strata = "study", formula = f, meta.method = "REML", family = binomial) ## End(Not run) # By default, metapred assumes the first column is the outcome. newdat <- data.frame(dvt=0, histdvt=0, ddimdich=0, sex=1, notraum=0) fitted <- predict(fit, newdata = newdat)
This function calculates (transformed versions of) the ratio of total number of observed versus expected events with the corresponding sampling variance.
oecalc( OE, OE.se, OE.cilb, OE.ciub, OE.cilv, EO, EO.se, citl, citl.se, N, O, E, Po, Po.se, Pe, data, slab, add = 1/2, g = NULL, level = 0.95, ... )oecalc( OE, OE.se, OE.cilb, OE.ciub, OE.cilv, EO, EO.se, citl, citl.se, N, O, E, Po, Po.se, Pe, data, slab, add = 1/2, g = NULL, level = 0.95, ... )
OE |
vector with the estimated ratio of total observed versus total expected events |
OE.se |
Optional vector with the standard errors of the estimated O:E ratios. |
OE.cilb |
Optional vector to specify the lower limits of the confidence interval for |
OE.ciub |
Optional vector to specify the upper limits of the confidence interval for |
OE.cilv |
Optional vector to specify the levels of aformentioned confidence interval limits. (default: 0.95, which corresponds to the 95% confidence interval). |
EO |
Optional vector with the estimated ratio of total expected versus total observed events |
EO.se |
Optional vector with the standard errors of the estimated E:O ratios |
citl |
Optional vector with the estimated calibration-in-the-large statistics |
citl.se |
Optional vector with the standard error of the calibration-in-the-large statistics |
N |
Optional vector to specify the sample/group sizes. |
O |
Optional vector to specify the total number of observed events. |
E |
Optional vector to specify the total number of expected events |
Po |
Optional vector to specify the (cumulative) observed event probabilities. |
Po.se |
Optional vector with the standard errors of |
Pe |
Optional vector to specify the (cumulative) expected event probabilites
(if specified, during time |
data |
Optional data frame containing the variables given to the arguments above. |
slab |
Optional vector with labels for the studies. |
add |
a non-negative number indicating the amount to add to zero counts. See ‘Details’ |
g |
a quoted string that is the function to transform estimates of the total O:E ratio; see the details below. |
level |
level for confidence interval, default |
... |
Additional arguments. |
An object of class c("mm_perf","data.frame") with the following columns:
The (transformed) O:E ratio.
Standard errors of the (transformed) O:E ratio.
Lower confidence interval of the (transformed) O:E ratios. The level is specified in
level. Intervals are calculated on the same scale as theta by assuming a Normal distribution.
Upper confidence interval of the (transformed) c-statistics. The level is specified in
level. Intervals are calculated on the same scale as theta by assuming a Normal distribution.
Method used for calculating the (transformed) O:E ratio.
Method used for calculating the standard error of the (transformed) O:E ratio.
Thomas Debray <[email protected]>
Debray TPA, Damen JAAG, Snell KIE, Ensor J, Hooft L, Reitsma JB, et al. A guide to systematic review and meta-analysis of prediction model performance. BMJ. 2017;356:i6460.
Debray TPA, Damen JAAG, Riley R, Snell KIE, Reitsma JB, Hooft L, et al. A framework for meta-analysis of prediction model studies with binary and time-to-event outcomes. Stat Methods Med Res. 2019 Sep;28(9):2768–86.
Snell KI, Ensor J, Debray TP, Moons KG, Riley RD. Meta-analysis of prediction model performance across multiple studies: Which scale helps ensure between-study normality for the C -statistic and calibration measures? Stat Methods Med Res. 2017.
######### Validation of prediction models with a binary outcome ######### data(EuroSCORE) # Calculate the total O:E ratio and its standard error est1 <- oecalc(O = n.events, E = e.events, N = n, data = EuroSCORE, slab = Study) est1 # Calculate the log of the total O:E ratio and its standard error est2 <- oecalc(O = n.events, E = e.events, N = n, data = EuroSCORE, slab = Study, g = "log(OE)") est2 # Display the results of all studies in a forest plot plot(est1)######### Validation of prediction models with a binary outcome ######### data(EuroSCORE) # Calculate the total O:E ratio and its standard error est1 <- oecalc(O = n.events, E = e.events, N = n, data = EuroSCORE, slab = Study) est1 # Calculate the log of the total O:E ratio and its standard error est2 <- oecalc(O = n.events, E = e.events, N = n, data = EuroSCORE, slab = Study, g = "log(OE)") est2 # Display the results of all studies in a forest plot plot(est1)
Obtain performance estimates from a model fit.
perf(object, ...) performance(object, ...)perf(object, ...) performance(object, ...)
object |
A model fit object, either a metapred or |
... |
By default, the final model is selected. This parameter allows other arguments to be
passed to subset.metapred such that the performance estimates of other steps/models may be
returned. Use |
Valentijn de Jong
Generates a funnel plot for a fitted fat object.
## S3 method for class 'fat' plot( x, ref, confint = TRUE, confint.level = 0.1, confint.col = "skyblue", confint.alpha = 0.5, confint.density = NULL, xlab = "Effect size", add.pval = TRUE, ... )## S3 method for class 'fat' plot( x, ref, confint = TRUE, confint.level = 0.1, confint.col = "skyblue", confint.alpha = 0.5, confint.density = NULL, xlab = "Effect size", add.pval = TRUE, ... )
x |
An object of class |
ref |
A numeric value indicating the fixed or random effects summary estimate. If no value is provided then it will be retrieved from a fixed effects meta-analysis (if possible). |
confint |
A logical indicator. If |
confint.level |
Significance level for constructing the confidence interval. |
confint.col |
The color for filling the confidence interval. Choose |
confint.alpha |
A numeric value between 0 and 1 indicating the opacity for the confidence region. |
confint.density |
The density of shading lines, in lines per inch. The default value of |
xlab |
A title for the x axis |
add.pval |
Logical to indicate whether a P-value should be added to the plot |
... |
Additional arguments. |
Thomas Debray <[email protected]>
Frantisek Bartos <[email protected]>
data(Fibrinogen) b <- log(Fibrinogen$HR) b.se <- ((log(Fibrinogen$HR.975) - log(Fibrinogen$HR.025))/(2*qnorm(0.975))) n.total <- Fibrinogen$N.total # A very simple funnel plot plot(fat(b=b, b.se=b.se), xlab = "Log hazard ratio") # Plot the funnel for an alternative test plot(fat(b=b, b.se=b.se, n.total=n.total, method="M-FIV"), xlab = "Log hazard ratio")data(Fibrinogen) b <- log(Fibrinogen$HR) b.se <- ((log(Fibrinogen$HR.975) - log(Fibrinogen$HR.025))/(2*qnorm(0.975))) n.total <- Fibrinogen$N.total # A very simple funnel plot plot(fat(b=b, b.se=b.se), xlab = "Log hazard ratio") # Plot the funnel for an alternative test plot(fat(b=b, b.se=b.se, n.total=n.total, method="M-FIV"), xlab = "Log hazard ratio")
Function to create forest plots for objects of class "mm_perf".
## S3 method for class 'mm_perf' plot(x, ...)## S3 method for class 'mm_perf' plot(x, ...)
x |
An object of class |
... |
Additional arguments which are passed to forest. |
The forest plot shows the performance estimates of each study with corresponding confidence intervals.
An object of class ggplot
Thomas Debray <[email protected]>
Lewis S, Clarke M. Forest plots: trying to see the wood and the trees. BMJ. 2001; 322(7300):1479–80.
data(EuroSCORE) # Calculate the c-statistic and its standard error est1 <- ccalc(cstat = c.index, cstat.se = se.c.index, cstat.cilb = c.index.95CIl, cstat.ciub = c.index.95CIu, N = n, O = n.events, data = EuroSCORE, slab = Study) plot(est1) # Calculate the total O:E ratio and its standard error est2 <- oecalc(O = n.events, E = e.events, N = n, data = EuroSCORE, slab = Study) plot(est2)data(EuroSCORE) # Calculate the c-statistic and its standard error est1 <- ccalc(cstat = c.index, cstat.se = se.c.index, cstat.cilb = c.index.95CIl, cstat.ciub = c.index.95CIu, N = n, O = n.events, data = EuroSCORE, slab = Study) plot(est1) # Calculate the total O:E ratio and its standard error est2 <- oecalc(O = n.events, E = e.events, N = n, data = EuroSCORE, slab = Study) plot(est2)
Generates a forest plot for each outcome of the bivariate meta-analysis.
## S3 method for class 'riley' plot(x, title, sort = "asc", xlim, refline, ...)## S3 method for class 'riley' plot(x, title, sort = "asc", xlim, refline, ...)
x |
An object of class |
title |
Title of the forest plot |
sort |
By default, studies are ordered by ascending effect size ( |
xlim |
The |
refline |
Optional numeric specifying a reference line |
... |
Additional parameters for generating forest plots |
Thomas Debray <[email protected]>
Riley RD, Thompson JR, Abrams KR. An alternative model for bivariate random-effects meta-analysis when the within-study correlations are unknown. Biostatistics 2008; 9: 172–186.
data(Scheidler) #Perform the analysis fit <- riley(Scheidler[which(Scheidler$modality==1),]) plot(fit) require(ggplot2) plot(fit, sort="none", theme=theme_gray())data(Scheidler) #Perform the analysis fit <- riley(Scheidler[which(Scheidler$modality==1),]) plot(fit) require(ggplot2) plot(fit, sort="none", theme=theme_gray())
Function to create forest plots for objects of class "uvmeta".
## S3 method for class 'uvmeta' plot(x, sort = "asc", ...)## S3 method for class 'uvmeta' plot(x, sort = "asc", ...)
x |
An object of class |
sort |
By default, studies are ordered by ascending effect size ( |
... |
Additional arguments which are passed to forest. |
The forest plot shows the performance estimates of each validation with corresponding confidence intervals. A polygon is added to the bottom of the forest plot, showing the summary estimate based on the model. A 95% prediction interval is added by default for random-effects models, the dotted line indicates its (approximate) bounds
Full lines indicate confidence intervals or credibility intervals (in case of a Bayesian meta-analysis). Dashed lines indicate prediction intervals. The width of all intervals is defined by the significance level chosen during meta-analysis.
Thomas Debray <[email protected]>
Lewis S, Clarke M. Forest plots: trying to see the wood and the trees. BMJ. 2001; 322(7300):1479–80.
Riley RD, Higgins JPT, Deeks JJ. Interpretation of random effects meta-analyses. BMJ. 2011 342:d549–d549.
data(Roberts) # Frequentist random-effects meta-analysis fit <- with(Roberts, uvmeta(r=SDM, r.se=SE, labels=rownames(Roberts))) plot(fit)data(Roberts) # Frequentist random-effects meta-analysis fit <- with(Roberts, uvmeta(r=SDM, r.se=SE, labels=rownames(Roberts))) plot(fit)
Function to create forest plots for objects of class "valmeta".
## S3 method for class 'valmeta' plot(x, ...)## S3 method for class 'valmeta' plot(x, ...)
x |
An object of class |
... |
Additional arguments which are passed to forest. |
The forest plot shows the performance estimates of each validation with corresponding confidence intervals. A polygon is added to the bottom of the forest plot, showing the summary estimate based on the model. A 95% prediction interval is added by default for random-effects models, the dotted line indicates its (approximate) bounds.
An object of class ggplot
Thomas Debray <[email protected]>
Debray TPA, Damen JAAG, Snell KIE, Ensor J, Hooft L, Reitsma JB, et al. A guide to systematic review and meta-analysis of prediction model performance. BMJ. 2017;356:i6460.
Lewis S, Clarke M. Forest plots: trying to see the wood and the trees. BMJ. 2001; 322(7300):1479–80.
Riley RD, Higgins JPT, Deeks JJ. Interpretation of random effects meta-analyses. BMJ. 2011 342:d549–d549.
When a Bayesian meta-analysis was conducted, the prior and posterior distribution can be visualized using dplot.valmeta.
Further, the running means and the presence of autocorrelation can be inspected using rmplot.valmeta and, respectively,
acplot.valmeta.
data(EuroSCORE) fit <- valmeta(cstat=c.index, cstat.se=se.c.index, cstat.cilb=c.index.95CIl, cstat.ciub=c.index.95CIu, N=n, O=n.events, data=EuroSCORE) plot(fit) library(ggplot2) plot(fit, theme=theme_grey())data(EuroSCORE) fit <- valmeta(cstat=c.index, cstat.se=se.c.index, cstat.cilb=c.index.95CIl, cstat.ciub=c.index.95CIu, N=n, O=n.events, data=EuroSCORE) plot(fit) library(ggplot2) plot(fit, theme=theme_grey())
Calculates a prediction interval for the summary parameters of Riley's alternative model for bivariate random-effects meta-analysis. This interval predicts in what range future observations will fall given what has already been observed.
## S3 method for class 'riley' predict(object, ...)## S3 method for class 'riley' predict(object, ...)
object |
A |
... |
Additional arguments (currently ignored) |
Prediction intervals are based on Student's t-distribution with (numstudies - 5) degrees of freedom. The width of the interval is specified by the significance level chosen during meta-analysis.
Data frame containing prediction intervals with the summary estimates beta1 and beta2
(for effect size data), or with the mean sensitivity and false positive rate (for diagnostic test accuracy data).
Thomas Debray <[email protected]>
recalibrate is used to recalibrate a prediction model of classes metapred, glm or lm.
recalibrate(object, newdata, f = ~1, estFUN = NULL, ...)recalibrate(object, newdata, f = ~1, estFUN = NULL, ...)
object |
A model fit object to be recalibrated, of class |
newdata |
data.frame containing new data set for updating. |
f |
formula. Which coefficients of the model should be updated? Default: intercept only. Left-hand side may be left out. See formula for details. |
estFUN |
Function for model estimation. If left |
... |
Optional arguments to pass to |
Currently only the coefficients are updated and the variances and other aspects are left untouched. For updating the entire model and all its statistics, see update.
Recalibrated model fit object, of the same class as object. Generally, updated coefficients can
be retrieved with coef().
data(DVTipd) DVTipd$cluster <- 1:4 # Add a fictional clustering to the data set. # Suppose we estimated the model in three studies: DVTipd123 <- DVTipd[DVTipd$cluster <= 3, ] mp <- metamisc:::metapred(DVTipd123, strata = "cluster", f = dvt ~ vein + malign, family = binomial) # and now want to recalibrate it for the fourth: DVTipd4 <- DVTipd[DVTipd$cluster == 4, ] metamisc:::recalibrate(mp, newdata = DVTipd4)data(DVTipd) DVTipd$cluster <- 1:4 # Add a fictional clustering to the data set. # Suppose we estimated the model in three studies: DVTipd123 <- DVTipd[DVTipd$cluster <= 3, ] mp <- metamisc:::metapred(DVTipd123, strata = "cluster", f = dvt ~ vein + malign, family = binomial) # and now want to recalibrate it for the fourth: DVTipd4 <- DVTipd[DVTipd$cluster == 4, ] metamisc:::recalibrate(mp, newdata = DVTipd4)
This function fits the alternative model for bivariate random-effects meta-analysis when the within-study correlations
are unknown. This bivariate model was proposed by Riley et al. (2008) and is similar to the general bivariate
random-effects model (van Houwelingen et al. 2002), but includes an overall correlation parameter rather than
separating the (usually unknown) within- and between-study correlation. As a consequence, the alternative model
is not fully hierarchical, and estimates of additional variation beyond sampling error (psi) are not
directly equivalent to the between-study variation (tau) from the general model. This model is particularly
useful when there is large within-study variability, few primary studies are available or the general model
estimates the between-study correlation as 1 or -1.
riley(X, slab, optimization = "Nelder-Mead", control = list(), pars, ...)riley(X, slab, optimization = "Nelder-Mead", control = list(), pars, ...)
X |
data frame containing integer variables |
slab |
Optional vector specifying the label for each study |
optimization |
The optimization method that should be used for minimizing the negative (restricted) log-likelihood function. The default method is an implementation of that of Nelder and Mead (1965), that uses only function values and is robust but relatively slow. Other methods are described in optim. |
control |
A list of control parameters to pass to optim. |
pars |
List with additional arguments. The width of confidence, credibility and prediction intervals is
defined by |
... |
Arguments to be passed on to other functions. See "Details" for more information. |
Parameters are estimated by iteratively maximizing the restriced log-likelihood using the Newton-Raphson procedure.
The results from a univariate random-effects meta-analysis with a method-of-moments estimator are used as starting
values for beta1, beta2, psi1 and psi2 in the optim command. Standard errors for all parameters are obtained
from the inverse Hessian matrix.
The following parameters are estimated by iteratively maximizing the restriced log-likelihood using the Newton-Raphson
procedure: pooled effect size for outcome 1 (beta1), pooled effect size for outcome 2 (beta2),
additional variation of beta1 beyond sampling error (psi1), additional variation of beta2
beyond sampling error (psi2) and the correlation rho between psi1 and psi2.
Although the model can also be used for diagnostic test accuracy data when substantial within-study correlations are expected, assuming zero within-study correlations (i.e. applying Reitsma's approach) is usually justified (Reitsma et al. 2005, Daniels and Hughes 1997, Korn et al. 2005, Thompson et al. 2005, Van Houwelingen et al. 2002).
A logit transformation is applied to the sensitivities ans false positive rates of each study, in order to meet the normality
assumptions. When zero cell counts occur, continuity corrections may be required. The correction value can be specified using
correction (defaults to 0.5). Further, when the argument correction.control is set to "all"
(the default) the continuity correction is added to the whole data if only one cell in one study is zero.
If correction.control="single" the correction is only applied to rows of the data which have a zero.
The following parameters are estimated: logit of sensitivity (beta1), logit of false positive rate (beta2),
additional variation of beta1 beyond sampling error (psi1), additional variation of beta2 beyond
sampling error (psi2) and the correlation (rho) between psi1 and psi2.
An object of the class riley for which many standard methods are available. A warning message is
casted when the Hessian matrix contains negative eigenvalues, which implies that the identified solution is a
saddle point and thus not optimal.
The overall correlation parameter rho is transformed during estimation to ensure that corresponding values
remain between -1 and 1. The transformed correlation rhoT is given as logit((rho+1)/2). During optimization,
the starting value for rhoT is set to 0. The standard error of rho is derived from rhoT using
the Delta method. Similarly, the Delta methods is used to derive the standard error of the sensitivity and false
positive rate from beta1 and, respectively, beta2.
Algorithms for dealing with missing data are currently not implemented, but Bayesian approaches will become available in later versions.
Thomas Debray <[email protected]>
Korn EL, Albert PS, McShane LM. Assessing surrogates as trial endpoints using mixed models. Statistics in Medicine 2005; 24: 163–182.
Nelder JA, Mead R. A simplex algorithm for function minimization. Computer Journal (1965); 7: 308–313.
Reitsma J, Glas A, Rutjes A, Scholten R, Bossuyt P, Zwinderman A. Bivariate analysis of sensitivity and specificity produces informative summary measures in diagnostic reviews. Journal of Clinical Epidemiology 2005; 58: 982–990.
Riley RD, Thompson JR, Abrams KR. An alternative model for bivariate random-effects meta-analysis when the within-study correlations are unknown. Biostatistics 2008; 9: 172–186.
Thompson JR, Minelli C, Abrams KR, Tobin MD, Riley RD. Meta-analysis of genetic studies using mendelian randomization–a multivariate approach. Statistics in Medicine 2005; 24: 2241–2254.
van Houwelingen HC, Arends LR, Stijnen T. Advanced methods in meta-analysis: multivariate approach and meta-regression. Statistics in Medicine 2002; 21: 589–624.
data(Scheidler) data(Daniels) data(Kertai) # Meta-analysis of potential surrogate markers data # The results obtained by Riley (2008) were as follows: # beta1 = -0.042 (SE = 0.063), # beta2 = 14.072 (SE = 4.871) # rho = -0.759 ## Not run: fit1 <- riley(Daniels) #maxit reached, try again with more iterations ## End(Not run) fit1 <- riley(Daniels, control=list(maxit=10000)) summary(fit1) # Meta-analysis of prognostic test studies fit2 <- riley(Kertai) fit2 # Meta-analysis of computed tomography data ds <- Scheidler[which(Scheidler$modality==1),] fit3 <- riley(ds) fit3data(Scheidler) data(Daniels) data(Kertai) # Meta-analysis of potential surrogate markers data # The results obtained by Riley (2008) were as follows: # beta1 = -0.042 (SE = 0.063), # beta2 = 14.072 (SE = 4.871) # rho = -0.759 ## Not run: fit1 <- riley(Daniels) #maxit reached, try again with more iterations ## End(Not run) fit1 <- riley(Daniels, control=list(maxit=10000)) summary(fit1) # Meta-analysis of prognostic test studies fit2 <- riley(Kertai) fit2 # Meta-analysis of computed tomography data ds <- Scheidler[which(Scheidler$modality==1),] fit3 <- riley(ds) fit3
Function to display running means of a fitted Bayesian meta-analysis model.
rmplot(...)rmplot(...)
... |
Additional arguments, which are currently ignored. |
This is a generic function.
A ggplot object.
Thomas Debray <[email protected]>
Function to display running means of a fitted Bayesian meta-analysis model.
## S3 method for class 'mcmc.list' rmplot(x, P, greek = FALSE, ...)## S3 method for class 'mcmc.list' rmplot(x, P, greek = FALSE, ...)
x |
An object of class |
P |
Optional dataframe describing the parameters to plot and their respective names |
greek |
Logical value indicating whether parameter labels have to be parsed to get Greek letters. Defaults to FALSE. |
... |
Additional arguments which are passed to ggs_running |
A ggplot object.
Thomas Debray <[email protected]>
Function to display running means of a fitted Bayesian meta-analysis model.
## S3 method for class 'uvmeta' rmplot(x, ...)## S3 method for class 'uvmeta' rmplot(x, ...)
x |
An object of class |
... |
Additional arguments which are currently not used |
Results are displayed for the estimated mean (mu) and standard-deviation (tau) of the meta-analysis model.
A ggplot object.
An object of class ggplot
Thomas Debray <[email protected]>
## Not run: data(Roberts) fit <- with(Roberts, uvmeta(r=SDM, r.se=SE, labels=rownames(Roberts), method="BAYES")) rmplot(fit) ## End(Not run)## Not run: data(Roberts) fit <- with(Roberts, uvmeta(r=SDM, r.se=SE, labels=rownames(Roberts), method="BAYES")) rmplot(fit) ## End(Not run)
Function to display running means of a fitted Bayesian meta-analysis model.
## S3 method for class 'valmeta' rmplot(x, ...)## S3 method for class 'valmeta' rmplot(x, ...)
x |
An object of class |
... |
Additional arguments which are currently not used |
Results are displayed for the estimated mean (mu) and standard-deviation (tau) of the meta-analysis model.
A ggplot object.
An object of class ggplot
Thomas Debray <[email protected]>
## Not run: data(EuroSCORE) fit <- valmeta(cstat=c.index, cstat.se=se.c.index, cstat.cilb=c.index.95CIl, cstat.ciub=c.index.95CIu, N=n, O=n.events, data=EuroSCORE, method="BAYES", slab=Study) rmplot(fit) ## End(Not run)## Not run: data(EuroSCORE) fit <- valmeta(cstat=c.index, cstat.se=se.c.index, cstat.cilb=c.index.95CIl, cstat.ciub=c.index.95CIu, N=n, O=n.events, data=EuroSCORE, method="BAYES", slab=Study) rmplot(fit) ## End(Not run)
Data frame with summary data from 14 comparative studies.
data("Roberts")data("Roberts")
One data frame with 2 variables.
Effect sizes (standardized differences in means)
Standard error of the effect sizes
The Roberts data set is a meta-analysis of 14 studies comparing 'set shifting' ability (the ability to move back and forth between different tasks) in people with eating disorders and healthy controls.
Roberts ME, Tchanturia K, Stahl D, Southgate L, Treasure J. A systematic review and meta-analysis of set-shifting ability in eating disorders. Psychological Medicine 2007, 37: 1075–1084.
Higgins JPT, Thompson SG, Spiegelhalter DJ. A re-evaluation of random-effects meta-analysis. Journal of the Royal Statistical Society. Series A (Statistics in Society) 2009, 172: 137–159.
Data frame with diagnostic accuracy data from three imaging techniques for the diagnosis of lymph node metastasis in women with cervical cancer.
data("Scheidler")data("Scheidler")
One data frame with 6 variables.
string . author of article
integer . type of test (1=CT, 2=LAG, 3=MRI)
integer. number of true positives
integer. number of false negatives
integer. number of false positives
integer. number of true negatives
The Scheidler data comprises the results from a meta-analysis where three imaging techniques for the diagnosis of lymph node metastasis in women with cervical cancer are compared. Forty-four studies in total were included: 17 studies evaluated lymphangiography, another 17 studies examined computed tomography and the remaining 10 studies focused on magnetic resonance imaging. Diagnosis of metastatic disease by lymphangiography (LAG) is based on the presence of nodal-filling defects, whereas computed tomography (CT) and magnetic resonance imaging (MRI) rely on nodal enlargement.
Scheidler J, Hricak H, Yu KK, Subak L, Segal MR. Radiological evaluation of lymph node metastases in patients with cervical cancer. A meta-analysis. Journal of the American Medical Association 1997; 278: 1096–1101.
Reitsma J, Glas A, Rutjes A, Scholten R, Bossuyt P, Zwinderman A. Bivariate analysis of sensitivity and specificity produces informative summary measures in diagnostic reviews. Journal of Clinical Epidemiology 2005; 58: 982–990.
Obtain standard errors or variances of a model fit
se(object, ...) variances(object, ...) tau(object, ...) tau2(object, ...)se(object, ...) variances(object, ...) tau(object, ...) tau2(object, ...)
object |
A model fit object |
... |
other arguments |
For se the standard errors of object, and for
variances the variances. For tau the heterogeneity of the coefficients.
Valentijn de Jong
This function combines one or more existing prediction models into a so/called meta-model.
stackedglm(models, family = binomial, data)stackedglm(models, family = binomial, data)
models |
a list containing the historical prediction models, which can be defined in several ways. For instance,
historical regression models can be specified using a named vector containing the regression coefficients of the
individual predictors (no need to include the intercept term). List items may also represent an object for
which the function |
family |
a description of the error distribution and link function to be used in the meta-model. This can be a character string naming a family function, a family function or the result of a call to a family function. (See family for details of family functions.) |
data |
an optional data frame, list or environment (or object coercible by as.data.frame to a data frame)
containing the variables in the model. If not found in |
Thomas Debray <[email protected]>
Return a model from the cross-validation procedure or the final 'global' model. Caution: This function is still under development.
## S3 method for class 'metapred' subset( x, select = "cv", step = NULL, model = NULL, stratum = NULL, add = TRUE, ... )## S3 method for class 'metapred' subset( x, select = "cv", step = NULL, model = NULL, stratum = NULL, add = TRUE, ... )
x |
metapred object |
select |
Which type of model to select: "cv" (default), "global", or (experimental) "stratified", or "stratum". |
step |
Which step should be selected? Defaults to the best step. numeric is converted to name of the step: 0 for an unchanged model, 1 for the first change... |
model |
Which model change should be selected? NULL (default, best change) or character name of variable or (integer) index of model change. |
stratum |
Experimental. Stratum to return if select = "stratum". |
add |
Logical. Add data, options and functions to the resulting object? Defaults to |
... |
For compatibility only. |
An object of class mp.cv for select = "cv" and an object of class mp.global for select = "global".
In both cases, additional data is added to the resulting object, thereby making it suitable for further methods.
Valentijn de Jong
data(DVTipd) DVTipd$cluster <- letters[1:4] # Add a fictional clustering to the data. mp <- metapred(DVTipd, strata = "cluster", formula = dvt ~ histdvt + ddimdich, family = binomial) subset(mp) # best cross-validated model subset(mp, select = "global") # Final model fitted on all strata. subset(mp, step = 1) # The best model of step 1 subset(mp, step = 1, model = "histdvt") # The model in which histdvt was removed, in step 1.data(DVTipd) DVTipd$cluster <- letters[1:4] # Add a fictional clustering to the data. mp <- metapred(DVTipd, strata = "cluster", formula = dvt ~ histdvt + ddimdich, family = binomial) subset(mp) # best cross-validated model subset(mp, select = "global") # Final model fitted on all strata. subset(mp, step = 1) # The best model of step 1 subset(mp, step = 1, model = "histdvt") # The model in which histdvt was removed, in step 1.
Parameter summaries Provides the summary estimates of the alternative model for bivariate random-effects meta-analysis by Riley et al. (2008) with their corresponding standard errors (derived from the inverse Hessian). For confidence intervals, asymptotic normality is assumed.
## S3 method for class 'riley' summary(object, ...)## S3 method for class 'riley' summary(object, ...)
object |
A |
... |
Arguments to be passed on to other functions (currently ignored) |
For meta-analysis of diagnostic test accuracy data, beta1 equals the logit sensitivity (Sens) and
beta2 equals the logit false positive rate (FPR).
array with confidence intervals for the estimated model parameters. For diagnostic test accuracy data, the resulting summary sensitivity and false positive rate are included.
For the overall correlation (rho) confidence intervals are derived using the transformation
logit((rho+1)/2). Similarly, the logit transformation is used to derive confidence intervals for the summary
sensitivity and false positive rate.
Thomas Debray <[email protected]>
Riley RD, Thompson JR, Abrams KR. An alternative model for bivariate random-effects meta-analysis when the within-study correlations are unknown. Biostatistics 2008; 9: 172–186.
This function provides summary estimates of a fitted univariate meta-analysis model.
## S3 method for class 'uvmeta' summary(object, ...)## S3 method for class 'uvmeta' summary(object, ...)
object |
An object of class |
... |
Optional arguments to be passed on to other functions |
Thomas Debray <[email protected]>
Borenstein M, Hedges LV, Higgins JPT, Rothstein HR. A basic introduction to fixed-effect and random-effects models for meta-analysis. Research Synthesis Methods 2010; 1: 97–111.
DerSimonian R, Laird N. Meta-analysis in clinical trials. Controlled Clinical Trials 1986; 7: 177–188.
Riley RD, Higgins JPT, Deeks JJ. Interpretation of random effects meta-analyses. British Medical Journal 2011; 342: d549.
Tzoulaki et al. (2009) reviewed studies that evaluated various candidate prognostic factors in their ability to improve prediction of coronary heart disease (CHD) or other outcomes beyond what the Framingham risk score (FRS) can achieve.
data("Tzoulaki")data("Tzoulaki")
A data frame containing data from 27 studies on the following 2 variables.
PubmedIDa character vector with the Pubmed ID of the study
Na numeric vector describing the study size
N.eventsa numeric vector describing the observed number of events
FRS.orig.refitteda boolean vector describing whether the coefficients of the original Framingham Risk Score (FRS) were re-estimated
FRS.modif.refitteda boolean vector describing whether the coefficients of the modified Framingham Risk Score were re-estimated
predictorsa character vector indicating which new risk factor(s) were included in the modified FRS
outcomea character vector indicating the primary outcome being predicted
AUC.origa numeric vector describing the Area under the ROC curve (AUC) of the original FRS model
AUC.orig.CIla numeric vector describing the lower boundary of the 95% confidence interval of the AUC of the original FRS model
AUC.orig.CIua numeric vector describing the upper boundary of the 95% confidence interval of the AUC of the original FRS model
AUC.modifa numeric vector describing the Area under the ROC curve (AUC) of the modified FRS model that includes one or more new risk factors
AUC.modif.CIla numeric vector describing the lower boundary of the 95% confidence interval of the AUC of the modified FRS model
AUC.modif.CIua numeric vector describing the upper boundary of the 95% confidence interval of the AUC of the modified FRS model
pval.AUCdiffa numeric vector with the p-value of the difference between AUC.orig and AUC.modif
sign.AUCdiffa boolean vector indicating whether the difference between AUC.orig and AUC.modif is below 0.05
Tzoulaki I, Liberopoulos G, Ioannidis JPA. Assessment of claims of improved prediction beyond the Framingham risk score. JAMA. 2009 Dec 2;302(21):2345–52.
data(Tzoulaki) ## maybe str(Tzoulaki) ; plot(Tzoulaki) ...data(Tzoulaki) ## maybe str(Tzoulaki) ; plot(Tzoulaki) ...
This function summarizes multiple estimates for a single parameter by assuming a fixed (i.e. common) effect or random effects across studies. The summary estimate is obtained by calculating a weighted mean that accounts for sample size and (in case random effects are assumed) for between-study heterogeneity.
uvmeta( r, r.se, r.vi, method = "REML", test = "knha", labels, na.action, n.chains = 4, pars, verbose = FALSE, ... )uvmeta( r, r.se, r.vi, method = "REML", test = "knha", labels, na.action, n.chains = 4, pars, verbose = FALSE, ... )
r |
Vector of numerics containing the effect size of each study |
r.se |
Vector of numerics containing the standard error of the effect sizes |
r.vi |
Vector of numerics containing the sampling variance of the effect sizes |
method |
Character string specifying whether a fixed-effect or a random-effects model should be fitted.
A fixed-effect model is fitted when using |
test |
Optional character string when |
labels |
Optional vector of characters containing the labels for the studies |
na.action |
A function which indicates what should happen when the data contain NAs.
Defaults to |
n.chains |
Optional numeric specifying the number of chains to use in the Gibbs sampler ( |
pars |
Optional list with additional arguments. The width of confidence, credibility and prediction intervals is
defined by |
verbose |
If TRUE then messages generated during the fitting process will be displayed. |
... |
Additional arguments that are passed to rma or runjags (if |
Unless specified otherwise, all meta-analysis models assume random effects and are fitted using restricted
maximum likelihood estimation with the metafor package (Viechtbauer 2010). Further, confidence intervals for
the average performance are based on the Hartung-Knapp-Sidik-Jonkman method, to better account for the uncertainty
in the estimated between-study heterogeneity (Debray 2016). A Bayesian meta-analysis can be performed by specifying
method="BAYES". In that case, the R packages runjags and rjags must be installed.]
For random-effects models, a prediction interval for the pooled effect size is displayed. This interval predicts in what
range future effect sizes will fall given what has already been observed (Higgins 2009, Riley 2011).
For Bayesian meta-analysis models that involve the Gibbs sampler (method="BAYES"), the R packages runjags
and rjags must be installed. The Bayesian approach uses an uninformative Normal prior for the mean and a
uniform prior for the between-study variance of the pooled effect size (Higgins 2009). By default, the Normal prior
has a mean of 0 and a variance of 1000. These hyperparameters can, however, be altered through the
variables hp.mu.mean and hp.mu.var in the argument pars. The prior distribution of the between-study
standard deviation is given by a uniform distribution, by default bounded between 0 and 100.
An object of the class uvmeta for which many standard methods are available.
array with (transformed) data used for meta-analysis, and method(s) used for restoring missing information.
character string specifying the meta-analysis method.
estimated performance statistic of the model. For Bayesian meta-analysis, the posterior median is returned.
standard error (or posterior standard deviation) of the summary estimate.
estimated amount of (residual) heterogeneity. Always 0 when method="FE". For Bayesian meta-analysis, the posterior median is returned.
estimated standard error (or posterior standard deviation) of the between-study variation.
lower bound of the confidence (or credibility) interval of the summary estimate
upper bound of the confidence (or credibility) interval of the summary estimate
lower bound of the (approximate) prediction interval of the summary estimate
upper bound of the (approximate) prediction interval of the summary estimate
the full results from the fitted model
vector specifying the label of each study.
Thomas Debray <[email protected]>
Biggerstaff BJ, Tweedie RL. Incorporating variability in estimates of heterogeneity in the random effects model in meta-analysis. Statistics in Medicine 1997; 16: 753–768.
Borenstein M, Hedges LV, Higgins JPT, Rothstein HR. A basic introduction to fixed-effect and random-effects models for meta-analysis. Research Synthesis Methods 2010; 1: 97–111. doi:10.1002/jrsm.12
DerSimonian R, Laird N. Meta-analysis in clinical trials. Controlled Clinical Trials 1986; 7: 177–188.
Graham PL, Moran JL. Robust meta-analytic conclusions mandate the provision of prediction intervals in meta-analysis summaries. Journal of Clinical Epidemiology 2012; 65: 503–510.
Higgins JPT, Thompson SG. Quantifying heterogeneity in a meta-analysis. Statistics in Medicine 2002; 21: 1539–1558.
Higgins JPT, Thompson SG, Spiegelhalter DJ. A re-evaluation of random-effects meta-analysis. J R Stat Soc Ser A Stat Soc. 2009;172:137–59. doi:10.1111/j.1467-985X.2008.00552.x
Riley RD, Higgins JPT, Deeks JJ. Interpretation of random effects meta-analyses. British Medical Journal 2011; 342: d549. doi:10.1136/bmj.d549
Viechtbauer W. Conducting Meta-Analyses in R with the metafor Package. Journal of Statistical Software. 2010; 36. doi:10.18637/jss.v036.i03
data(Roberts) # Frequentist random-effects meta-analysis fit1 <- with(Roberts, uvmeta(r=SDM, r.se=SE, labels=rownames(Roberts))) summary(fit1) plot(fit1) #show a forest plot fit1 ## Not run: # Bayesian random effects meta-analysis fit2 <- with(Roberts, uvmeta(r=SDM, r.se=SE, labels=rownames(Roberts), method="BAYES")) plot(fit2) ## End(Not run)data(Roberts) # Frequentist random-effects meta-analysis fit1 <- with(Roberts, uvmeta(r=SDM, r.se=SE, labels=rownames(Roberts))) summary(fit1) plot(fit1) #show a forest plot fit1 ## Not run: # Bayesian random effects meta-analysis fit2 <- with(Roberts, uvmeta(r=SDM, r.se=SE, labels=rownames(Roberts), method="BAYES")) plot(fit2) ## End(Not run)
This class encapsulates results of a univariate meta-analysis.
Objects can be created by calls of the form uvmeta.
call:(language) The call to uvmeta.
data:(data frame) The data used for the meta-analysis.
results:(data frame) Contains the pooled effect size (mu), the between-study variability (tausq), Cochran's Q statistic (Q) and Higgins' and Thompson's I square statistic (Isq). For each estimate, error variances are provided with predefined confidence (method="MOM") or credibility (method="bayes") intervals.
model:(character) The meta-analysis model used.
method:(character) The estimator used.
na.action:(character) Information from the action which was applied to object if NAs were handled specially, or NULL.
df:(numeric) Degrees of freedom.
numstudies:(numeric) The amount of studies used in the meta-analysis.
pred.int:(data frame) A prediction interval, predicting in what range future effect sizes will fall given what has already been observed (based on a Student's t-distribution, cfr. Riley 2011)
formula:(character) If a formula was specified, a character vector giving the formula and parameter specifications.
signature(object = "uvmeta"): Print object summary.
signature(object = "uvmeta"): Plot a forest plot with the summary estimate.
signature(object = "uvmeta"): Generate object summary.
data(Collins) #Extract effect size and error variance r <- Collins$logOR vars <- Collins$SE**2 #Frequentist random-effects meta-analysis fit1 <- uvmeta(r,vars) #Extract results fit1$resultsdata(Collins) #Extract effect size and error variance r <- Collins$logOR vars <- Collins$SE**2 #Frequentist random-effects meta-analysis fit1 <- uvmeta(r,vars) #Extract results fit1$results
This function provides summary estimates for the concordance statistic, the total observed-expected ratio or the calibration slope. Where appropriate, data transformations are applied and missing information is derived from available quantities. Unless specified otherwise, all meta-analysis models assume random effects and are fitted using restricted maximum likelihood estimation with the metafor package (Viechtbauer 2010). Further, confidence intervals for the average performance are based on the Hartung-Knapp-Sidik-Jonkman method. When conducting a Bayesian meta-analysis, the R packages runjags and rjags must be installed.
valmeta( measure = "cstat", cstat, cstat.se, cstat.cilb, cstat.ciub, cstat.cilv, sd.LP, OE, OE.se, OE.cilb, OE.ciub, OE.cilv, citl, citl.se, N, O, E, Po, Po.se, Pe, data, method = "REML", test = "knha", verbose = FALSE, slab, n.chains = 4, pars, ... )valmeta( measure = "cstat", cstat, cstat.se, cstat.cilb, cstat.ciub, cstat.cilv, sd.LP, OE, OE.se, OE.cilb, OE.ciub, OE.cilv, citl, citl.se, N, O, E, Po, Po.se, Pe, data, method = "REML", test = "knha", verbose = FALSE, slab, n.chains = 4, pars, ... )
measure |
A character string indicating which summary performance
measure should be calculated. Options are
|
cstat |
Optional vector with the estimated c-statistic for each valiation |
cstat.se |
Optional vector with the standard error of the estimated c-statistics |
cstat.cilb |
Optional vector to specify the lower limits of the confidence interval. |
cstat.ciub |
Optional vector to specify the upper limits of the confidence interval. |
cstat.cilv |
Optional vector to specify the levels of aformentioned confidence interval limits. (default: 0.95, which corresponds to the 95% confidence interval). |
sd.LP |
Optional vector with the standard deviation of the linear predictor (prognostic index) |
OE |
Optional vector with the estimated ratio of total observed versus total expected events |
OE.se |
Optional vector with the standard errors of the estimated O:E ratios |
OE.cilb |
Optional vector to specify the lower limits of the confidence interval for |
OE.ciub |
Optional vector to specify the upper limits of the confidence interval for |
OE.cilv |
Optional vector to specify the levels of aformentioned confidence interval limits. (default: 0.95, which corresponds to the 95% confidence interval). |
citl |
Optional vector with the estimated calibration-in-the-large for each valiation |
citl.se |
Optional vector with the standard error of the estimated calibration-in-the-large statistics |
N |
Optional vector with the total number of participants for each valiation |
O |
Optional vector with the total number of observed events for each valiation
(if specified, during time |
E |
Optional vector with the total number of expected events for each valiation
(if specified, during time |
Po |
Optional vector with the (cumulative) observed event probability for each valiation
(if specified, during time |
Po.se |
Optional vector with the standard errors of |
Pe |
Optional vector with the (cumulative) expected event probability for each validation
(if specified, during time |
data |
optional data frame containing the variables given to the arguments above. |
method |
Character string specifying whether a fixed- or a random-effects model should be fitted.
A fixed-effects model is fitted when using |
test |
Optional character string specifying how test statistics and confidence intervals for the fixed effects
should be computed. By default ( |
verbose |
If TRUE then messages generated during the fitting process will be displayed. |
slab |
Optional vector specifying the label for each study |
n.chains |
Optional numeric specifying the number of chains to use in the Gibbs sampler
(if |
pars |
A list with additional arguments. See 'Details' for more information. The following parameters configure the MCMC sampling procedure:
|
... |
Additional arguments that are passed to rma or runjags (if |
A summary estimate for the concorcance (c-) statistic can be obtained by specifying measure="cstat".
The c-statistic is a measure of discrimination, and indicates the ability of a prediction model to
distinguish between patients developing and not developing the outcome. The c-statistic typically ranges
from 0.5 (no discriminative ability) to 1 (perfect discriminative ability). When missing, the c-statistic
and/or its standard error are derived from other reported information.
See ccalc for more information.
By default, it is assumed that the logit of the c-statistic is Normally distributed within and across studies
(pars$model.cstat = "normal/logit"). Alternatively, it is possible to assume that the raw c-statistic
is Normally distributed across studies pars$model.cstat = "normal/identity".
A summary estimate for the total observed versus expected (O:E) ratio can be
obtained by specifying measure="OE". The total O:E ratio provides a
rough indication of the overall model calibration (across the entire range
of predicted risks). When missing, the total O:E ratio and/or its standard
error are derived from other reported information. See oecalc
for more information.
For frequentist meta-analysis, within-study variation can either be modeled
using a Normal (model.oe = "normal/log" or
model.oe = "normal/identity") or a Poisson distribution
(model.oe = "poisson/log").
When performing a Bayesian meta-analysis, all data are modeled using a
one-stage random effects (hierarchical related regression) model.
In particular, a binomial distribution (if O, E and
N is known), a Poisson distribution (if only O and
E are known) or a Normal distribution (if OE and
OE.se or OE.95CI are known) is selected separately for
each study.
All Bayesian meta-analysis models assume the presence of random effects. Summary estimates are based on the posterior mean. Credibility and prediction intervals are directly obtained from the corresponding posterior quantiles.
The prior distribution for the (transformed) performance estimate is modeled
using a Normal distribution, with mean hp.mu.mean (defaults to 0)
and variance hp.mu.var (defaults to 1000).
For meta-analysis of the total O:E ratio, the maximum value for
hp.mu.var is 100.
By default, the prior distribution for the between-study standard deviation
is modeled using a uniform distribution (hp.tau.dist="dunif"),
with boundaries hp.tau.min and hp.tau.max. Alternative choices
are a truncated Student-t distribution (hp.tau.dist="dhalft") with a
mean of hp.tau.mean, a standard deviation of hp.tau.sigma and
hp.tau.df degrees of freedom. This distribution is again restricted
to the range hp.tau.min to hp.tau.max.
An object of class valmeta with the following elements:
array with (transformed) data used for meta-analysis, and method(s) used for restoring missing information.
character string specifying the performance measure that has been meta-analysed.
character string specifying the meta-analysis method.
character string specifying the meta-analysis model (link function).
summary estimate for the performance statistic. For Bayesian meta-analysis, the posterior median is returned.
lower bound of the confidence (or credibility) interval of the summary performance estimate.
upper bound of the confidence (or credibility) interval of the summary performance estimate.
lower bound of the (approximate) prediction interval of the summary performance estimate.
upper bound of the (approximate) prediction interval of the summary performance estimate.
the full results from the fitted model.
vector specifying the label of each study.
The width of calculated confidence, credibility and prediction
intervals can be specified using level in the pars argument
(defaults to 0.95).
Debray TPA, Damen JAAG, Snell KIE, Ensor J, Hooft L, Reitsma JB, et al. A guide to systematic review and meta-analysis of prediction model performance. BMJ. 2017;356:i6460. doi:10.1136/bmj.i6460
Debray TPA, Damen JAAG, Riley R, Snell KIE, Reitsma JB, Hooft L, et al. A framework for meta-analysis of prediction model studies with binary and time-to-event outcomes. Stat Methods Med Res. 2019;28:2768–86. doi:10.1177/0962280218785504
Riley RD, Tierney JF, Stewart LA. Individual participant data meta-analysis: a handbook for healthcare research. Hoboken, NJ: Wiley; 2021. ISBN: 978-1-119-33372-2.
Steyerberg EW, Nieboer D, Debray TPA, van Houwelingen HC. Assessment of heterogeneity in an individual participant data meta-analysis of prediction models: An overview and illustration. Stat Med. 2019; 38:4290–309. doi:10.1002/sim.8296
Viechtbauer W. Conducting Meta-Analyses in R with the metafor Package. Journal of Statistical Software. 2010; 36. doi:10.18637/jss.v036.i03
ccalc to calculate concordance statistics
and corresponding standard errors, oecalc to calculate the
total O:E ratio and corresponding standard errors,
plot.valmeta to generate forest plots
######### Validation of prediction models with a binary outcome ######### data(EuroSCORE) # Meta-analysis of the c-statistic (random effects) fit <- valmeta(cstat=c.index, cstat.se=se.c.index, cstat.cilb=c.index.95CIl, cstat.ciub=c.index.95CIu, cstat.cilv=0.95, N=n, O=n.events, slab=Study, data=EuroSCORE) plot(fit) # Nearly identical results when we need to estimate the SE valmeta(cstat=c.index, N=n, O=n.events, slab=Study, data=EuroSCORE) # Two-stage meta-analysis of the total O:E ratio (random effects) valmeta(measure="OE", O=n.events, E=e.events, N=n, slab=Study, data=EuroSCORE) valmeta(measure="OE", O=n.events, E=e.events, data=EuroSCORE) valmeta(measure="OE", Po=Po, Pe=Pe, N=n, data=EuroSCORE) ## Not run: # One-stage meta-analysis of the total O:E ratio (random effects) valmeta(measure="OE", O=n.events, E=e.events, data=EuroSCORE, method="ML", pars=list(model.oe="poisson/log")) # Bayesian random effects meta-analysis of the c-statistic fit2 <- valmeta(cstat=c.index, cstat.se=se.c.index, cstat.cilb=c.index.95CIl, cstat.ciub=c.index.95CIu, cstat.cilv=0.95, N=n, O=n.events, data=EuroSCORE, method="BAYES", slab=Study) # Bayesian one-stage random effects meta-analysis of the total O:E ratio # Consider that some (but not all) studies do not provide information on N # A Poisson distribution will be used for studies 1, 2, 5, 10 and 20 # A Binomial distribution will be used for the remaining studies EuroSCORE.new <- EuroSCORE EuroSCORE.new$n[c(1, 2, 5, 10, 20)] <- NA pars <- list(hp.tau.dist="dhalft", # Prior for the between-study standard deviation hp.tau.sigma=1.5, # Standard deviation for 'hp.tau.dist' hp.tau.df=3, # Degrees of freedom for 'hp.tau.dist' hp.tau.max=10, # Maximum value for the between-study standard deviation seed=5) # Set random seed for the simulations fit3 <- valmeta(measure="OE", O=n.events, E=e.events, N=n, data=EuroSCORE.new, method="BAYES", slab=Study, pars=pars) plot(fit3) print(fit3$fit$model) # Inspect the JAGS model print(fit3$fit$data) # Inspect the JAGS data ## End(Not run) ######### Validation of prediction models with a time-to-event outcome ######### data(Framingham) # Meta-analysis of total O:E ratio after 10 years of follow-up valmeta(measure="OE", Po=Po, Pe=Pe, N=n, data=Framingham)######### Validation of prediction models with a binary outcome ######### data(EuroSCORE) # Meta-analysis of the c-statistic (random effects) fit <- valmeta(cstat=c.index, cstat.se=se.c.index, cstat.cilb=c.index.95CIl, cstat.ciub=c.index.95CIu, cstat.cilv=0.95, N=n, O=n.events, slab=Study, data=EuroSCORE) plot(fit) # Nearly identical results when we need to estimate the SE valmeta(cstat=c.index, N=n, O=n.events, slab=Study, data=EuroSCORE) # Two-stage meta-analysis of the total O:E ratio (random effects) valmeta(measure="OE", O=n.events, E=e.events, N=n, slab=Study, data=EuroSCORE) valmeta(measure="OE", O=n.events, E=e.events, data=EuroSCORE) valmeta(measure="OE", Po=Po, Pe=Pe, N=n, data=EuroSCORE) ## Not run: # One-stage meta-analysis of the total O:E ratio (random effects) valmeta(measure="OE", O=n.events, E=e.events, data=EuroSCORE, method="ML", pars=list(model.oe="poisson/log")) # Bayesian random effects meta-analysis of the c-statistic fit2 <- valmeta(cstat=c.index, cstat.se=se.c.index, cstat.cilb=c.index.95CIl, cstat.ciub=c.index.95CIu, cstat.cilv=0.95, N=n, O=n.events, data=EuroSCORE, method="BAYES", slab=Study) # Bayesian one-stage random effects meta-analysis of the total O:E ratio # Consider that some (but not all) studies do not provide information on N # A Poisson distribution will be used for studies 1, 2, 5, 10 and 20 # A Binomial distribution will be used for the remaining studies EuroSCORE.new <- EuroSCORE EuroSCORE.new$n[c(1, 2, 5, 10, 20)] <- NA pars <- list(hp.tau.dist="dhalft", # Prior for the between-study standard deviation hp.tau.sigma=1.5, # Standard deviation for 'hp.tau.dist' hp.tau.df=3, # Degrees of freedom for 'hp.tau.dist' hp.tau.max=10, # Maximum value for the between-study standard deviation seed=5) # Set random seed for the simulations fit3 <- valmeta(measure="OE", O=n.events, E=e.events, N=n, data=EuroSCORE.new, method="BAYES", slab=Study, pars=pars) plot(fit3) print(fit3$fit$model) # Inspect the JAGS model print(fit3$fit$data) # Inspect the JAGS data ## End(Not run) ######### Validation of prediction models with a time-to-event outcome ######### data(Framingham) # Meta-analysis of total O:E ratio after 10 years of follow-up valmeta(measure="OE", Po=Po, Pe=Pe, N=n, data=Framingham)
Returns the variance-covariance matrix of the main parameters of a fitted model object.
## S3 method for class 'riley' vcov(object, ...)## S3 method for class 'riley' vcov(object, ...)
object |
a |
... |
arguments to be passed on to other functions |
The variance-covariance matrix is obtained from the inverse Hessian as provided by optim.
A matrix of the estimated covariances between the parameter estimates in the Riley model: logit of sensitivity (mu1), logit of false positive rate (mu2), additional variation of mu1 beyond sampling error (psi1), additional variation of mu2 beyond sampling error (psi2) and a transformation of the correlation between psi1 and psi2 (rhoT). The original correlation is given as inv.logit(rhoT)*2-1.
A warning message is casted when the Hessian matrix contains negative eigenvalues. This implies that the identified minimum for the (restricted) negative log-likelihood is a saddle point, and that the solution is therefore not optimal.
Thomas Debray <[email protected]>
Riley, RD., Thompson, JR., & Abrams, KR. (2008). “An alternative model for bivariate random-effects meta-analysis when the within-study correlations are unknown.” Biostatistics, 9, 172–186.
This dataset comprises the results from 16 studies assessing the prognostic role of human epidermal growth factor receptor 2 (HER2) in endometrial cancer. These studies were previously identified in a systematic review by Zhang et al. to evaluate the overall risk of several hormone receptors for endometrial cancer survival.
data("Zhang")data("Zhang")
A data frame with 20 observations on the following 10 variables.
Studya factor with 16 levels to indicate the study
PrimaryAuthora factor indicating the first author's last name
yeara numeric vector indicating the publication year
Countrya factor indicating the source country of the study data
Diseasea factor indicating the studied disease. Possible levels are EC (endometrial cancer), EEC (endometrioid endometrial cancer) and UPSC (uterine papillary serous carcinoma)
Na numeric vector describing the total sample size of each study
HRa numeric vector describing the estimated hazard ratio of each study
HR.025a numeric vector describing the lower boundary of the 95% confidence interval of HR
HR.975a numeric vector describing the upper boundary of the 95% confidence interval of HR
outcomea factor indicating the studied outcome. Possible levels are OS (overall survival) and PFS (progression-free survival)
Eligible studies were identified by searching the PubMed and EMBASE databases for publications from 1979 to May 2014. Data were collected from studies comparing overall survival or progression-free survival in patients with elevated levels of human epidermal growth factor receptor 2 with those in patients with lower levels.
Zhang Y, Zhao D, Gong C, Zhang F, He J, Zhang W, et al. Prognostic role of hormone receptors in endometrial cancer: a systematic review and meta-analysis. World J Surg Oncol. 2015 Jun 25;13:208.
Riley RD, Jackson D, Salanti G, Burke DL, Price M, Kirkham J, et al. Multivariate and network meta-analysis of multiple outcomes and multiple treatments: rationale, concepts, and examples. BMJ. 2017 13;358:j3932.
data(Zhang) # Display the hazard ratios for overall survival in a forest plot ds <- subset(Zhang, outcome=="OS") with(ds, forest(theta = HR, theta.ci.lb = HR.025, theta.ci.ub = HR.975, theta.slab = Study, xlab = "Hazard ratio of HER2 versus OS", refline = 1))data(Zhang) # Display the hazard ratios for overall survival in a forest plot ds <- subset(Zhang, outcome=="OS") with(ds, forest(theta = HR, theta.ci.lb = HR.025, theta.ci.ub = HR.975, theta.slab = Study, xlab = "Hazard ratio of HER2 versus OS", refline = 1))