---
title: "rema"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{rema}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
library(Rdpack)
```
```{r, include = FALSE}
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")
```
## Introduction
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.
## The `rema` package
```{r load package, eval=TRUE}
library(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
```{r Data introduction/background, eval=FALSE}
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.
```{r Data introduction/vectors}
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.
### Example with default values
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.
```{r Example with default values, eval=FALSE}
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.
```{r Example with distr, eval=TRUE}
antibiotics.rema$dist
antibiotics.rema$tstat
```
```{r fig.height = 5, fig.width = 5, fig.align = "center"}
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.
### Example with the one-sided p-value
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.
```{r Example with one.sided.p, eval=TRUE}
antibiotics.rema.one.pval <- rema(antibiotics.rema, one.sided.p = TRUE)
antibiotics.rema.one.pval$pval.one.sided
```
### Example with no mid-p adjustment
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.
```{r Example with mid.p, eval=TRUE}
rema(antibiotics.rema, mid.p = FALSE)
```
### Example with alpha set to 0.1
This example sets the `alpha` argument to 0.1, resulting in a 90% confidence
interval. As expected, increasing alpha decreases the confidence interval width.
```{r Example with alpha, eval=TRUE}
rema(antibiotics.rema, alpha = 0.1)
```