[R-meta] R-sig-meta-analysis Digest, Vol 35, Issue 7

[Viechtbauer, Wolfgang (SP)]

Also, vcov(res) returns a matrix and if the model includes moderators, then this won't be a 1x1 matrix. Another problem is that it may happen that a moderator gets dropped from the model in a bootstrapped sample due to multicollinearity. Here is an example that works:

library(metafor)

dat <- dat.collins1985b[,1:7]
dat <- escalc(measure="OR", ai=pre.xti, n1i=pre.nti, ci=pre.xci, n2i=pre.nci, data=dat)
res <- rma(yi, vi, mods = ~ year, data=dat)
res

library(boot)

boot.func <- function(dat, indices) {
   sub <- dat[indices,]
   res <- try(suppressWarnings(rma(yi, vi, mods = ~ year, data=sub)), silent=TRUE)
   if (inherits(res, "try-error") || any(res$coef.na)) NA else c(coef(res), diag(vcov(res)))
}

set.seed(1234)
res.boot <- boot(dat, boot.func, R=1000)
boot.ci(res.boot, index=c(1,3))
boot.ci(res.boot, index=c(2,4))

I just noticed that rma() was returning a bit too verbose output in a particular case, so to avoid that, you might want to install the 'devel' version of metafor:

https://github.com/wviechtb/metafor#installation

Best,
Wolfgang

>-----Original Message-----
>From: Ken Beath [mailto:ken using kjbeath.com.au]
>Sent: Tuesday, 21 April, 2020 18:20
>To: Tarun Khanna
>Cc: Michael Dewey; Viechtbauer, Wolfgang (SP); r-sig-meta-analysis using r-
>project.org
>Subject: Re: [R-meta] R-sig-meta-analysis Digest, Vol 35, Issue 7
>
>one obvious error is that your res <- rma.mv( statement uses data=dat and it
>should be data=sub and the statement above it should be sub <- df[indices,]
>
>This gives:
>
>boot.func <- function(df, indices) {
>  sub <- df[indices,]
>  res <-rma.mv(yi, vi, random = ~ 1 | trial, data=sub, mods = ~ alloc)
>  if (is.element("try-error", class(res))) NA else c(coef(res), vcov(res))
>}
>It is essential when using R that you understand what each statement is
>doing, and then you can start modifying code. One way is to drop into the
>browser and see what the commands actually do.
>
>Ken
>
>On 22 Apr 2020, at 2:06 am, Tarun Khanna <khanna using hertie-school.org> wrote:
>
>Hi Michael,
>
>Sorry. I copied the wrong link. I am using data from this link. Attaching
>the data as a csv as well. Below is the code and output. http://www.metafor-
>project.org/doku.php/tips:rma.uni_vs_rma.mv
>
>Perhaps just to be clear - I am trying to bootstrap errors in estimation
>which uses rma.mv and moderator variables.
>
>##Code
>dat <- read_csv("test.csv")
>#RE
>rma(yi, vi, data=dat)
>#CRVE
>rma.mv(yi, vi, random = ~ 1 | trial, data=dat)
>#Bootstrapping
>res <-rma.mv(yi, vi, random = ~ 1 | trial, data=dat, mods = ~ alloc)
>
>boot.func <- function(df, indices) {
>  sub <- dat[indices,]
>  res <-rma.mv(yi, vi, random = ~ 1 | trial, data=dat, mods = ~ alloc)
>  if (is.element("try-error", class(res))) NA else c(coef(res), vcov(res))
>}
>
>## Output
>Bootstrap Statistics :
>        original  bias    std. error
>t1*  -0.51792376       0           0
>t2*  -0.44786464       0           0
>t3*   0.08899601       0           0
>t4*   0.19465035       0           0
>t5*  -0.19465035       0           0
>t6*  -0.19465035       0           0
>t7*  -0.19465035       0           0
>t8*   0.26607085       0           0
>t9*   0.19465035       0           0
>t10* -0.19465035       0           0
>t11*  0.19465035       0           0
>t12*  0.31363498       0           0
>
>Error in sort.int(x, na.last = na.last, decreasing = decreasing, ...) :
>  index 0 outside bounds
>In addition: Warning messages:
>1: In sqrt(tv[, 2L]) : NaNs produced
>2: In norm.inter(z, (1 + c(conf, -conf))/2) :
>  extreme order statistics used as endpoints
>
>Tarun Khanna
>PhD Researcher
>
>Hertie School
>
>Friedrichstraße 180
>10117 Berlin ∙ Germany
>khanna using hertie-school.org ∙ www.hertie-school.org