There is a simple explanation:
1) The command:
factor(ftype)
does not actually turn 'ftype' permanently into a factor, since you are
not re-assigning it back to the object 'ftype'. You have to use:
ftype <- factor(ftype)
2) If you want to use the formula interface for specifying moderators, you have
to use mods = ~ <formula>, so in other words:
rma(KL, VL, mods = ~ ftype)
after you have made 'ftype' a factor (see 1).
Or you can simply use:
rma(KL, VL, mods = ~ factor(ftype))
which does the conversion of 'ftype' into a factor within the model
formula.
Best,
Wolfgang
--
Wolfgang Viechtbauer, Ph.D., Statistician
Department of Psychiatry and Psychology
School for Mental Health and Neuroscience
Faculty of Health, Medicine, and Life Sciences
Maastricht University, P.O. Box 616 (VIJV1)
6200 MD Maastricht, The Netherlands
+31 (43) 388-4170 | http://www.wvbauer.com
> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at
r-project.org]
> On Behalf Of John Hodgson
> Sent: Thursday, August 09, 2012 12:56
> To: r-help at r-project.org
> Subject: [R] Factor moderators in metafor
>
> I'm puzzled by the behaviour of factors in rma models, see example and
> comments below. I'm sure there's a simple explanation but can't
see it...
>
> Thanks for any input
>
> John Hodgson
>
>
> ------------------------------------- code/selected output ---------------
> --
>
>
> library(metafor)
>
> ## Set up data (from Lenters et al A Meta-analysis of Asbestos and
> Lung
> Cancer...
> ## Environmental Health Perspectives ? volume 119 | number 11 |
> November
> 2011)
>
> KL = c(0.02905, 0.06929, -0.1523, 1.6441, 0.1215, 0.3975, 1.0566, 0.1257,
> 0.2277, 0.06791, 0.08164, 0.2526, 0.07577, 0.03266, 0.1141, 0.1836,
> 1.8276,
> 0.4149, 15.4974)
> SE = c(0.006633, 0.09335, 0.08909, 0.4297, 0.07858, 0.1753, 0.3679,
> 0.1837,
> 0.2172, 0.2775, 0.4201, 0.1976, 0.7688, 0.06507, 0.06239, 0.09061, 0.9509,
> 0.2181, 7.331)
>
> VL = SE*SE
>
> amph = c(0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
> mix = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
> ftype = c(0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2)
>
>
> factor(amph)
> factor(ftype)
> factor(mix)
>
> ## Fit ftype...
>
> > rma(KL,VL,mods=ftype)
>
>
> Mixed-Effects Model (k = 19; tau^2 estimator: REML)
>
> tau^2 (estimate of residual amount of heterogeneity): 0.0111 (SE = 0.0095)
> tau (sqrt of the estimate of residual heterogeneity): 0.1054
>
> Test for Residual Heterogeneity:
> QE(df = 17) = 43.0937, p-val = 0.0005
>
> Test of Moderators (coefficient(s) 2):
> QM(df = 1) = 1.1069, p-val = 0.2928
>
> Model Results:
>
> estimate se zval pval ci.lb ci.ub
> intrcpt 0.0811 0.0606 1.3380 0.1809 -0.0377 0.2000
> mods 0.0473 0.0449 1.0521 0.2928 -0.0408 0.1353
>
>
> Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
>
>
> ## why only one df for the 2-level factor?
> ## in other words, why isn't the above model the same as the
> following...
>
>
> >
> >
> >
> > rma(KL,VL,mods=cbind(amph,mix))
>
> Mixed-Effects Model (k = 19; tau^2 estimator: REML)
>
> tau^2 (estimate of residual amount of heterogeneity): 0.0030 (SE = 0.0046)
> tau (sqrt of the estimate of residual heterogeneity): 0.0549
>
> Test for Residual Heterogeneity:
> QE(df = 16) = 37.5762, p-val = 0.0017
>
> Test of Moderators (coefficient(s) 2,3):
> QM(df = 2) = 6.9220, p-val = 0.0314
>
> Model Results:
>
> estimate se zval pval ci.lb ci.ub
> intrcpt 0.0380 0.0402 0.9448 0.3447 -0.0408 0.1169
> amph 0.2879 0.1163 2.4754 0.0133 0.0599 0.5158 *
> mix 0.0888 0.0625 1.4199 0.1556 -0.0338 0.2114
>
> ---
> Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1