-----Original Message-----
From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces at r-project.org]
On Behalf Of Gram, Gil (IITA)
Sent: Tuesday, 14 July, 2020 11:38
To: r-sig-meta-analysis at r-project.org
Subject: [R-meta] extracting variances
Dear all,
I have a question regarding extracting the variances from my model.
Say I want to analyse the yields (tonnes per hectare) of 4 treatments
(control, OR, MR, ORMR) running across different sites and times. A
simplified version of my model would then be:
dat = escalc(measure="MN", mi=yield, sdi=sdYield, ni=nRep, data=temp)
dat$yi = sqrt(dat$yi) # sqrt transformation
dat$vi = dat$vi/(4*dat$yi) # variance adjustment to sqrt transformation
mod = rma.mv(yi, as.matrix(vi), method = 'REML', struct="HCS", sparse=TRUE,
data=dat,
mods = ~ rateORone + kgMN + I(rateORone^2) +
I(kgMN^2)
+ rateORone:kgMN + I(rateORone^2):I(kgMN^2) +
[?],
random = list(~1|ref, ~1|idRow, ~ treatment |
idSite, ~ treatment | idSite.time))
I?m interested in the yield variance responses over time, of OR and ORMR
versus control. So I extract the variance-covariance matrix H = mod$H:
Control MR OR ORMR
Control 0.1098190 0.1179042 0.1055471 0.1216751
MR 0.1179042 0.1360579 0.1174815 0.1354332
OR 0.1055471 0.1174815 0.1090329 0.1212389
ORMR 0.1216751 0.1354332 0.1212389 0.1449001
The variance responses I then calculate with e.g. responseOR = varianceOR +
varianceControl - 2*covar(OR, Control):
resOR
= (H['OR','OR'] + H['Control','Control'] - 2*H['Control','OR'])
= 0.1090329 + 0.1098190 - 2* 0.1055471
~ 0.00775
resORMR
~ 0.0114
I understand therefore that the variance responses over time for treatments
OR and ORMR are about 0.77% and 1.1%. These values are extremely small,
hence my questions:
- Am I correct that this was the correct way to estimate the yield
variability (responses) over time?
If this is all correct, then this means that there is hardly any variability
associated with these components. And one could start wondering what the
point is of even looking at this. I tried looking at the values of the other
components, and see whether these are larger.
- Keeping in mind the original data was sqrt transformed, can these values
still be considered as variances? or as standard deviations instead?
- If this makes up so little variance, then where is the variance coming
from? How much variability is associated with the error term? Or the other
components. Are these then magnitudes larger? How do I check if the sum of
all variance components equals 100% with the model output below?
I hope my questions are clear?
Thanks a lot in advance for your help,
Gil
------
Multivariate Meta-Analysis Model (k = 1161; method: REML)
Variance Components:
estim sqrt nlvls fixed factor
sigma^2.1 0.0604 0.2458 40 no ref
sigma^2.2 0.0285 0.1688 1161 no idRow
outer factor: idSite (nlvls = 71)
inner factor: treatment (nlvls = 4)
estim sqrt k.lvl fixed level
tau^2.1 0.1285 0.3584 275 no Control
tau^2.2 0.0952 0.3086 374 no MR
tau^2.3 0.1217 0.3488 234 no OR
tau^2.4 0.0711 0.2666 278 no ORMR
rho 0.7172 no
outer factor: idSite.time (nlvls = 271)
inner factor: treatment (nlvls = 4)
estim sqrt k.lvl fixed level
gamma^2.1 0.1098 0.3314 275 no Control
gamma^2.2 0.1361 0.3689 374 no MR
gamma^2.3 0.1090 0.3302 234 no OR
gamma^2.4 0.1449 0.3807 278 no ORMR
phi 0.9646 no
Test for Residual Heterogeneity:
QE(df = 1151) = 501266.0717, p-val < .0001
Test of Moderators (coefficients 2:10):
QM(df = 9) = 441.0373, p-val < .0001
Model Results:
estimate se zval pval ci.lb
ci.ub
intrcpt 1.2855 0.0691 18.6010 <.0001 1.1501
1.4210 ***
rateORone 0.0059 0.0007 8.5224 <.0001 0.0045
0.0072 ***
kgMN 0.0096 0.0009 10.5108 <.0001 0.0078
0.0114 ***
I(rateORone^2) -0.0000 0.0000 -5.2103 <.0001 -0.0000 -
0.0000 ***
I(kgMN^2) -0.0000 0.0000 -6.6753 <.0001 -0.0000 -
0.0000 ***
[?]
rateORone:kgMN -0.0000 0.0000 -3.7035 0.0002 -0.0000 -
0.0000 ***
I(rateORone^2):I(kgMN^2) 0.0000 0.0000 2.5775 0.0100 0.0000
0.0000 **