library(tidyverse)
library(metafor)
library(devtools)
library(clubSandwich)
load_all() # Load dev package

orchard_leave1out()

orchard_leave1out() creates an orchard plot with the results from a leave-one-out analysis. For each element left-out, the plot shows effect sizes, the model estimate, 95% confidence intervals and prediction intervals. Also, it shows the left-out effect sizes as ‘ghost points’ (this is experimental, of course).

Example 1: Eklof data

data(eklof)

# Calculate the effect size
eklof <- escalc(measure = "ROM",
                n1i = N_control,
                sd1i = SD_control,
                m1i = mean_control, 
                n2i = N_treatment,
                sd2i = SD_treatment,
                m2i = mean_treatment,
                var.names = c("lnRR", "vlnRR"),
                data = eklof)

# Add the observation level factor
eklof$Datapoint <- as.factor(seq(1, dim(eklof)[1], 1))

# Also, we can get the sample size, which we can use for weighting if we would like
eklof$N <- rowSums(eklof[, c("N_control", "N_treatment")])

Create Authors column with first author and year. This is going to be used as the group from which leave-one-out.

eklof <- eklof %>% 
  mutate(Author = paste(First.author, Publication.year, sep = ", "))
eklof_MR0 <- rma.mv(yi = lnRR,
                    V = vlnRR,
                    mods = ~ Grazer.type,
                    random = list(~1 | ExptID,
                                  ~1 | Datapoint),
                    data = eklof)

results <- mod_results(eklof_MR0, group = "Author")
results
##      name   estimate   lowerCL    upperCL   lowerPR  upperPR
## 1 Intrcpt -0.6355998 -1.106768 -0.1644312 -2.814296 1.543096

You can create the basic leave-one-out plot like this:

orchard_leave1out(model = eklof_MR0,
                  group = "Author",
                  xlab  = "lnRR")

But you may want to make changes to the plot. You can customize it using the arguments from orchard_plot():

orchard_leave1out(eklof_MR0,
                  group = "Author",
                  xlab = "lnRR",
                  trunk.size = 1.2,
                  branch.size = 1.5,
                  alpha = 0.2,
                  legend.pos = "top.out")

There are a few arguments for aesthetic changes that are specific to orchard_leave1out():

  • ylab: String. Label for the y-axis.

  • ci_lines: Logical. If TRUE, the 95% CI lines are added to the plot. Default is TRUE.

  • ci_lines_color: String. Color for the 95% CI lines. Default is “red”.

  • ghost_points: Logical. If TRUE, the points from the study left out are added to the plot. Default is TRUE.

orchard_leave1out(eklof_MR0,
                  group = "Author",
                  xlab = "lnRR",
                  trunk.size = 1.2,
                  branch.size = 1.5,
                  alpha = 0.2,
                  k.pos = "none",
                  legend.pos = "top.out",
                  ylab = "Study left out",
                  ci_lines_color = "darkorchid")

You can remove ghost_points or ci_lines if you don’t want them:

orchard_leave1out(eklof_MR0,
                  group = "Author",
                  xlab = "lnRR",
                  trunk.size = 1.2,
                  branch.size = 1.5,
                  alpha = 0.1,
                  legend.pos = "top.out",
                  ylab = "Study left out",
                  ci_lines = FALSE,
                  ghost_points = FALSE)

Example 2: BCG vaccine data from metafor

dat <- metafor::escalc(measure = "RR", ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat.bcg)

dat <- dat %>% 
  mutate(reference = paste(author, year, sep = ", ")) 
  
res <- metafor::rma(yi, vi, data = dat)

You can use transformations in orchard_leave1out() using the transfm argument, just like in orchard_plot():

# Simple orchard plot just for comparison
orchard_plot(res, group = "reference", transfm = "percentr", xlab = "%")

# Leave one out using transformation and lot of args
orchard_leave1out(res,
                  group = "reference",
                  xlab = "%",
                  ylab = "Study left out",
                  legend.pos = "bottom.out", 
                  trunk.size = 1.2,
                  branch.size = 1.3,
                  twig.size = 0,
                  transfm = "percentr",
                  k.pos = "none",
                  alpha = 0.2,
                  ci_lines_color = "black")

