[R-meta] Differences in calculation of CVR in escalc()
Michael Dewey
lists at dewey.myzen.co.uk
Mon Oct 9 16:17:32 CEST 2017
Dear Samuel
Not sure what the issue is but the code from escalc is
if (measure == "CVR") {
yi <- log(sd1i/m1i) + 1/(2 * (n1i - 1)) - log(sd2i/m2i) -
1/(2 * (n2i - 1))
vi <- 1/(2 * (n1i - 1)) + sd1i^2/(n1i * m1i^2) +
1/(2 * (n2i - 1)) + sd2i^2/(n2i * m2i^2)
}
Note you can obtain this by going
library(metafor)
sink("escalc.txt")
escalc.default
sink()
and examining escalc.txt with your favourite text editor
Michael
On 09/10/2017 13:12, Samuel Knapp wrote:
> Dear all,
>
> I am conducting a meta-analysis on the stability of crop yields. I now
> follow the approach suggeted by Nakagawa et al. (2015) approach and its
> implementation in the metafor package, which helps me a lot!
>
> As I first step I compared the estimates of the escalc function for ROM,
> VR and CVR to the actual formulas (actually I used the functions in the
> supplement of Nakagawa). Fortunately, they all yielded the same
> estimates, except for the variance estimate of CVR. I did the
> calculations on the gibson example data. The respective code (only for
> CVR) is:
>
> data <- get(data(dat.gibson2002))
>
> metadat <- escalc(measure="CVR", m1i=m1i, m2i=m2i, sd1i=sd1i, sd2i=sd2i,
> n1i=n1i, n2i=n2i, data=data)
>
> # functions from Nakagawa et al. (2015)
> Calc.lnCVR<-function(CMean, CSD, CN, EMean, ESD, EN){
> ES<-log(ESD) - log(EMean) + 1 / (2*(EN - 1)) - (log(CSD) - log(CMean)
> + 1 / (2*(CN - 1)))
> return(ES)
> }
> Calc.var.lnCVR<-function(CMean, CSD, CN, EMean, ESD, EN, Equal.E.C.Corr=T){
> if(Equal.E.C.Corr==T){
> mvcorr<-cor.test(log(c(CMean, EMean)), log(c(CSD, ESD)))$estimate
> S2<- CSD^2 / (CN * (CMean^2)) + 1 / (2 * (CN - 1)) - 2 * mvcorr *
> sqrt((CSD^2 / (CN * (CMean^2))) * (1 / (2 * (CN - 1)))) +
> ESD^2 / (EN * (EMean^2)) + 1 / (2 * (EN - 1)) - 2 * mvcorr *
> sqrt((ESD^2 / (EN * (EMean^2))) * (1 / (2 * (EN - 1))))
> } else{
> Cmvcorr<-cor.test(log(CMean), log(CSD))$estimate # corrected
> (missing log()), was "cor.test(log(EMean), (ESD))$estimate"
> Emvcorr<-cor.test(log(EMean), log(ESD))$estimate
> S2<- CSD^2 / (CN * (CMean^2)) + 1 / (2 * (CN - 1)) - 2 *Cmvcorr *
> sqrt((CSD^2 / (CN * (CMean^2))) * (1 / (2 * (CN - 1)))) +
> ESD^2 / (EN * (EMean^2)) + 1 / (2 * (EN - 1)) - 2 *Emvcorr *
> sqrt((ESD^2 / (EN * (EMean^2))) * (1 / (2 * (EN - 1))))
> }
> return(S2)
> }
>
> # compare
>
> with(data,Calc.lnCVR(m2i,sd2i,n2i,m1i,sd1i,n1i))
> metadat$yi # is the same
> # with pooled correlation
> with(data,Calc.var.lnCVR(m2i,sd2i,n2i,m1i,sd1i,n1i,Equal.E.C.Corr = T))
> metadat$vi # NOT THE SAME!!!
> # with separate correlations for E and C
> with(data,Calc.var.lnCVR(m2i,sd2i,n2i,m1i,sd1i,n1i,Equal.E.C.Corr = F))
> metadat$vi # ALSO NOT THE SAME!!!
>
>
> I checked all the equations in the Nakagawa functions and couldn't find
> any error. Also, I tried the pooled and separate correlation.
> Unfortunately, I didn't manage to access the code behind the escalc
> function in order to check the underlying calculations.
>
> Does anybody have a suggestion, what this difference could be due to?
>
> (Versions: R 3.4.2, metafor 2.0)
>
> Many thanks,
>
> Sam
>
> Reference: Nakagawa et al. , 2015. Meta-analysis of variation:
> ecological and evolutionary applications and beyond. Methods Ecol Evol
> 6, 143–152. doi:10.1111/2041-210X.12309
>
>
--
Michael
http://www.dewey.myzen.co.uk/home.html
More information about the R-sig-meta-analysis
mailing list