[R-meta] Differences in calculation of CVR in escalc()

[Samuel Knapp]

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


-- 
Samuel Knapp

Lehrstuhl für Pflanzenernährung
Technische Universität München
(Chair of Plant Nutrition
Technical University of Munich)

Emil-Ramann-Strasse 2
D-85354 Freising

Tel. +49 8161 71-3578	
samuel.knapp at tum.de
www.researchgate.net/profile/Samuel_Knapp


	[[alternative HTML version deleted]]