Read Ella
Comments in-line below
On 14/08/2015 00:36, mcknight e. (em8g14) wrote:
Hello,
I am working on ecological data covering a meta-analysis on invasive species
traits.
I am not very skilled in R and would love if someone could assist me in my
production of multi-line forest plots.
The data I have is: random effects mixed model and is further divided into
subsets, I have used metafor package and here is a breakdown of my code =
#this is my main data calculation into effect sizes
MA <- escalc (measure="SMD", m1i=Invasive..mean, sd1i=sd.invasive,
n1i=N.invasive, m2i=Control.mean, sd2i=sd.control, n2i=N.control, data=dataset)
res.MA <- rma(yi,vi,data=MA,method="REML");res.MA #random-effects models ;
"HS" Viechtbauer (2005)
Why does your comment say "HS" when you have in fact asked for it to
estimate tau^2 using REML (which is the default)?
#separating data
lab <- subset (x=MA, Type.of.ex=="Lab")
field <- subset (x=MA, Type.of.ex=="Field")
You do not need to do that since rma.uni has a subset parameter as you
use later
res.MAlab <- rma(yi,vi,data=lab,method="REML");res.MAlab #random-effects models ;
"HS" Viechtbauer (2005)
res.MAfield <- rma(yi,vi,data=field,method="REML");res.MAfield #random-effects models ;
"HS" Viechtbauer (2005)
Could you not get the estimates you want by using Type.of.ex as a
moderator here with the full dataset?
From here on it is difficult to be sure what to say without more
knowledge about your dataset and the underlying science but I wonder
whether in fact you need to fit a multivariate model instead if you have
(say) k studies each giving 6 trait estimates, and so on. If that is the
way to go you need to read
?rma.mv
Please set your mailer to not send HTML as it mangles your messages and
it is more readable if you (a) put each command on a separate line (b)
use more spaces. Disk space is cheap.
res.traitlab <- rma(yi,vi,mods= ~ factor(Trait)-1,data=lab);res.traitlab #model
for traits
res.traitfield <- rma(yi,vi,mods= ~ factor(Trait)-1,data=field);res.traitfield
#model for each lab traits
res.labct <- rma(yi,vi,subset=Trait=="Consumption",data=lab);res.labct
res.labec <- rma(yi,vi,subset=Trait=="Exploitative
competition",data=lab);res.labec
res.labgr <- rma(yi,vi,subset=Trait=="Growth",data=lab);res.labgr
res.labic <- rma(yi,vi,subset=Trait=="Interference
competition",data=lab);res.labic
res.labpa <- rma(yi,vi,subset=Trait=="Predator avoidance",data=lab);res.labpa
res.labpe <- rma(yi,vi,subset=Trait=="Predator escape",data=lab);res.labpe
#model for each field traits
res.fieldct <- rma(yi,vi,subset=Trait=="Consumption",data=field);res.labct
res.fieldec <- rma(yi,vi,subset=Trait=="Exploitative
competition",data=field);res.labec
res.fieldgr <- rma(yi,vi,subset=Trait=="Growth",data=field);res.labgr
res.fieldic <- rma(yi,vi,subset=Trait=="Interference
competition",data=field);res.labic
res.fieldpa <- rma(yi,vi,subset=Trait=="Predator
avoidance",data=field);res.labpa
res.fieldpe <- rma(yi,vi,subset=Trait=="Predator escape",data=field);res.labpe
#producing a graph for lab data
estimateslab <- c(coef(res.labct), coef(res.labec), coef(res.labgr),
coef(res.labic), coef(res.labpa),coef(res.labpe))
varianceslab <- c(vcov(res.labct), vcov(res.labec), vcov(res.labgr),
vcov(res.labic), vcov(res.labpa),vcov(res.labpe))
labelslab <- c("Consumption (109)","Exploitative competition (21)","Growth (33)","Interference
competition (31)","Predator avoidance (4)","Predator escape (61)")
forest(estimateslab, varianceslab, slab=labelslab, digit=0, annotate=F, xlab="Mean
effect size",ylim=c(0,11))
#producing a graph for field data
estimatesfield <- c(coef(res.fieldct),coef(res.fieldec), coef(res.fieldgr),
coef(res.fieldic),coef(res.fieldpa),coef(res.fieldpe))
variancesfield <- c(vcov(res.fieldct),vcov(res.fieldec), vcov(res.fieldgr),
vcov(res.fieldic), vcov(res.fieldpa), vcov(res.fieldpe))
labelsfield <- c("Consumption (7)","Exploitative competition (19)","Growth (2)","Interference
competition (34)","Predator avoidance (2)","Predator escape (15)")
forest(estimatesfield, variancesfield, slab=labelsfield,digit=0,annotate=F,xlab="Mean effect
size",ylim=c(0,11)) # "psize=1" size of mean box on forest plot
addpoly(res.MAfield, row=0.2, cex=1, atransf=F, mlab="RE Model Field Studies
(79)",annotate=F)
OK, so sorry for code overload... I hope you can understand what i have done.
What i need it to produce one graph with both data sets Lab and field showing
effect sizes for each of the mentioned traits. Im not super up to scratch on R
and some of the current code was shared through a colleague, however this
person isnt great at plots.
Please please can someone help me. Im currently wasting heaps of my time and
getting no where.
Sincerely grateful
Ella
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
Michael
http://www.dewey.myzen.co.uk/home.html
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
