Git Product home page Git Product logo

orchard's Introduction

Introducing the Orchard Plot for Meta-analysis

R-CMD-check test-coverage codecov Ask Us Anything\ ! Open Source Love DOI

Citing orchaRd

To cite orchaRd 2.0 in publications one can use the following reference:

Shinichi Nakagawa, Malgorzata Lagisz, Rose E. O'Dea, Patrice Pottier, Joanna Rutkowska, Alistair M. Senior, Yefeng Yang, Daniel W.A. Noble. 2023. orchaRd 2.0: An R package for visualizing meta-analyses with orchard plots. Methods in Ecology and Evolution, https://doi.org/10.1111/2041-210X.14152 (preprint = EcoEvoRxiv, https://doi.org/10.32942/X2QC7).

For earlier versions please cite:

Nakagawa, S., Lagisz, M., O'Dea, R. E., Rutkowska, J., Yang, Y., Noble, D. W., & Senior, A. M. (2019). The Orchard Plot: Cultivating Forest Plots for Use in Ecology, Evolution and Beyond. Research Synthesis Methods https://doi.org/10.1002/jrsm.1424 12: 4-12 (preprint = EcoEvoRxiv https://doi.org/10.32942/osf.io/epqa7)

Installation

To install orchaRd use the following code in R:

install.packages("pacman")
pacman::p_load(devtools, tidyverse, metafor, patchwork, R.rsp, emmeans)

devtools::install_github("daniel1noble/orchaRd", force = TRUE)
library(orchaRd)

How to use?

We detail how to use the orchaRd package in the vignette.

Issues with orchaRd 2.0?

Please note that orchaRd 2.0 is still under active development and testing. If you use it, you should check that the results are what you expect. We do have a number of tests already in place, but there may still be situations where it fails. If you find a bug or a situation that doesn't match your expectations let us know by lodging an issue on GitHub.

orchard's People

Contributors

befriendabacterium avatar daniel1noble avatar itchyshin avatar roseodea avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar  avatar

orchard's Issues

R2 bootstrapping

Some potential issues with bootstrapping R2 for meteor models. Possible updates in metafor to mods_formula

Reordering moderators in Orchard Plots

I am trying to change the order of how the moderators are presented in the orchard plot. For example, in the following image, I want to move Chordata up, then Mollusca, etc. How can I move the moderators in the order that I want in an orchard plot?
image

i2_ml()

Yefeng: "using i2_ml() to calculate confidence intervals for an intercept-only model (rma object), we need to fit intercept-only model with argument “mods = ~ 1”. Otherwise, CIs can not be returned by i2_ml"

Interaction models - ref_grid

Some interaction models, mainly between two categorical moderators, may fail when metafor drops 'redundant predictors'. This creates a disconnect between data and model object such that emmeans cannot create a reference grid. Should add more checks in for this.

For now, a very simple solution, is to create a new categorical moderator that interacts the two moderators of interest and fit this new model. Keeping as character string will mean that factor levels don't cause issues.

'k' labels don't reorder when reordering groups/trees via +scale_x_discrete()

Describe the bug
When you reorder the groups/trees using +scale_x_discrete(), the k label ordering gets left behind/doesn't reorder to match the new order accordingly.

To Reproduce

library(orchaRd)

data(lim)

# Add in the sampling variance
lim$vi <- (1/sqrt(lim$N - 3))^2

# Lets fit a meta-regression - I will do Article non-independence. The
# phylogenetic model found phylogenetic effects, however, instead we could fit
# Phylum as a fixed effect and explore them with an Orchard Plot
lim_MR <- metafor::rma.mv(yi = yi, V = vi, mods = ~Phylum - 1, random = list(~1 |
                                                                               Article, ~1 | Datapoint), data = lim)
summary(lim_MR)

set.seed(123)
#make a vector to reorder the groups/trees
neworder<-sample(unique(lim$Phylum))
  
# Plot the meta-regression model
orchaRd::orchard_plot(lim_MR, mod = "Phylum", group = "Article", xlab = "Correlation coefficient",
                      alpha = 0.5, transfm = "tanh", angle = 45, N = "N", cb = FALSE)+
                      scale_x_discrete(limits = neworder) #reorder using the neworder vector

