knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) library(Rdpack)
trt.events <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 5, 0, 2, 2, 26) trt.total <- c(121, 59, 358, 454, 238, 62, 186, 87, 369, 257, 798, 978, 605, 277, 157, 650) ctrl.events <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 17, 35, 2, 5, 1, 12) ctrl.total <- c(118, 58, 164, 216, 268, 59, 97, 94, 386, 109, 804, 996, 608, 198, 50, 220) mid.p <- TRUE alpha <- 0.05 arguments <- list(as.name("rema"), "rema.obj" = as.name("antibiotics.rema")) TE <- 0.5708234 CI <-c(0.1947373, 1.5984156) CMLE <- "CMLE" pval <- 0.2555818 test.stat <- c(27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53) norm.probs <- c(6.716051e-12, 6.243005e-10, 4.605874e-08, 1.178779e-06, 1.810190e-05, 1.900967e-04, 1.361220e-03, 6.418245e-03, 2.290489e-02, 5.697199e-02, 1.160254e-01, 1.694432e-01, 2.084688e-01, 1.724891e-01, 1.360043e-01, 6.433921e-02, 3.335175e-02, 8.452742e-03, 3.156814e-03, 3.155256e-04, 8.395551e-05, 2.683688e-06, 7.058532e-07, 6.078317e-09, 1.895755e-09, 3.216023e-12, 5.180658e-14) dist <- data.frame(test.stat, norm.probs) tstat <- 37 antibiotics.rema <- list("trt.events" = trt.events, "trt.total" = trt.total, "ctrl.events" = ctrl.events, "ctrl.total" = ctrl.total, "mid.p" = mid.p, "alpha" = alpha, "arguments"= arguments, "TE" = TE, "CI" = CI, "method"= CMLE, "pval" = pval, "dist" = dist, "tstat" = tstat) class(antibiotics.rema) <- c("rema", "list")
The rema
(rare event meta-analysis) package performs a permutation-based
meta-analysis for heterogeneous, rare event data. Rare events and sparse data
often cause conventional approaches to fail. While recent work in this area has
provided options for sparse data, problems can still arise when heterogeneity
across the available studies differs based on treatment group. This package
provides a meta-analysis approach that often outperforms traditional methods
when such patterns of rare events and heterogeneity are observed.
It is important to note that rema
is designed for small data sets with rare
events. Computation time and memory usage increase rapidly as the number of
observed events increases. Particularly, rema
should not be used on a machine
with a standard amount of RAM if the total number of observations is relatively
large and the observed event rates are relatively large, or if the total
number of observed events is relatively large.
rema
packagelibrary(rema)
The following table contains information from a data set regarding the use of
antibiotics, which we will use to illustrate this package (data provided by
r insert_citeOnly("Spinks2013;textual", "rema")
,
Analysis 4.1 and used in
r insert_citeOnly("Friedrich2007;textual", "rema")
).
Sixteen studies, conducted from 1950 to 2000, compare the use of antibiotics to
prevent acute rheumatic fever, a sore throat resulting from inadequately treated
strep throat or scarlet fever, compared to a placebo. The event of interest is
the occurrence of acute rheumatic fever.
Study | Antibiotics Events | Antibiotics Total | Placebo Events | Placebo Total :-----:|:------------------:|:-----------------:|:--------------:|:-------------: 1 | 0 | 121 | 0 | 118 2 | 0 | 59 | 0 | 58 3 | 0 | 358 | 0 | 164 4 | 0 | 454 | 0 | 216 5 | 0 | 238 | 0 | 268 6 | 0 | 62 | 0 | 59 7 | 0 | 186 | 0 | 97 8 | 0 | 87 | 0 | 94 9 | 0 | 369 | 0 | 386 10 | 0 | 257 | 2 | 109 11 | 2 | 798 | 17 | 804 12 | 5 | 978 | 35 | 996 13 | 0 | 605 | 2 | 608 14 | 2 | 277 | 5 | 198 15 | 2 | 157 | 1 | 50 16 | 26 | 650 | 12 | 220
We will use the rema
package to analyze this data. We note that the event of
acute rheumatic fever is rare, and the data is sparse.
We also notice the relatively large number of events in some of the cells
(particularly, in study 16 and studies 11 and 12 for the placebo group), which
will increase computation time and memory usage. Accordingly, we ran this data
set on a system with an Intel Xeon processor E5-2689 v4 @ 3.10GHz with 192GB of
RAM, and it took about 1 hour and 20 minutes to run.
The default arguments of the rema
function are
rema(trt.events = NULL, trt.total = NULL, ctrl.events = NULL, ctrl.total = NULL, rema.obj, mid.p = TRUE, distr = TRUE, one.sided.p = FALSE, alpha = 0.05)
where trt.events
, trt.total
, ctrl.events
, and ctrl.total
are numeric
vectors containing, for each independent study, the number of observed events in
the treatment group, the total number of observations in the treatment group,
the number of observed events in the control group, and the total number of
observations in the control group, respectively. The remaining arguments specify
whether the p-values and confidence interval
should be adjusted using the mid-p correction to reduce conservatism (mid.p
),
whether the permutational distribution of the test statistic is returned
(distr
), whether the one-sided p-value is returned (one.sided.p
), and the
level of confidence for a (1 - alpha)% confidence interval for the overall
treatment effect (alpha
).
For the antibiotics data example, we will create the four needed vectors as follows.
anti.events <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 5, 0, 2, 2, 26) anti.total <- c(121, 59, 358, 454, 238, 62, 186, 87, 369, 257, 798, 978, 605, 277, 157, 650) plac.events <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 17, 35, 2, 5, 1, 12) plac.total <- c(118, 58, 164, 216, 268, 59, 97, 94, 386, 109, 804, 996, 608, 198, 50, 220)
It is important to note that no continuity correction should be applied to
studies with zero observed events to data sets passed into rema
.
rema
will return a list with elements that include the combined
odds ratio, its associated (1 - alpha)% confidence interval, and the two-sided p-value.
The following examples illustrate how the returned object and resulting output
differ depending on the values passed to the various arguments listed above.
This example runs rema
with the default values. The output contains the
combined odds ratio, its associated 95% confidence interval, and the mid-p
corrected two-sided p-value. Under the details section, the method used when
computing the odds ratio and its confidence interval is reported. If the observed
test statistic is not at either extreme (the minimum or maximum value) of the
distribution, the conditional maximum likelihood estimate (CMLE) will be used;
otherwise, the median unbiased estimate (MUE) will be employed. Here, the CMLE
was used. Overall, this method indicates no decrease the risk of acute
rheumatic fever from the use of antibiotics.
antibiotics.rema <- rema(anti.events, anti.total, plac.events, plac.total) summary(antibiotics.rema) # Call: # rema(trt.events = anti.events, trt.total = anti.total, ctrl.events = plac.events, # ctrl.total = plac.total) # # OR 95%-CI p-value # 0.5708 [0.1947; 1.5984] 0.2556 # # Details on meta-analytical method: # - Rare event, heterogeneous meta-analysis method # - Two-sided p-value returned (mid.p = TRUE) # - Conditional Maximum Likelihood Estimate (CMLE) used when computing the odds ratio
Some may find it interesting to view the permutational distribution. Since the
distribution is saved by default, it can be called and printed or visualized
with the plot
function. The asterisk in the plot marks the observed value of
the test statistic, and the dark grey bars mark the area that contributes to the
two-sided p-value. The observed value of the test statistics is also saved in
the rema
object.
antibiotics.rema$dist antibiotics.rema$tstat
plot(antibiotics.rema)
For the remaining examples, we will use the antibiotics.rema
object as input instead of
the first four vectors. Running rema
once with the first four vectors allows
us to obtain the permutational distribution. Once we have that distribution, the various
arguments pertaining to inference (e.g. mid.p
, alpha
) can be changed
without the need to recompute the distribution, saving computation time.
This example requests the one-sided p-value be returned. The summary output
does not change, but the returned object now has the element pval.one.sided
that can be accessed. This p-value is the sum of the probabilities from test
statistics less than or equal to the observed test statistic (the direction of
the treatment effect), with the mid-p adjustment.
antibiotics.rema.one.pval <- rema(antibiotics.rema, one.sided.p = TRUE) antibiotics.rema.one.pval$pval.one.sided
This example applies no correction to the confidence bounds and the two-sided p-value. As expected, the 95% confidence interval is wider, and the two-sided p-value is larger when the mid-p correction to reduce conservatism is not applied.
rema(antibiotics.rema, mid.p = FALSE)
This example sets the alpha
argument to 0.1, resulting in a 90% confidence
interval. As expected, increasing alpha decreases the confidence interval width.
rema(antibiotics.rema, alpha = 0.1)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.