[R-meta] Covariance-variance matrix when studies share multiple treatment x control comparison - R-SIG-meta-analysis

Ju Lee

Mon, Oct 14, 2019 3:48 PM #

Dear Wolfgang, James,

Thank you for all your helpful suggestions and feedback.
Following up on previous posts, I am wondering what would be the best way to visualize the funnel plot when I am using modified egger test suggested by Pustejovsky and Rodgers (2018).

Based on our previous conversation, I am wondering if I should plot modified precision sqrt(1 / n1+ 1 / n2) against  Hedges' d (x axis), instead of SE or inverse SE against residuals....

Is there a way to do this through funnel() or funnel.rma() function? Below is my current code with egger.full being the model being run with Hedges' d and "sqrt(1 / n1+ 1 / n2)" as modified precision.

Would specifying yaxis as "yaxis="sqrtninv" equivalent to this approach? or maybe I need an alternative approach here?

Thank you very much,
Best
JU

prec.<-function(CN,TN){
  pr<-sqrt((1 / CN) + (1/TN))
  return(pr)
}
precision.2<-prec.(MHF$n.t, MHF$n.c)
head(precision.2)
egger.full<-rma.mv(hedged~precision.2, CV, data=MHF, random = list(~ 1 | Study, ~1|Id))
egger.full

funnel(egger.full, level=c(90, 95, 99),
       shade=c("white", "gray", "gray40"),back="white",
       yaxis="     ", xlab="", ylab="",
       digits=c(2,2))

From: Viechtbauer, Wolfgang (SP) <wolfgang.viechtbauer at maastrichtuniversity.nl>
Sent: Friday, September 27, 2019 8:27 AM
To: Ju Lee <juhyung2 at stanford.edu>; James Pustejovsky <jepusto at gmail.com>
Cc: r-sig-meta-analysis at r-project.org <r-sig-meta-analysis at r-project.org>
Subject: RE: Covariance-variance matrix when studies share multiple treatment x control comparison

You are still using 'residuals' as the outcome. Don't do that. Just use the Hedges' g values as the outcome. Also, you should specify the correct V matrix (I think you called it VCV in an earlier post). So, for example:

rma.mv(hedged ~ precision.2, VCV, data=MHF, random = list(~ 1 | Study, ~1 | Id),  subset=(Spatial.scale.2=="Fragmentation scale"))

I haven't looked at your other post in detail, but 'random = ~ factor(x) | Study/Id' doesn't actually work (at least not in the way you think it does). Please update your metafor installation to get an error that will inform you of this. Instead, random = list(~factor(x)|Study, ~factor(x)|Id) is indeed the correct way to specify two '~ inner | outer' terms.

Best,
Wolfgang

-----Original Message-----
From: Ju Lee [mailto:juhyung2 at stanford.edu]
Sent: Friday, 27 September, 2019 16:46
To: James Pustejovsky
Cc: Viechtbauer, Wolfgang (SP); r-sig-meta-analysis at r-project.org
Subject: Re: Covariance-variance matrix when studies share multiple treatment x control comparison

Dear James, Wolfgang,

Thank you very much for this information.
I have one question extending from this is: While I run my main mixed modes always using var-covar. matrix (to account for shared study groups within each study),
it is acceptable that my egger-like regression does not incorporate this structure, but rather just use sqrt(1 / n1 + 1 / n2)  as precision (instead of sqrt(diag(v.c.matrix)) like Wolfgang suggested as one possibility) and use p-value for precision term (precision.2 which is p=0.2382) to determine the asymmetry?

prec.<-function(CN,TN){
  pr<-sqrt((1 / CN) + (1/TN))
  return(pr)
}
precision.2<-prec.(MHF$n.t, MHF$n.c)
> egger.pr2.frag<-rma.mv(residuals~precision.2,var,data=MHF,random =list( ~ 1 | Study, ~1|Id),  subset=(Spatial.scale.2=="Fragmentation scale"))
> egger.pr2.frag

Multivariate Meta-Analysis Model (k = 285; method: REML)
Variance Components:
  estim    sqrt  nlvls  fixed  factor
