[R-meta] extracting variances
Gram, Gil (IITA)
G@Gr@m @end|ng |rom cg|@r@org
Tue Jul 14 11:38:09 CEST 2020
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 **
[[alternative HTML version deleted]]
More information about the R-sig-meta-analysis
mailing list