[R-meta] Differences in calculation of CVR in escalc()
Samuel Knapp
samuel.knapp at tum.de
Mon Oct 9 18:34:59 CEST 2017
Dear Michael and others,
thanks for your reply.
I discovered the difference: in escalc(), in the calculation of the
variance of lnCVR (vi), the subtraction of term including the
correlation (Eqn 12 in Nakagawa et al 2015) is ommited. Does anybody
know, if this Is this due due to any mathematical reasons, or just to
keep it simple? Or should/could this be adjusted in the function?
The values for vi estimated by escalc() differ from the values based on
the original equation, on average by a factor of 2.6 assessed on the
gibson example dataset.
Best regards,
Samuel
On 09/10/17 16:17, Michael Dewey wrote:
> 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
>>
>>
>
More information about the R-sig-meta-analysis
mailing list