The plot is build on ggplot2, so there are some thing that you can modify using ggplot2 functions:

orchard_leave1out(res,
                  group = "reference",
                  xlab = "%",
                  ylab = "Study left out",
                  legend.pos = "bottom.out", 
                  trunk.size = 1.2,
                  branch.size = 1.3,
                  twig.size = 0,
                  transfm = "percentr",
                  k.pos = "none",
                  alpha = 0.2,
                  ci_lines_color = "black") +
  theme_minimal() +
  labs(title = "BCG vaccine data",
       subtitle = "Leave-one-out analysis") +
  scale_fill_discrete() + 
  scale_color_discrete()
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.

Note on colors: By default, orchard_leave1out() uses a color-blind-friendly palette. But this palette has only 20 colors. If there are more than 20 elements in ‘group’, the palette is set to Viridis.

leave_one_out()

Very often the first plot is not good enough and you want to make changes. The workflow is: run orchard_leave1out(), make changes, run it again, and repeat. But orchard_leave1out() can be slow because it fits all the meta-analytic models for the leave-one-out.

You can avoid this by running the leave-one-out analysis once, storing the results, and then passing them to orchard_leave1out().

res_loo <- leave_one_out(res, group = "reference")
res_loo
##                          name   estimate    lowerCL    upperCL   lowerPR
## 1               Aronson, 1948 -0.7070838 -1.0794010 -0.3347666 -1.902891
## 2      Ferguson & Simes, 1949 -0.6540312 -1.0081881 -0.2998743 -1.771786
## 3       Rosenthal et al, 1960 -0.6855544 -1.0495324 -0.3215765 -1.853709
## 4     Hart & Sutherland, 1977 -0.6284094 -0.9745766 -0.2822422 -1.691052
## 5  Frimodt-Moller et al, 1973 -0.7641701 -1.1400674 -0.3882728 -1.947616
## 6       Stein & Aronson, 1953 -0.7108766 -1.1033678 -0.3183853 -1.949935
## 7      Vandiviere et al, 1973 -0.6552476 -1.0089693 -0.3015259 -1.773625
## 8            TPT Madras, 1980 -0.7947676 -1.1473219 -0.4422133 -1.878150
## 9      Coetzee & Berjak, 1968 -0.7411973 -1.1266735 -0.3557211 -1.962263
## 10      Rosenthal et al, 1961 -0.6530329 -1.0141910 -0.2918748 -1.783509
## 11       Comstock et al, 1974 -0.7578586 -1.1415943 -0.3741230 -1.964187
## 12   Comstock & Webster, 1969 -0.7598121 -1.1167043 -0.4029199 -1.904874
## 13       Comstock et al, 1976 -0.7775362 -1.1411804 -0.4138919 -1.917281
##      upperPR
## 1  0.4887237
## 2  0.4637236
## 3  0.4826001
## 4  0.4342333
## 5  0.4192760
## 6  0.5281817
## 7  0.4631301
## 8  0.2886153
## 9  0.4798681
## 10 0.4774432
## 11 0.4484696
## 12 0.3852500
## 13 0.3622090

The output from leave_one_out() can be passed to orchard_leave1out with the argument loo_output. Note that orchard_leave1out() still needs the model and group arguments:

orchard_leave1out(model = res,
                  loo_output = res_loo,
                  group = "reference",
                  xlab = "lnRR",
                  ylab = "Study left out",
                  trunk.size = 1.2,
                  branch.size = 1.5,
                  alpha = 0.1,
                  legend.pos = "top.out",
                  k.pos = "none")

Using leave_one_out() alone

The function leave_one_out() returns the results, imitating the output from mod_results(). So it can be used for any other purpose, like custom plots or tables.

res_loo <- leave_one_out(res, group = "reference")

res_loo$mod_table %>% 
  ggplot(aes(x = reorder(name, -estimate),
             y = estimate,
             ymin = lowerCL,
             ymax = upperCL)) +
  geom_pointrange() +
  coord_flip() +
  geom_hline(yintercept = res$beta, linetype = 2, color = "blue") +
  geom_hline(yintercept = 0, linetype = 2, color = "black") +
  xlab("Study left out") +
  ylab("lnRR") +
  ylim(-1.5, 0.5) +
  theme_minimal() + 
  scale_color_discrete()