Expected behavior
K labels should reorder to match new groups/trees order.

Screenshots
image

Desktop (please complete the following information):

  • OS: Windows 10
  • Browser: Firefox
  • Version: 115.0.2

Additional context
Sorry to raise another issue - no real rush, I am just playing around with the package a lot and so thought I may as well report bugs/double up as a tester as I do it. Hope I am being more of a help than a thorn in your side! I will try to fix myself too just logging it here.

When flip=FALSE, 'angle' doesn't behave as expected

When flip=FALSE, the rotation confuses the 'angle' argument because it now changes the angle of the effect size axis tick marks (which are now on the y axis), rather than the grouping labels. As a workaround you can override it by using +theme() lines to manually specify the axes, but longer term it might be worth fixing this in the orchaRd function itself. For example, could have two arguments 'angle.group' and 'angle.effect' to specify them specifically in a way that is robust to flipping the orientation of the plot? Sorry don't mean to add to your plate just flagging it up for future!

Originally posted by @befriendabacterium in #28 (comment)

Prediction intervals given by mod_results() don't work with random effects specified with an '~ inner | outer' term

Describe the bug
When specifying a random effect with an '~ inner | outer' term, the prediction intervals outputted are identical to the confidence intervals (which are correct). You can see the issue by comparing to the predict() function output for models without and with '~ inner | outer' terms (see code below).

To Reproduce

library(metadat)

# PREPARE DATA (ISHAK 2007) -----------------------------------------------

### copy data into 'dat' and examine data
dat <- dat.ishak2007
head(dat, 5)
## Not run:
### load metafor package
library(metafor)
### create long format dataset
dat <- reshape(dat, direction="long", idvar="study", v.names=c("yi","vi"),
               varying=list(c(2,4,6,8), c(3,5,7,9)))
dat <- dat[order(study, time),]
### remove missing measurement occasions from dat.long
dat <- dat[!is.na(yi),]
rownames(dat) <- NULL
head(dat, 8)

# PLOT DATA ---------------------------------------------------------------

### plot data
with(dat, interaction.plot(time, study, yi, type="b", pch=19, lty="solid", xaxt="n",
                           legend=FALSE, xlab="Time Point", ylab="Mean Difference", bty="l"))
axis(side=1, at=1:4, lab=c("1 (3 months)", "2 (6 months)", "3 (12 months)", "4 (12+ months)"))

# MODEL 1: SIMPLE RANDOM EFFECT -------------------------------------------

### intercept-only model with study as a random effect
res1 <- rma.mv(yi, vi, random =  list(~1 | study),
              data = dat)
print(res1, digits=2)

#expected result
predict(res1)
 pred     se    ci.lb    ci.ub    pi.lb    pi.ub 
 -27.1854 0.9450 -29.0375 -25.3333 -37.8370 -16.5339 
#actual (orchaRd) result - identical to above expected, and prediction interval is different from confidence interval
mod_results1<-orchaRd::mod_results(model = res1, mod = "1", group='study')
mod_results1
name  estimate   lowerCL   upperCL   lowerPR   upperPR
1 Intrcpt -27.18543 -29.03754 -25.33333 -37.83702 -16.53384
# MODEL 2: RANDOM EFFECT WITH INNER/OUTER TERM-------------------------------------------

### intercept-only model with heteroscedastic AR(1) structure for the true (random) effects
res2 <- rma.mv(yi, vi, random =  list(~ time | study),
              struct = "CAR", data = dat)
print(res2, digits=2)

#expected result
predict(res2)
   pred     se    ci.lb    ci.ub    pi.lb    pi.ub 
 -27.1854 0.9450 -29.0375 -25.3333 -37.8370 -16.5338 
#actual (orchaRd) result - different from above expected, and prediction interval is identical to confidence interval
mod_results2<-orchaRd::mod_results(model = res2, mod = "1", group='study')
mod_results2
     name  estimate   lowerCL   upperCL   lowerPR   upperPR
1 Intrcpt -27.18543 -29.03754 -25.33333 -29.03754 -25.33333

