[R] Multi-line forestplots, how to in R

[Michael Dewey]

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 at 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