[R-meta] weight in rmv metafor
Viechtbauer, Wolfgang (SP)
wo||g@ng@v|echtb@uer @end|ng |rom m@@@tr|chtun|ver@|ty@n|
Mon Jun 15 13:52:16 CEST 2020
Dear Normand,
Did you read this?
http://www.metafor-project.org/doku.php/tips:weights_in_rma.mv_models
It answers 1) and 2). Let's take an even simpler example:
library(metafor)
dat <- data.frame(study = c(1,1,2,2,3,4),
id = 1:6,
yi = c(.1,.3,.2,.4,.6,.8),
vi = rep(.01,6))
dat
# studies 1 and 2 contribute 2 estimates each, studies 3 and 4 a single estimate
res <- rma.mv(yi, vi, random = ~ 1 | study/id, data=dat)
res
weights(res)
The output of this is:
> weights(res)
1 2 3 4 5 6
20.485124 20.485124 20.485124 20.485124 9.029753 9.029753
So maybe this is what you are observing - that the values along the diagonal of the weight matrix are larger for the studies with 2 estimates. But those weights are just based on the diagonal of the weight matrix. One really needs to take the whole weight matrix into consideration. As described on that page, the actual weights assigned to the estimates when computing the weight average of the estimates are the row sums of the weight matrix. You can get this with:
W <- weights(res, type="matrix")
rowSums(W) / sum(W) * 100
This yields:
> rowSums(W) / sum(W) * 100
1 2 3 4 5 6
13.34116 13.34116 13.34116 13.34116 23.31768 23.31768
So actually, the weight assigned to the first and second estimate of the studies with 2 estimates is smaller than the weight assigned to the single estimate of the studies that contribute a single estimate.
Not sure what the problem is with 3) but I don't think your data are needed to clarify the issue.
Best,
Wolfgang
>-----Original Message-----
>From: Norman DAURELLE [mailto:norman.daurelle using agroparistech.fr]
>Sent: Monday, 15 June, 2020 10:10
>To: Viechtbauer, Wolfgang (SP)
>Cc: r-sig-meta-analysis
>Subject: RE: [R-meta] weight in rmv metafor
>
>Dear all, Dr Viechtbauer, Dr Del Ponte,
>
>I read what you suggested, and went through the archives of the mailing
>list, and even though some threads were dealing with similar questions, I
>didn't really find the answers I was looking for there.
>
>I have multiple questions : 1) How exactly is calculated the weight of each
>estimate in the function rma.mv() ?
>
>2) (this question might be answered too when answering the first one but...)
>Why does that function (rma.mv) attribute more weight to estimates coming
>from studies that have multiple estimates relatively to ones that come from
>studies that only have one estimate ?
>
>3) I tried to produce a reproducible example of my code, would it be
>possible to tell me if I got things right ?
>
>my dataset is this table (I didn't know if I could attach anything) :
>
>> dat
> author year estimate std_error p_value r
>number_of_obs
>1 Ballinger et al 1988 -19.3450 5.8220 4.638e-03 -0.6511331
>17
>2 Matthieu & Norman 2020 -2.7398 0.8103 7.424e-04 -0.0920000
>1352
>3 Khangura 2011 -1.9610 1.3470 1.536e-01 -0.2298655
>40
>4 Kutcher 1990 -21.1350 3.9330 1.420e-05 -0.7321321
>27
>5 Steed 2007 -13.9700 3.2500 7.355e-04 -0.7542836
>16
>6 Steed 2007 -31.9440 4.7650 2.765e-04 -0.9301781
>9
>7 Steed 2007 -4.2780 3.0940 1.883e-01 -0.3466844
>16
>8 Sprague et al 2010 -21.9880 5.3010 7.562e-04 -0.7198487
>18
>9 Upadhaya et al 2019 -43.6170 3.3440 1.133e-06 -0.9772861
>10
>10 Khangura et al 2005 -5.3500 NA 5.000e-02 -0.4200000
>29
>11 Khangura et al 2005 -10.5600 NA 1.000e-03 -0.5700000
>29
>12 Khangura et al 2005 -9.9700 NA 1.300e-02 -0.4600000
>29
>13 Khangura et al 2005 -5.4500 NA 1.100e-02 -0.4700000
>29
>14 Khangura et al 2005 -22.7500 NA 2.800e-02 -0.4200000
>29
>15 Khangura et al 2005 -16.8300 NA 2.200e-02 -0.4300000
>29
>16 Khangura et al 2005 -9.2100 NA 3.900e-02 -0.3900000
>29
>
>V1 <- c("Ballinger et al", "Matthieu & Norman", "Khangura", "Kutcher",
>"Steed", "Steed", "Steed", "Sprague et al",
> "Upadhaya et al", "Khangura et al", "Khangura et al", "Khangura et
>al", "Khangura et al", "Khangura et al",
> "Khangura et al", "Khangura et al" )
>V2 <- c(1988, 2020, 2011, 1990, 2007, 2007, 2007, 2010, 2019, 2005, 2005,
>2005, 2005, 2005, 2005, 2005)
>V3 <- c(-19.3450, -2.7398, -1.9610, -21.1350, -13.9700, -31.9440, -4.2780, -
>21.9880, -43.6170, -5.3500,
> -10.5600, -9.9700, -5.4500, -22.7500, -16.8300, -9.2100)
>V4 <- c(5.8220, 0.8103, 1.3470, 3.9330, 3.2500, 4.7650, 3.0940, 5.3010,
>3.3440, NA, NA, NA, NA, NA, NA, NA)
>V5 <- c(4.638e-03, 7.424e-04, 1.536e-01, 1.420e-05, 7.355e-04, 2.765e-04,
>1.883e-01, 7.562e-04, 1.133e-06,
> 5.000e-02, 1.000e-03, 1.300e-02, 1.100e-02, 2.800e-02, 2.200e-02,
>3.900e-02)
>V6 <- c(-0.6511331, -0.0920000, -0.2298655, -0.7321321, -0.7542836, -
>0.9301781, -0.3466844, -0.7198487,
> -0.9772861, -0.4200000, -0.5700000, -0.4600000, -0.4700000, -
>0.4200000, -0.4300000, -0.3900000)
>V7 <- c(17, 1352, 40, 27, 16, 9, 16, 18, 10, 29, 29, 29, 29, 29, 29, 29)
>dat <- cbind(V1, V2, V3, V4, V5, V6, V7)
>dat <- as.data.frame(dat)
>dat$V1 <- as.character(dat$V1)
>dat$V2 <- as.integer(as.character(dat$V2))
>dat$V3 <- as.numeric(as.character(dat$V3))
>dat$V4 <- as.numeric(as.character(dat$V4))
>dat$V5 <- as.numeric(as.character(dat$V5))
>dat$V6 <- as.numeric(as.character(dat$V6))
>dat$V7 <- as.numeric(as.character(dat$V7))
>str(dat)
>dat <- dat %>% rename(author = "V1",
> year = "V2",
> estimate = "V3",
> std_error = "V4",
> p_value = "V5",
> r = "V6",
> number_of_obs = "V7")
>
>for (i in 1:nrow(dat)){
> if (is.na(dat$std_error[i]) == TRUE ){
> dat$std_error[i] <- abs(dat$estimate[i]) / qt(dat$p_value[i]/2,
> df=dat$number_of
>_obs-2, lower.tail=FALSE)
> }
>}
>
>res <- rma.mv(yi = dat$estimate, V = (dat$std_error)**2, random = ~ 1 |
>author, data=dat)
>coef_test(res, vcov="CR2")
>forest(res, addcred = TRUE, showweights = TRUE, header = TRUE,
> order = "obs", col = "blue", slab = dat$author)
>
>funnel(res)
>ranktest(res)
>
>The for loop that I use (and in which I use a formula that Dr Viechtbauer
>gave me in a previous answer) to calculate the standard errors for estimates
>that I didn't calculate myself gives me an error message, but still gives me
>values.
>
>The error message is :
>
>Warning messages:
>1: In data$std_error[i] <- abs(data$estimate[i])/qt(data$p_value[i]/2, :
> number of items to replace is not a multiple of replacement length
>
>and it is repeated on 7 lines.
>I am not entirely sure, but I think this comes from the fact that values are
>set to NAs before the loop is run and are replaced by values, so replacing
>an item of length 0 with an item of length 1 I believe.
>
>I hope this example can be run simply with a copy / paste, I think it
>should.
>
>Did I do things correctly ? If not what should I modify ?
>
>Thank you !
>
>Norman
More information about the R-sig-meta-analysis
mailing list