sigma^2.1  0.4255  0.6523     73     no   Study
sigma^2.2  0.3130  0.5595    285     no      Id
Test for Residual Heterogeneity:
  QE(df = 283) = 1041.1007, p-val < .0001
Test of Moderators (coefficient(s) 2):
  QM(df = 1) = 1.3909, p-val = 0.2382

Model Results:
  estimate      se     zval    pval    ci.lb   ci.ub
intrcpt        0.0529  0.1617   0.3274  0.7433  -0.2640  0.3699
precision.2   -0.0668  0.0567  -1.1794  0.2382  -0.1779  0.0442
---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Thank you very much, both of you.
Best,
JU

p.s. Wolfgang, I think I figured out what went wrong with how I specified my random effects in my previous e-mail. Specifying it as random=list(~factor(x)|Study, ~factor(x)|Id) instead of random= ~factor(x)|Study/Id generates results that makes sense to me now. Please let me know if this is correct way I should be coding.

James Pustejovsky

Wed, Oct 16, 2019 7:25 AM #

Ju,

Yes, I agree that it would be best to create a funnel plot of Hedges' d
estimates versus the modified measure of precision. You can do this using
the basic funnel plot function where you feed in the vector of effect sizes
and vector of modified precisions:
funnel(x = MHF$hedged, sei = precision.2, <...additional formatting
options...>)

James

On Mon, Oct 14, 2019 at 5:48 PM Ju Lee <juhyung2 at stanford.edu> wrote:

Dear Wolfgang, James,

Thank you for all your helpful suggestions and feedback.
Following up on previous posts, I am wondering what would be the best way
to visualize the funnel plot when I am using modified egger test suggested
by Pustejovsky and Rodgers (2018).

Based on our previous conversation, I am wondering if I should plot
modified precision sqrt(1 / n1+ 1 / n2) against  Hedges' d (x axis),
instead of SE or inverse SE against residuals....

Is there a way to do this through funnel() or funnel.rma() function? Below
is my current code with egger.full being the model being run with Hedges' d
and "sqrt(1 / n1+ 1 / n2)" as modified precision.

Would specifying yaxis as "yaxis="sqrtninv" equivalent to this approach?
or maybe I need an alternative approach here?

Thank you very much,
Best
JU

prec.<-function(CN,TN){
  pr<-sqrt((1 / CN) + (1/TN))
  return(pr)
}
precision.2<-prec.(MHF$n.t, MHF$n.c)
head(precision.2)
egger.full<-rma.mv(hedged~precision.2, CV, data=MHF, random = list(~ 1 |
Study, ~1|Id))
egger.full

funnel(egger.full, level=c(90, 95, 99),
       shade=c("white", "gray", "gray40"),back="white",
       yaxis="     ", xlab="", ylab="",
       digits=c(2,2))


------------------------------
*From:* Viechtbauer, Wolfgang (SP) <
wolfgang.viechtbauer at maastrichtuniversity.nl>
*Sent:* Friday, September 27, 2019 8:27 AM
*To:* Ju Lee <juhyung2 at stanford.edu>; James Pustejovsky <jepusto at gmail.com

*Cc:* r-sig-meta-analysis at r-project.org <r-sig-meta-analysis at r-project.org

*Subject:* RE: Covariance-variance matrix when studies share multiple
treatment x control comparison

You are still using 'residuals' as the outcome. Don't do that. Just use
the Hedges' g values as the outcome. Also, you should specify the correct V
matrix (I think you called it VCV in an earlier post). So, for example:

rma.mv(hedged ~ precision.2, VCV, data=MHF, random = list(~ 1 | Study, ~1
| Id),  subset=(Spatial.scale.2=="Fragmentation scale"))

I haven't looked at your other post in detail, but 'random = ~ factor(x) |
Study/Id' doesn't actually work (at least not in the way you think it
does). Please update your metafor installation to get an error that will
inform you of this. Instead, random = list(~factor(x)|Study, ~factor(x)|Id)
is indeed the correct way to specify two '~ inner | outer' terms.

Best,
Wolfgang