Expected behavior

  • Prediction interval should not equal confidence interval
  • Predictions should be identical to those outputted by predict() - I am less sure about this as I'm not sure if the underlying statistical methods between predict() and emmeans() (used in mod_results) differ.

Desktop (please complete the following information):

  • OS: Windows
  • Browser: Firefox
  • Version: 115.0.2

Documentation/function to make an orchard plot like a forest plot

Is your feature request related to a problem? Please describe.
It may sound like an odd one, but I'm wanting to create an orchard plot that substitutes for a forest plot i.e. an information-richer plot with per-study estimates (but individual effects within studies shown) and an overall estimate at the bottom.

Describe the solution you'd like
I imagine this is possible in orchaRd already, but I haven't found specific documentation or a function to guide me on how to do this, which I'd like.

Describe alternatives you've considered
My workaround at the moment is to create an intercept-only random effects model (model 1) and a random effects model with 'study' as moderator (model 2), and then plot the two alongside one another in a similar way to a forest plot. However, I'm not sure if this is appropriate. Metafor::forest takes a random effects model without moderators and produces a forest plot with per study estimates and overall estimates - is there a way to do this in metafor.

Additional context
NA

Thanks!

`bubble_plot` error

@daniel1noble - I am putting it here - I will see what is happening.

It was working till I updated all R and RStudio. One other member reported the same error.

I will show the error when we meet

bubble_plot error

> bubble_plot(mod_ord,
+              mod = "Treat_No",
+              group = "RecNo",
+              xlab = "The number of simuli",
+              g = TRUE)
Error in `$<-.data.frame`(`*tmp*`, "condition", value = integer(0)) : 
  replacement has 0 rows, data has 640

Desktop (please complete the following information):

  • OS: Ventura 13.5

Potential Error in `pred_interval_esmeans()`

Describe the bug
In comparing the prediction intervals on some random intercept vs random slope models, I think there may be an error in the pred_interval_esmeans() function.

In the code of the function, the following is listed:

