Perform a Copas selection model analysis for selection bias in
meta-analysis.
The program takes an object of class meta, which is most
easily created by an analysis using one of the functions
metabin, metacont and metagen in the
package meta, performs a 'Copas selection model analysis' and presents
a graphical and tabular summary of the results. An object of class
copas is created and this can be used to recreate the results
table and graphs subsequently, without re-running the analysis, using
the print, summary and plot function.
An object of class meta, obtained from one of the
functions metabin, metacont and metagen in
the package meta.
gamma0.range
(Advanced users only) A numerical vector of
length two specifying the range of gamma0 values the program will
explore.
The parameter gamma0 is the constant in the probit selection
model for study publication. Thus, the cumulative normal of
gamma0 is approximately the probability that a small study is
published (in non-technical terms gamma0 relates to the probability of
publishing a small study, although its values are not restricted to
the range [0,1]; larger values correspond to higher probabilities of
publishing a small study). Most users will not need to specify a
range for this parameter. When no argument is specified, the program
uses an algorithm to determine a suitable range. This is based on
the range of treatment effect standard errors in the meta-analysis,
and is described in more detail below.
gamma1.range
(Advanced users only) A numerical vector of length
two specifying the range of gamma1 values the program will explore.
The parameter gamma1 is the coefficient of study precision
(1/standard error) in the probit selection model for study
publication (in non-technical terms gamma1 relates to the
rate at which the probability of publishing a study
increases as the standard error of the treatment effect it reports
decreases; larger values correspond to higher probabilities of
publishing a small study). Most users will not need to specify a
range for this parameter. When no argument is specified,
the program uses an algorithm to determine a suitable range.
This is based on the range of treatment effect standard errors
in the meta-analysis, and is described in more detail below.
ngrid
The program fits the Copas selection model over a grid
defined by the range of values of gamma0 and gamma1 specified in the
previous two arguments. This parameter fixes the square-root of the
number of points in the grid.
nlevels
(Advanced users only). Fitting the Copas model over the
grid specified by the previous three arguments results in a
treatment estimate at every point in the grid. These can then be
displayed on a contour plot where contours of treatment effect
(z-axis) are shown by gamma0 (x-axis) and gamma1 (y-axis). This
argument specifies the number of contour lines that will be drawn.
Note
(i) Calculations for the contour plot are performed by the function
copas, so this argument has no effect in the plot
function.
(ii) If a large number of contour lines are desired, then you may
wish to consider increasing the grid size (argument ngrid
above).
Leave this option unspecified if you are using the option
levels below.
levels
A numerical vector of treatment values for which contour
lines will be drawn. In more detail, fitting the Copas model over the
grid specified by the arguments gamma0.range,
gamma1.range and ngrid results in a treatment estimate
at every point in the grid. These are then displayed on a contour plot
where contours of treatment effect (z-axis) are shown by gamma0
(x-axis) and gamma1 (y-axis). This argument is a numerical vector
which specifies the treatment effects for which contour lines will be
drawn.
It is usually not a good idea to set this argument for initial runs,
as one does not know the range of treatment values that the contour
plot will cover, and treatment values which do not correspond to
values in the contour plot (defined by the range of gamma0 and
gamma1) will not be plotted.
Note
(i) Calculations for the contour plot are performed by the function
copas, so this argument has no effect in the plot
function.
(ii) Contours will not be drawn if a large number of contour lines
are desired, then you may wish to consider increasing the grid size
(argument ngrid above).
Leave this option unspecified if you are using the option
nlevels above.
slope
A numeric providing the slope of the line approximately
orthogonal to contours in the contour plot. If the argument
slope is NULL (default) the program seeks to estimate
the slope of the contours in the region of the maximum, which are
usually approximately parallel. Most users will leave the argument
slope unspecified, at least for the first analysis of a data
set, but in certain cases setting it manually can improve the
results.
left
A logical indicating whether the cause of any selection
bias is due to missing studies on the left or right of the funnel
plot: left hand side if left=TRUE, right hand side if
left=FALSE. This information is needed in order to be sure
the test for presence of residual selection bias is calculated
correctly. If not set, the linear regression test for funnel plot
asymmetry (i.e., function metabias(..., meth="linreg")) is
used to determine whether studies are missing on the left or right
hand side. In the majority of cases this will work correctly.
rho.bound
(Advanced users only) A number giving the upper bound
for the correlation parameter rho (see details below). This
must be < 1, and usually > 0.95. The lower bound is calculated as
-(the upper bound).
sign.rsb
The significance level for the test of residual
selection bias (between 0 and 1).
backtransf
A logical indicating whether results should be
back transformed in printouts and plots. If backtransf=TRUE
(default), results for sm="OR" are printed as odds ratios
rather than log odds ratio, for example.
silent
A logical indicating whether information on progress in
fitting the Copas selection model should be printed:
silent=TRUE, do not print information (the default);
silent=FALSE, print information.
warn
A number setting the handling of warning messages. It is
not uncommon for numerical problems to be encountered during
estimation over the grid of (gamma0, gamma1) values. Usually this
does not indicate a serious problem. This option specifies what to
do with warning messages. warn=-1: ignore all warnings;
warn=0 (the default): store warnings till function finishes;
if there are less than 10, print them, otherwise print a message
saying warning messages were generated; warn=1: print
warnings as they occur; warn=2: stop the function when the
first warning is generated. For further details see
help(options).
Details
Conduct a Copas selection model analysis to investigate, and
attempt to correct for, selection/publication bias in a meta-analysis.
The Copas selection model consists of two models, which are fitted
jointly. The first is the usual random effects meta-analysis model,
and the second is a selection model, where study i is selected for
publication if Z>0, where
Z = gamma0 + gamma1/(SE(i)) + delta(i)
The error delta(i) is correlated with the error in the random effects
meta-analysis, with correlation rho. If rho=0, the model corresponds
to the usual random effects meta-analysis. As rho moves from 0 to 1,
studies with larger treatment estimates are more likely to be
selected/published.
The software chooses a grid of gamma0 and gamma1 values, corresponding
to a range of selection/publication probabilities for the study with
the largest treatment effect standard error (often the smallest
study). For each value in this grid, the treatment effect is
estimated using the function optim. This information is used to
produce the contour plot (top right panel of output from
plot.copas).
Contours of constant treatment effect are usually locally
parallel. The software estimates the slope of these contours, and
combines this information with other parameter estimates from the
model to explore (i) how the treatment estimate, and its standard
error, change with increasing selection (bottom left panel,
plot.copas) and (ii) how much selection needs to be accounted
for before any remaining asymmetry in the funnel plot is likely to
have occurred by chance (bottom right panel, plot.copas).
A table of results can be produced by the function
summary.copas. A more detail output is provided by the function
print.copas.
For a fuller description of the model, our implementation and
specifically our approach to estimating the locally parallel
contours, see Carpenter et al. (2009) and Schwarzer et al. (2010).
Value
An object of class copas with corresponding print,
summary, plot function. The object is a
list containing the following components:
TE
Vector of treatment effects plotted in treatment effect
plot
seTE
Vector of standard error of TE
TE.random
Usual random effects estimate of treatment effect
seTE.random
Usual standard error of TE.random
left
Whether selection bias expected on left or right
rho.bound
Bound on rho
gamma0.range
Range of gamma0 (see help on copas
arguments above)
gamma1.range
Range of gamma1 (see help on copas
arguments above)
slope
Slope of line approximately orthogonal to contours in
contour plot
regr
A list containing information on regression lines fitted
to contours in contour plot
ngrid
Square root of grid size
nlevels
Number of contour lines
gamma0
Vector of gamma0 values at which model fitted
(determined by gamma0.range and grid). x-axis values for contour
plot
gamma1
vector of gamma1 values at which model fitted
(determined by gamma1.range and grid). y-axis values for contour
plot
TE.contour
Treatment values (ie z-axis values) used to draw
contour plot.
x.slope
x coordinates for 'orthogonal line' in contour plot
y.slope
y coordinates for 'orthogonal line' in contour plot
TE.slope
Vector of treatment values plotted in treatment effect
plot
seTE.slope
Standard error of TE.slope
rho.slope
Vector of estimated rho values corresponding to
treatment estimates in TE.slope
tau.slope
Vector of estimated heterogeneity values
corresponding to treatment estimates in TE.slope
loglik1
Vector of log-likelihood values corresponding to
treatment estimates in TE.slope
conv1
Numerical vector indicating convergence status for each
treatment estimate in TE.slope - see parameter
convergence in function optim
message1
Character vector - translation of conv1
loglik2
Vector of log-likelihoods from fitting model to evaluate
presence of residual selection bias
conv2
Numerical vector indicating convergence status for models
to evaluate presence of residual selection bias - see parameter
convergence in function optim
message2
Character vector - translation of conv2
publprob
Vector of probabilities of publishing the
smallest study, used in x-axis of bottom two panels in function
plot.copas
pval.rsb
P-values for tests on presence of residual selection
bias, plotted in bottom right panel in plot.copas
sign.rsb
The significance level for the test of residual
selection bias
N.unpubl
Approximate number of studies the model suggests
remain unpublished
Carpenter JR, Schwarzer G, R<c3><83><c2><bc>cker G, K<c3><83><c2><bc>nstler R (2009),
Empirical evaluation showed that the Copas selection model provided
a useful summary in 80% of meta-analyses.
Journal of Clinical Epidemiology, 62, 624–631.
Copas J (1999),
What works?: Selectivity models and meta-analysis.
Journal of the Royal Statistical Society, Series A,
162, 95–109.
Copas J, Shi JQ (2000),
Meta-analysis, funnel plots and sensitivity analysis.
Biostatistics, 1, 247–262.
Copas JB, Shi JQ (2001),
A sensitivity analysis for publication bias in systematic reviews.
Statistical Methods in Medical Research, 10, 251–265.
Schwarzer G, Carpenter J, R<c3><83><c2><bc>cker G (2010),
Empirical evaluation suggests Copas selection model preferable to
trim-and-fill method for selection bias in meta-analysis.
Journal of Clinical Epidemiology, 63, 282–288.
##
## Basic example
##
## Load data
##
data(Fleiss93)
##
## Perform meta-analysis
## (Note event.e indicates events, n.e total in exposed arm;
## event.c indicates events, n.c total in control arm)
##
meta1 <- metabin(event.e, n.e, event.c, n.c, data=Fleiss93, sm="OR")
summary(meta1)
##
## To perform a basic Copas-selection model analysis
##
cop1 <- copas(meta1)
plot(cop1)
summary(cop1)
##
## Interpretation:
##
## a. The initial meta-analysis shows the fixed and random effects pooled
## ORs differ; consistent with asymmetry in the funnel plot and
## possible selection bias. Both fixed effect and random effects model
## show a significant treatment effect in this dataset.
##
## b. Plotting the copas analysis shows
##
## (i) funnel plot: asymmetry indicates possible selection bias.
##
## (ii) contour plot treatment effect declines steadily as selection
## increases (no selection, top right, log OR < -0.12; increasing
## selection as move to left of plot, log OR rises to -0.03.
##
## (iii) Treatment effect plot suggests that even with no selection,
## p-value for treatment effect is larger than 0.05 which is
## different from the result of the usual random effects model
## (see output of summary(cop1). This difference is due to the
## use of different methods to estimate the between-study
## variance: maximum-likelihood in Copas analysis compared to
## method-of-moments in usual random effects model.
## The p-value for treatment effect is increasing with
## increasing selection.
##
## (iv) P-value for residual selection bias plot: this shows that even
## with no selection bias, the p-value for residual selection bias
## is non-significant at the 10% level. As expected, as selection
## increases the p-value for residual selection bias increases too.
## Repeat the same example, setting all the arguments of the copas
## function:
##
cop2 <- copas(meta1,
gamma0.range=c(-0.5,2.1), # range of gamma0 parameter
gamma1.range=c(0, 0.08), # range of gamma1 parameter
ngrid=20, # specify a 20X20 grid (finer than default)
levels=c(-0.13, -0.12, -0.1, -0.09, -0.07, -0.05, -0.03),# specify contour lines
slope=0.2, # specify slope of 'orthogonal' line in contour plot
left=FALSE, # as any selection bias due to missing studies on right
rho.bound=0.998, # constrain rho between [-0.998, 0.998]
silent=FALSE, # update user on progress
warn=-1 # suppress warning messages
)
plot(cop2)
##
## Print table of results used to draw treatment effect plot:
##
summary(cop2)
providing 
.