For more complex models

For more complex meta-analytic models, it is common to use approximate variance-covariance matrix and/or cluster robust variance estimation. Both leave_one_out() and orchard_leave1out() can handle this:

It is a little bit more complicated, but it is a solution. For example:

# Didn't try to be a correct analysis, just an example
VCV <- vcalc(vi = vlnRR,
             cluster = ExptID,
             obs = Datapoint,
             rho = 0.5,
             data = eklof)

eklof_res2 <- rma.mv(yi = lnRR,
                     V = VCV,
                     mods = ~ Grazer.type,
                     random = list(~1 | ExptID,
                                   ~1 | Datapoint),
                     data = eklof)

rob.eklof <- robust(eklof_res2, cluster = ExptID, clubSandwich = TRUE)
# Define arguments for metafor::vcalc. They are passed as strings, except for rho.
vcalc_args = list(vi = "vlnRR",
                  cluster = "ExptID",
                  obs = "Datapoint",
                  rho = 0.5)

# Define arguments for metafor::robust.
robust_args = list(cluster = "ExptID", clubSandwich = TRUE)


orchard_leave1out(rob.eklof, 
                  group = 'Author',
                  xlab = "lnRR",
                  ylab = "Study left out",
                  vcalc_args = vcalc_args,    # vcalc arguments
                  robust_args = robust_args,    # robust arguments
                  k.pos = "none",
                  alpha = 0.2,
                  legend.pos = "top.out")

The same works for leave_one_out():

loo_outputs <- leave_one_out(rob.eklof,
                             group = "Author",
                             vcalc_args = vcalc_args,
                             robust_args = robust_args)

loo_outputs
##                name   estimate   lowerCL     upperCL   lowerPR  upperPR
## 1      Nelson, 1981 -0.3037829 -1.083352  0.47578600 -3.001946 2.394380
## 2 Richardsson, 1998 -0.6107525 -1.168225 -0.05328034 -3.078728 1.857223
## 3    Eriksson, 2009 -0.6182176 -1.173935 -0.06250045 -3.070017 1.833582
## 4    Lombardo, 1997 -0.6847343 -1.301984 -0.06748459 -3.332614 1.963146
## 5       Baden, 2010 -0.6346871 -1.192214 -0.07716020 -3.030516 1.761142
## 6     Moksnes, 2008 -0.6769109 -1.095287 -0.25853520 -2.580818 1.226996
## 7      Sieben, 2010 -0.6025050 -1.164736 -0.04027394 -3.079952 1.874942
## 8     Sieben , 2010 -0.6382638 -1.203170 -0.07335789 -3.151201 1.874673

This is an internal function, but just to check that it is actually using robust with clubSandwich (there are tests to check that robust is working):

raw_output <- .run_leave1out(rob.eklof,
              group = "Author",
              vcalc_args = vcalc_args,
              robust_args = robust_args)

# show the results when leaving out Richardsson 1998
raw_output$`Richardsson, 1998`
## 
## Multivariate Meta-Analysis Model (k = 32; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed     factor 
## sigma^2.1  0.4205  0.6484     17     no     ExptID 
## sigma^2.2  0.8374  0.9151     32     no  Datapoint 
## 
## Test for Residual Heterogeneity:
## QE(df = 30) = 195.3362, p-val < .0001
## 
## Number of estimates:   32
## Number of clusters:    17
## Estimates per cluster: 1-2 (mean: 1.88, median: 2)
## 
## Test of Moderators (coefficient 2):¹
## F(df1 = 1, df2 = 14.77) = 0.4371, p-val = 0.5187
## 
## Model Results:
## 
##                       estimate      se¹     tval¹     df¹    pval¹    ci.lb¹ 
## intrcpt                -0.7434  0.3971   -1.8722   14.08   0.0821   -1.5946  
## Grazer.typegastropod    0.2497  0.3778    0.6611   14.77   0.5187   -0.5565  
##                        ci.ub¹    
## intrcpt               0.1078   . 
## Grazer.typegastropod  1.0560     
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 1) results based on cluster-robust inference (var-cov estimator: CR2,
##    approx t/F-tests and confidence intervals, df: Satterthwaite approx)

Notes and to-do