if(length(model$tau2) <= 1){ #including gamma2

I believe this ifelse() statement is supposed to determine the way to calculate the prediction intervals based on the ~inner | outer structure of the random effects in rma.mv. Because tau2 will always return a number rather than NULL, the result of length(model$tau2) will always be <=1, which sends it down the first "path" of the ifelse() statement which only takes into account the sigma(s):

sigmas <- sum(model$sigma2)
PI <- test.stat * base::sqrt(tmp$SE^2 + sigmas)

The comment next to (i.e., "including gamma2) it is what makes me think this isn't the desired function, as if multiple ~ inner arguments are included, the length of tau2 and gamma2 arguments together would be >1, which would send it down a different "path" that includes tau and gamma in the calculations which are otherwise ignored.

To Reproduce
Steps to reproduce the behavior:

library(metafor)
library(orchaRd)
library(tidyverse)

data("dat.konstantopoulos2011")
df=dat.konstantopoulos2011%>%
  as.data.frame()

#fit random intercept model
m1=rma.mv(yi ~ year, vi, random = list(~1|study, ~1|district, ~1|school), data = df, method = "REML")
e1=emmeans(emmprep(m1),~year)
p1=pred_interval_esmeans(m1,e1,mod = "year")%>%
  as.data.frame()%>%
  mutate(model="rma.mv (Intercept)")

#fit single random slope model
m2=rma.mv(yi ~ year, vi, random = list(~year|study, ~1|district, ~1|school), data = df, method = "REML")
e2=emmeans(emmprep(m2),~year)
p2=pred_interval_esmeans(m2,e2,mod = "year")%>%
  as.data.frame()%>%
  mutate(model="rma.mv (Slope x 1)")

#fit double random slope model
m3=rma.mv(yi ~ year, vi, random = list(~year|study, ~year|district, ~1|school), data = df, method = "REML")
e3=emmeans(emmprep(m3),~year)
p3=pred_interval_esmeans(m3,e3,mod = "year")%>%
  as.data.frame()%>%
  mutate(model="rma.mv (Slope x 2)")

#combine dfs
prediction.intervals=rbind(p1,p2,p3)

#proofs
tmp <- summary(e1)
tmp <- tmp[ , ]
test.stat <- stats::qt(0.975, tmp$df[[1]])

sigmas <- sum(m1$sigma2)
PI <- test.stat * base::sqrt(tmp$SE^2 + sigmas)
tmp$lower.PI <- tmp$emmean - PI
tmp$upper.PI <- tmp$emmean + PI
p1b=tmp%>%
  as.data.frame()%>%
  mutate(model="proof (Intercept)")

tmp <- summary(e2)
tmp <- tmp[ , ]
test.stat <- stats::qt(0.975, tmp$df[[1]])

sigmas <- sum(m2$sigma2)
PI <- test.stat * base::sqrt(tmp$SE^2 + sigmas)
tmp$lower.PI <- tmp$emmean - PI
tmp$upper.PI <- tmp$emmean + PI
p2b=tmp%>%
  as.data.frame()%>%
  mutate(model="proof (Slope x 1)")

tmp <- summary(e3)
tmp <- tmp[ , ]
test.stat <- stats::qt(0.975, tmp$df[[1]])

sigmas <- sum(m3$sigma2)
PI <- test.stat * base::sqrt(tmp$SE^2 + sigmas)
tmp$lower.PI <- tmp$emmean - PI
tmp$upper.PI <- tmp$emmean + PI
p3b=tmp%>%
  as.data.frame()%>%
  mutate(model="proof (Slope x 2)")

#compare proofs

final.df=rbind(p1,p1b,p2,p2b,p3,p3b)
view(final.df) # all PIs are the same

Expected behavior
I'd expect the tau and gamma values to be taken into account for models with more complex ~inner | outer structures

Desktop (please complete the following information):

  • OS: IOS
  • Rstudio 2023.6.1.524

Features to make caterpillar plot more equivalent to orchaRd plot

Is your feature request related to a problem? Please describe.
I really like that orchaRd also has a caterpillars() function to plot this alternative way of displaying the same information, but going further with this by adding some options to colour the caterpillar plot by study and change the orientation (so effect size is on y axis) would a nice feature to facilitate closer comparison.

Describe the solution you'd like

  1. Option to colour by study as in orchard_plot()
  2. Option to flip orientation of plot as in orchard_plot()

Describe alternatives you've considered
Could manually recolour and reorientate the orchard plot object I suppose but can be tricky esp. with re-colouring as the syntax required is a bit different.

Issues with df and emmeans

Noting an unanticipated problem. When models are fit without test = "t" then the model$ddf slot contains NA. This causes unanticipated issues with marginal means function as upper.CI and lower.CI cannot be computed, but are rather changed to l.asymp.CI and u.asymp.CI and predictions intervals are not calculated. This is an issue with how metafor models are working

Error in model.frame.default(formula, data = data, ...) : 'data' must be a data.frame, environment, or list

Hi,
I ran the following code (note: data was loaded, and effect sizes were calculated before running this line):


overall_all_MA <- rma.mv(yi = lnRR, V = vlnRR, random = list(~1 | experimental_id, ~1 | study),
data = data_es)
summary(overall_all_MA)
orchaRd::orchard_plot(overall_all_MA, mod = "1", data = data_es, group = "experimental_id", xlab = "log Response Ratio",
transfm = "none")


and received the following error from only the last line - running the orchard_plot() function:


Error in model.frame.default(formula, data = data, ...) :
'data' must be a data.frame, environment, or list


To check if data_es was a data frame, I used the class() function:


class(data_es)
[1] "escalc" "data.frame"


So, data_es is a data_frame, but I am getting an error that data_es is not a data frame.
Data_es contains my originally collected data and the effect size calculations.

I used the code from "Example 1" in the orchaRd vignette. I want to get an orchard plot for all my data - not just the moderator analysis. Another observation: when I ran the code from "Example 2" in the orchaRd vignette, I was able to produce sub-group analysis plots. However, I need the intercept plot (i.e. the plot for the whole dataset, without any moderators) similar to the picture below. Could you please help to resolve this issue or provide another approach to plot the data?

image

Ability to plot just the data (without the model fit) in bubble_plot()

Is your feature request related to a problem? Please describe.
Good to have this option so can easily produce equivalent plots without model fit, especially where a model wasn't fit/able to be fit/was poor/non-sig for certain subsets of the data you're presenting, and you therefore may not want to show the model fit on some plots.

Describe the solution you'd like
Option to provide a dataframe as bubble_plot()'s 'object' argument, and it to plot without this.

Describe alternatives you've considered
Could remove model fits manually via +'ing a ggplot line maybe? But doesn't solve issue where I model wasn't fit/couldn't be fit (since there's no model object to fit).

Additional context
I already made some preliminary edits to bubble_plot() to solve this, just making this feature request for housekeeping purposes/in case edits are needed.

mod_results fails because of emmeans

Problem slightly unclear:
tab1_trait_rma <- mod_results(trait_rma, mod = "Category", group = "Study.ID", data = data)
Error in ref_grid(result, ...) : Something went wrong:
Non-conformable elements in reference grid.

Fails at L100 of mod_results: grid <- emmeans::qdrg....

This could also be related to the data added to model object although data and model have same dim. Needs more investigating

error in i2_ml() with random slopes

For some reason I haven't been able to get i2_ml() to work with a rma.mv model fit with random slopes even with the struct="GEN" argument included (I am referencing the code of a paper that included this argument to presumably get around this error). The error I receive is:

"Sorry. At the moment i2_ml cannot take models with heterogeneous variance."

When the model is refit with only random intercepts, it works fine. Any help here is appreciated!

Code:
model<-rma.mv(yi, vi, mods = ~ xi, random = list(~ xi | study, ~1 | group, ~1 | es), data = df, method = "REML", test = "t",struct = "GEN")

i2_ml(model)

creating cv_ml and m1_ml

@itchyshin to create

  • cv_ml - mulitlevel version of CV (heterogeneity)
  • m1_ml - mulitlevel version of M1 (heterogeneity)

extra stuff

  • check i2_ml - it fails when the sample size is small - investigate how small

Installing orchard issue

Hi! I tried to install orchaRd using the code from the vignette, but it didn't work. I received the following error.


devtools::install_github("daniel1noble/orchaRd", force = TRUE)
Downloading GitHub repo daniel1noble/orchaRd@HEAD
Error in utils::download.file(url, path, method = method, quiet = quiet, :
download from 'https://api.github.com/repos/daniel1noble/orchaRd/tarball/HEAD' failed


Any suggestions for how to fix this issue? I'm using Windows, by the way.

make cont x cont interaction plot (bubble plot) fancy

Alistair says this

Can also draw surfaces … I did this in the BMC Bio paper - could describe as an alternative. Probably don’t want to include in the package this late in the day though

Screen Shot 2022-12-22 at 4 30 02 pm

His is referring to Panel C

Issue with orchard_plot

I'm not sure why I am having issues with orchard plots.
This is my code:
REM1=rma.mv(yi=Fishers.Zr, V=Variance.ZR, random=list(~1|Effect.size.code,~1|Experiment.Code,~1|Paper.code,~1|GENUS.SPECIES), method="REML",data=TEMP)
summary(REM1)

REM1plot <- orchard_plot(REM1, xlab = "Effect size (Zr)", transfm = "none")
REM1plot

I am getting this error

REM1plot <- orchard_plot(REM1, xlab = "Effect size (Zr)", transfm = "none")
Error in orchaRd::mod_results(object, mod = "1", group, data, N, by = by, :
Please specify the 'group' argument by providing the name of the grouping variable. See ?mod_results
Thanks

bubble_plot not working

I need to figure out why this does not work
the example in Help does work - so something else

 mod_results(decline, mod = "cYear", group = "Study_ID", data = dat)
Error in ref_grid(result, ...) : Something went wrong:
 Non-conformable elements in reference grid.
> # TODO this does not work - will get this working 
> mod_results(decline, mod = "cYear", group = "Study_ID", data = dat)
Error in ref_grid(result, ...) : Something went wrong:
 Non-conformable elements in reference grid.

Option to flip orientation of orchard plot to vertical

Is your feature request related to a problem? Please describe.
I want to make a vertically-oriented orchard plot

Describe the solution you'd like
Option in the orchard_plot() function to flip the orientation of the plot vertically (i.e. effects on Y, studies or other groupings on X axis), keeping font in standard readable orientation.

Describe alternatives you've considered
Could manually flip the whole plot outside of R, but not codable.

Additional context
It's often useful to orient effect size plots this way to align better with meta-regression bubble plots, where effects are typically on Y and an explanatory moderator is on the X. Useful to keep things in consistent orientation for the reader.

Pottier (2021) example bubble plot in vignette no longer reproducible

Describe the bug
The Pottier (2021) example bubble plot code doesn't seem to work anymore.

To Reproduce

Run the following code from https://daniel1noble.github.io/orchaRd/#example-5-bubble-plots-with-continuous-moderators or lines 556-572 in the source code at https://github.com/daniel1noble/orchaRd/blob/main/orchaRd_vignette.qmd

library(orchaRd)

# extract complete cases
data(pottier)
pottier2 <- pottier[complete.cases(pottier$body_mass), ]

# Model body mass effects
mod.body_mass <- metafor::rma.mv(yi = dARR, V = Var_dARR, mods = ~body_mass, method = "REML",
    test = "t", dfs = "contain", random = list(~1 | species_ID, ~1 | es_ID), data = pottier2)

orchaRd::bubble_plot(mod.body_mass, mod = "body_mass", group = "species_ID", xlab = "Body mass (g)",
    ylab = "Acclimation Response Ratio (ARR)", legend.pos = "top.left")

Throws error:

Error in `$<-.data.frame`(`*tmp*`, "condition", value = integer(0)) : 
  replacement has 0 rows, data has 494

Expected behavior
Produces bubble plot shown in screenshot.

Screenshots
Untitled

Desktop (please complete the following information):

  • OS: Windows
  • Browser: Firefox
  • Version: 115.0.2

Additional options for 'k.pos' in orchard_plot()

Is your feature request related to a problem? Please describe.
For the 'k.pos' argument in orchard_plot(), there are only three options (right, left and none - could do with making this clearer in the help file under Arguments too) and these aren't always sufficient to produce legible k.pos labels (see picture).

Describe the solution you'd like
Maybe some more options like 'farright' or 'farleft' (though maybe not that precise wording lol), and/or an option to specify the exact coordinates of the labels (at least on the x-axis).

Describe alternatives you've considered
Could label manually but would be good to have it integrated into the plotting function.

Additional context
See demonstrative image - argument here is set to 'right'.
image

Thanks once again, really great package for better meta-analysis plots!

How to move k away from the dashed vertical line

Hello everyone,
I have two issues with the orchard plot function. First, I am getting an error message "Error: Discrete value supplied to continuous scale" if I add xlim() to orchard_plot(). However, I then set ylim() instead then the code worked, but the code used to work smoothly with xlim() before. See the pictures below. Now I do not understand how to specify axis limits with orchard_plot(). Can anyone help me understand? Second, how do I move k away from the dashed vertical line?
Screenshot (78)
Screenshot (79)
ashed vertical line?

missing values for the intercept model may be an issue

This issue may be due to the intercept model, not meta-regression models - @daniel1noble - we can probably look at this next time together?

> modelxs <- rma.mv(yi = logX_SET,
+                   V = var_X_SET,
+                   random = ~ 1 | Study_ID,
+                   data = dat,
+                   method = "REML")
Warning message:
Rows with NAs omitted from model fitting. 
> #summary(modelxs)
> orchard_plot(modelxs, mod = "1", data = dat, group = "Study_ID", 
+              xlab = "log risk ratio")
Error in data.frame(yi, vi, moderator, stdy, type) : 
  arguments imply differing number of rows: 46, 1, 50

We have the same issue for mod_results

> modelws <- rma.mv(yi = logW_SET,
+                   V = var_W_SET,
+                   random = ~ 1 | Study_ID,
+                   data = dat,
+                   method = "REML")
Warning message:
Rows with NAs omitted from model fitting. 
> #summary(modelws)
> m2 <- mod_results(modelws, group = "Study_ID", data = dat)
Error in data.frame(yi, vi, moderator, stdy, type) : 
  arguments imply differing number of rows: 46, 1, 50

Gets fixed if you manually trim data

> dat_ws <- dat[complete.cases(dat$logW_SET), ]
> m2 <- mod_results(modelws, group = "Study_ID", data = dat_ws)
> p2 <- orchard_plot(m2, 
+              xlab = "log risk ratio  (singltons: SET)",
+              g = FALSE, k.pos = "left") + 
+   ylim(-1, 1)

i2_ml() fails when `boot` arg is not NULL if optional arg 'mods' not specified during rma model fitting

When attempting to calculate $I^2$ with confidence intervals, i.e. with a non-null boot arg in call to i2_ml(), and when model was fitted without specifying the mods argument (which is optional according to rma documentation), i2_ml() fails.

Here's a minimal example using the (not-run) English example from the documentation for i2_ml():

library(metafor)
#> Loading required package: Matrix
#> Loading required package: metadat
#> Loading required package: numDeriv
#> 
#> Loading the 'metafor' package (version 4.0-0). For an
#> introduction to the package please type: help(metafor)
library(orchaRd)
#> 
#> Loading the 'orchaRd' package (version 2.0). For an
#> introduction and vignette to the package please see: https://daniel1noble.github.io/orchaRd/

# English example
data(english)
english <- escalc(measure = "SMD", n1i = NStartControl,
                  sd1i = SD_C, m1i = MeanC, n2i = NStartExpt, sd2i = SD_E,
                  m2i = MeanE, var.names=c("SMD","vSMD"),data = english)
english_MA <- rma.mv(yi = SMD, V = vSMD,
                     random = list( ~ 1 | StudyNo, ~ 1 | EffectID), data = english)
I2_eng_1 <- i2_ml(english_MA, data = english, boot = 10)
#> Error: metafor::rma.mv(yi = ysim, V = vi, mods = mods_formula, random = random_formula, :
#> The object/variable ('mods_formula') specified for the 'mods' argument is NULL.

Created on 2023-04-17 with reprex v2.0.2

Below is a suggested re-write that allows i2_ml() to be applied when metafor::formula.rma(model, type = "mods") evaluates to NULL. I've applied it the reprex below as proof-of-concept:

Local .Rprofile detected at /Users/egould/Documents/code/ManyAnalysts/.Rprofile

library(metafor)
#> Loading required package: Matrix
#> Loading required package: metadat
#> Loading required package: numDeriv
#> 
#> Loading the 'metafor' package (version 4.0-0). For an
#> introduction to the package please type: help(metafor)
library(orchaRd)
#> 
#> Loading the 'orchaRd' package (version 2.0). For an
#> introduction and vignette to the package please see: https://daniel1noble.github.io/orchaRd/

# English example
data(english)
english <- escalc(measure = "SMD", n1i = NStartControl,
                  sd1i = SD_C, m1i = MeanC, n2i = NStartExpt, sd2i = SD_E,
                  m2i = MeanE, var.names=c("SMD","vSMD"),data = english)
english_MA <- rma.mv(yi = SMD, V = vSMD,
                     random = list( ~ 1 | StudyNo, ~ 1 | EffectID), data = english)
I2_eng_1 <- i2_ml(english_MA, data = english, boot = 10)
#> Error: metafor::rma.mv(yi = ysim, V = vi, mods = mods_formula, random = random_formula, :
#> The object/variable ('mods_formula') specified for the 'mods' argument is NULL.

# Define new fun:

i2_ml <- function (model, method = c("ratio", "matrix"), data, boot = NULL) {
  if (all(class(model) %in% c("robust.rma", "rma.mv", "rma", 
                              "rma.uni")) == FALSE) {
    stop("Sorry, you need to fit a metafor model of class robust.rma, rma.mv, rma, rma.uni")
  }
  if (any(model$tau2 > 0)) {
    stop("Sorry. At the moment i2_ml cannot take models with heterogeneous variance.")
  }
  method <- match.arg(method)
  if (method == "matrix") {
    I2s <- matrix_i2(model)
  }
  else {
    I2s <- ratio_i2(model)
  }
  if (!is.null(boot)) {
    sim <- metafor::simulate.rma(model, nsim = boot)
    mods_formula <- metafor::formula.rma(model, type = "mods")
    random_formula <- model$random
    vi <- model$vi
    pb <- progress::progress_bar$new(total = boot, format = "Bootstrapping [:bar] :percent ETA: :eta", 
                                     show_after = 0)
    if (is.null(mods_formula)) {
      I2_each <- sapply(sim, function(ysim) {
        tmp <- tryCatch(metafor::rma.mv(ysim, vi,
                                        random = random_formula, data = data))
        pb$tick()
        Sys.sleep(1/boot)
        if (method == "matrix") {
          I2 <- matrix_i2(tmp)
        }
        else {
          I2 <- ratio_i2(tmp)
        }
        return(I2)
      })
    } else{
      I2_each <- sapply(sim, function(ysim) {
        tmp <- tryCatch(metafor::rma.mv(ysim, vi, mods = mods_formula,
                                        random = random_formula, data = data))
        pb$tick()
        Sys.sleep(1/boot)
        if (method == "matrix") {
          I2 <- matrix_i2(tmp)
        }
        else {
          I2 <- ratio_i2(tmp)
        }
        return(I2)
      })
    }
    
    I2s_each_95 <- data.frame(t(apply(I2_each, 1, stats::quantile, 
                                      probs = c(0.5, 0.025, 0.975))))
    I2s <- round(I2s_each_95, digits = 3)
    colnames(I2s) = c("Est.", "2.5%", "97.5%")
  }
  return(I2s)
}

# Apply updated fun:

I2_eng_1 <- i2_ml(english_MA, data = english, boot = 10)
I2_eng_1
#>               Est.   2.5%  97.5%
#> I2_Total    74.643 69.327 81.686
#> I2_StudyNo  31.608 13.877 51.806
#> I2_EffectID 40.196 19.659 57.785

# Now check that function works when 'mods' is specified, using fish example from
# i2_ml() documentation:

## Fish example
data(fish)
warm_dat <- fish
model <- metafor::rma.mv(yi = lnrr, V = lnrr_vi,
                         random = list(~1 | group_ID, ~1 | es_ID),
                         mods = ~ experimental_design + trait.type + deg_dif + treat_end_days,
                         method = "REML", test = "t", data = warm_dat,
                         control=list(optimizer="optim", optmethod="Nelder-Mead"))
I2_fish_1 <- i2_ml(model, data = warm_dat, boot = 10)
I2_fish_1
#>               Est.   2.5%  97.5%
#> I2_Total    99.943 99.936 99.958
#> I2_group_ID 34.485 27.059 47.810
#> I2_es_ID    65.457 52.147 72.878

Created on 2023-04-17 with reprex v2.0.2

N is not working in orchard_plot

Hi, I am trying to make an orchard plot by representing sampling size instead of the default 1/SE, but it is not working and R gives me an error like: Error in [.data.frame(data, , N) : selezionate colonne non definite. I tried to run even the example given in the help pages but I obtain the same error. Any idea?

Showing the overall effect with the orchard plots

Hi! In a sub group analysis orchard plot, can you add an overall analysis in the same graph? The closest I got was the figure below, but I would like the overall analysis (which is currently in panel 2) to be in the same panel as the moderator analysis (in panel 1).

image

Option to plot non-linear meta-regressions in bubble_plot()

Is your feature request related to a problem? Please describe.
Currently there's no way to plot commonly used non-linear meta-regressions via bubble_plot().

Describe the solution you'd like
Would be good if bubble_plot() could at least handle the four non-linear options given in https://www.metafor-project.org/doku.php/tips:non_linear_meta_regression (polynomial, restricted cubic spline, natural cubic spline, or thin plate spline).

Describe alternatives you've considered
You can do this outside the function and then add your predictions a blank bubble plot without any model predictions (I added this option here: #51), but it's much preferable if the model could automatically handle it.

Additional context
I've made a start on doing this for restricted cubic splines for my own purposes here (https://github.com/befriendabacterium/orchaRd/tree/nonlinear), but needs improvement and the other non-linear options done too.

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.