-----Original Message-----
From: Ju Lee [mailto:juhyung2 at stanford.edu <juhyung2 at stanford.edu>]
Sent: Friday, 27 September, 2019 16:46
To: James Pustejovsky
Cc: Viechtbauer, Wolfgang (SP); r-sig-meta-analysis at r-project.org
Subject: Re: Covariance-variance matrix when studies share multiple
treatment x control comparison

Dear James, Wolfgang,

Thank you very much for this information.
I have one question extending from this is: While I run my main mixed
modes always using var-covar. matrix (to account for shared study groups
within each study),
it is acceptable that my egger-like regression does not incorporate this
structure, but rather just use sqrt(1 / n1 + 1 / n2)  as precision (instead
of sqrt(diag(v.c.matrix)) like Wolfgang suggested as one possibility) and
use p-value for precision term (precision.2 which is p=0.2382) to determine
the asymmetry?

prec.<-function(CN,TN){
  pr<-sqrt((1 / CN) + (1/TN))
  return(pr)
}
precision.2<-prec.(MHF$n.t, MHF$n.c)

egger.pr2.frag<-rma.mv(residuals~precision.2,var,data=MHF,random =list(

~ 1 | Study, ~1|Id),  subset=(Spatial.scale.2=="Fragmentation scale"))

egger.pr2.frag

Multivariate Meta-Analysis Model (k = 285; method: REML)
Variance Components:
  estim    sqrt  nlvls  fixed  factor
sigma^2.1  0.4255  0.6523     73     no   Study
sigma^2.2  0.3130  0.5595    285     no      Id
Test for Residual Heterogeneity:
  QE(df = 283) = 1041.1007, p-val < .0001
Test of Moderators (coefficient(s) 2):
  QM(df = 1) = 1.3909, p-val = 0.2382

Model Results:
  estimate      se     zval    pval    ci.lb   ci.ub
intrcpt        0.0529  0.1617   0.3274  0.7433  -0.2640  0.3699
precision.2   -0.0668  0.0567  -1.1794  0.2382  -0.1779  0.0442
---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Thank you very much, both of you.
Best,
JU

p.s. Wolfgang, I think I figured out what went wrong with how I specified
my random effects in my previous e-mail. Specifying it
as random=list(~factor(x)|Study, ~factor(x)|Id) instead of random=
~factor(x)|Study/Id generates results that makes sense to me now. Please
let me know if this is correct way I should be coding.

James Pustejovsky

Wed, Oct 16, 2019 7:28 AM #

Forgot to add:
I think it's also reasonable to still include the contour-enhanced options
that depict significance levels, even with the modified measure of
precision. The reason is that these contours actually *do* correspond to
the significance levels for the test of (unstandardized) mean differences,
which is exactly the test that would be reported in a primary study and
might be used as the basis for judging whether a study/result should be
published.

On Wed, Oct 16, 2019 at 9:25 AM James Pustejovsky <jepusto at gmail.com> wrote:

Ju,

Yes, I agree that it would be best to create a funnel plot of Hedges' d
estimates versus the modified measure of precision. You can do this using
the basic funnel plot function where you feed in the vector of effect sizes
and vector of modified precisions:
funnel(x = MHF$hedged, sei = precision.2, <...additional formatting
options...>)

James

On Mon, Oct 14, 2019 at 5:48 PM Ju Lee <juhyung2 at stanford.edu> wrote:

Dear Wolfgang, James,

Thank you for all your helpful suggestions and feedback.
Following up on previous posts, I am wondering what would be the best way
to visualize the funnel plot when I am using modified egger test suggested
by Pustejovsky and Rodgers (2018).

Based on our previous conversation, I am wondering if I should plot
modified precision sqrt(1 / n1+ 1 / n2) against  Hedges' d (x axis),
instead of SE or inverse SE against residuals....

Is there a way to do this through funnel() or funnel.rma() function?
Below is my current code with egger.full being the model being run with
Hedges' d and "sqrt(1 / n1+ 1 / n2)" as modified precision.

Would specifying yaxis as "yaxis="sqrtninv" equivalent to this approach?
or maybe I need an alternative approach here?

Thank you very much,
Best
JU

prec.<-function(CN,TN){
  pr<-sqrt((1 / CN) + (1/TN))
  return(pr)
}
precision.2<-prec.(MHF$n.t, MHF$n.c)
head(precision.2)
egger.full<-rma.mv(hedged~precision.2, CV, data=MHF, random = list(~ 1 |
Study, ~1|Id))
egger.full

funnel(egger.full, level=c(90, 95, 99),
       shade=c("white", "gray", "gray40"),back="white",
       yaxis="     ", xlab="", ylab="",
       digits=c(2,2))


------------------------------
*From:* Viechtbauer, Wolfgang (SP) <
wolfgang.viechtbauer at maastrichtuniversity.nl>
*Sent:* Friday, September 27, 2019 8:27 AM
*To:* Ju Lee <juhyung2 at stanford.edu>; James Pustejovsky <
jepusto at gmail.com>
*Cc:* r-sig-meta-analysis at r-project.org <
r-sig-meta-analysis at r-project.org>
*Subject:* RE: Covariance-variance matrix when studies share multiple
treatment x control comparison

You are still using 'residuals' as the outcome. Don't do that. Just use
the Hedges' g values as the outcome. Also, you should specify the correct V
matrix (I think you called it VCV in an earlier post). So, for example:

rma.mv(hedged ~ precision.2, VCV, data=MHF, random = list(~ 1 | Study,
~1 | Id),  subset=(Spatial.scale.2=="Fragmentation scale"))

I haven't looked at your other post in detail, but 'random = ~ factor(x)
| Study/Id' doesn't actually work (at least not in the way you think it
does). Please update your metafor installation to get an error that will
inform you of this. Instead, random = list(~factor(x)|Study, ~factor(x)|Id)
is indeed the correct way to specify two '~ inner | outer' terms.

Best,
Wolfgang

-----Original Message-----
From: Ju Lee [mailto:juhyung2 at stanford.edu <juhyung2 at stanford.edu>]
Sent: Friday, 27 September, 2019 16:46
To: James Pustejovsky
Cc: Viechtbauer, Wolfgang (SP); r-sig-meta-analysis at r-project.org
Subject: Re: Covariance-variance matrix when studies share multiple
treatment x control comparison

Dear James, Wolfgang,

Thank you very much for this information.
I have one question extending from this is: While I run my main mixed
modes always using var-covar. matrix (to account for shared study groups
within each study),
it is acceptable that my egger-like regression does not incorporate this
structure, but rather just use sqrt(1 / n1 + 1 / n2)  as precision (instead
of sqrt(diag(v.c.matrix)) like Wolfgang suggested as one possibility) and
use p-value for precision term (precision.2 which is p=0.2382) to determine
the asymmetry?

prec.<-function(CN,TN){
  pr<-sqrt((1 / CN) + (1/TN))
  return(pr)
}
precision.2<-prec.(MHF$n.t, MHF$n.c)

egger.pr2.frag<-rma.mv(residuals~precision.2,var,data=MHF,random

=list( ~ 1 | Study, ~1|Id),  subset=(Spatial.scale.2=="Fragmentation
scale"))

egger.pr2.frag

Multivariate Meta-Analysis Model (k = 285; method: REML)
Variance Components:
  estim    sqrt  nlvls  fixed  factor
sigma^2.1  0.4255  0.6523     73     no   Study
sigma^2.2  0.3130  0.5595    285     no      Id
Test for Residual Heterogeneity:
  QE(df = 283) = 1041.1007, p-val < .0001
Test of Moderators (coefficient(s) 2):
  QM(df = 1) = 1.3909, p-val = 0.2382

Model Results:
  estimate      se     zval    pval    ci.lb   ci.ub
intrcpt        0.0529  0.1617   0.3274  0.7433  -0.2640  0.3699
precision.2   -0.0668  0.0567  -1.1794  0.2382  -0.1779  0.0442
---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Thank you very much, both of you.
Best,
JU

p.s. Wolfgang, I think I figured out what went wrong with how I specified
my random effects in my previous e-mail. Specifying it
as random=list(~factor(x)|Study, ~factor(x)|Id) instead of random=
~factor(x)|Study/Id generates results that makes sense to me now. Please
let me know if this is correct way I should be coding.