[R-meta] Extracting meta regression elements automatically
Reza Norouzian
rnorouz|@n @end|ng |rom gm@||@com
Mon Jun 26 00:41:24 CEST 2023
Hi Danial,
I should admit that I didn't go through your message in its entirety.
However, I'm hoping that a few simple tips will help you get started
on extracting the quantities that you want from an rma.mv object
(let's call that object *fit*).
Generally, you can obtain the regression table (including the
estimate, se, zval or tval, pval, ci.lb, ci.ub) from an rma.mv object
by using:
coef(summary(fit))
Regarding the 'Test of Moderators' section, you can get the QM results
as follows (looks like you have a typo here):
data.frame(Estimate = fit$QM, Df = fit$QMdf[1], pval = fit$QMp,
row.names = "QM")
If needed, you can get the QE results in a similar fashion:
data.frame(Estimate = fit$QE, Df = nobs.rma(fit), pval = fit$QEp,
row.names = "QE")
For your type of models (additive compound symmetric), you can also
get the estimates of random variance in your true effects at each
level of hierarchy by using:
data.frame(Sigma2 = fit$sigma2, row.names = fit$s.names)
Knowing these quantities should help you shape the output of the model
in your desired format. If you come across programming questions along
the way, you can always consult R programming forums such as:
https://stackoverflow.com/
Kind regards,
Reza
On Sun, Jun 25, 2023 at 4:34 PM Daniel Gucciardi via
R-sig-meta-analysis <r-sig-meta-analysis using r-project.org> wrote:
>
> Hi all,
>
> * Reposting as I forgot to use plain text rather than HTML first time around.
>
> I'm analysing data for a 3-level meta-analysis with metafor in which
> I'd like to examine (explore) moderation effects for 6 outcome
> variables across 11 different moderators variables. I’d like to
> extract key elements of the results of these moderator analyses
> automatically rather than manually, given the large quantity of
> output.
>
>
>
> First, I’ve created the following function to execute the
> meta-regression analyses:
>
>
>
> run_meta.moderator <- function(df_list, df_names, moderators,
> scale_mods = NULL) {
>
> output_list <- list()
>
>
>
> for (i in seq_along(df_list)) {
>
> df <- df_list[[i]]
>
> df_name <- df_names[i]
>
> moderator_results <- list()
>
>
>
> for (j in seq_along(moderators)) {
>
> moderator <- moderators[[j]]
>
>
>
> # Check the number of levels in the moderator variable
>
> if (length(unique(df[[moderator]])) <= 1) {
>
> # Skip the moderator if there is only one level
>
> next
>
> }
>
>
>
> # Scale continuous moderators if specified
>
> if (moderator %in% scale_mods) {
>
> df[[moderator]] <- scale(df[[moderator]])
>
> }
>
> # Perform the meta-analytic computations for the current moderator
>
> # Wrap the computations in a tryCatch block to handle convergence issues
>
> tryCatch({
>
> anova_result <- rma.mv(yi, vi,
>
> data = df,
>
> level = 95,
>
> method = "REML",
>
> tdist = TRUE,
>
> mods = as.formula(paste("~", moderator)),
>
> random = ~ 1 | study_id/esid,
>
> control=list(rel.tol=1e-8))
>
>
>
> mods_result <- rma.mv(yi, vi,
>
> data = df,
>
> level = 95,
>
> method = "REML",
>
> tdist = TRUE,
>
> mods = as.formula(paste("~", moderator, "- 1")),
>
> random = ~ 1 | study_id/esid,
>
> control=list(rel.tol=1e-8))
>
>
>
> moderator_results[[moderator]] <- list(anova_result = anova_result,
>
> mods_result = mods_result)
>
> }, error = function(e) {
>
> # Print a warning message if convergence fails for the current moderator
>
> warning(paste("Convergence failed for moderator:", moderator,
> "in data frame:", df_name))
>
> })
>
>
>
> }
>
>
>
> # Assign the moderator results for the current data frame to the output list
>
> output_list[[df_name]] <- moderator_results
>
> }
>
>
>
> return(output_list)
>
> }
>
>
>
> # run the function for each moderator individually
>
> df_list <- list(df1 = df.psychosocial.protective,
>
> df2 = df.psychosocial.risk,
>
> df3 = df.beh_inadvertent.protective,
>
> df4 = df.beh_inadvertent.risk,
>
> df5 = df.use.protective,
>
> df6 = df.use.risk)
>
> df_names <- c("df.psychosocial.protective",
>
> "df.psychosocial.risk",
>
> "df.beh_inadvertent.protective",
>
> "df.beh_inadvertent.risk",
>
> "df.use.protective",
>
> "df.use.risk")
>
>
>
> moderators <- list("data_adjusted", "data_type", "study_design",
> "participant_category",
>
> "sport_type", "sport_level", "gender",
> "substance_type", "stage_1", "stage_2", "age")
>
> scale_mods <- c("age") # Specify the moderators to be scaled
>
> results.meta.moderator <- run_meta.moderator(df_list, df_names,
> moderators, scale_mods)
>
>
>
> This function executes successfully.
>
>
>
> Second, I'd like to extract the following elements of the output for
> results.meta.moderator to an excel file (for easy manipulation,
> formatting, etc):
>
> 1. From the 'Test of Moderators' section of 'anova_result':
>
> a. F value,
>
> b. df1 and df2,
>
> c. p-val
>
> 2. From the 'Model Results' section of 'mods_result' for each
> level of each moderator:
>
> a. estimate
>
> b. se
>
> c. ci.lb
>
> d. ci.ub
>
>
>
> Here is where I was hoping you could help, as I encounter issues when
> attempting to extract the information from ‘anova_result’ and
> ‘mods_result’. For example, here is a function which I thought could
> work but doesn’t (e.g., “Error in anova_result$QM$df1 : $ operator is
> invalid for atomic vectors”):
>
>
>
> extract_moderator_results <- function(results) {
>
> extracted_results <- list()
>
>
>
> for (df_name in names(results)) {
>
> moderator_results <- results[[df_name]]
>
> extracted_df <- data.frame()
>
>
>
> for (moderator_name in names(moderator_results)) {
>
> anova_result <- moderator_results[[moderator_name]]$anova_result
>
> mods_result <- moderator_results[[moderator_name]]$mods_result
>
>
>
> # Extract elements from 'Test of Moderators' section
>
> anova_df <- data.frame(
>
> F_value = anova_result$QM,
>
> df1 <- anova_result$QM$df1,
>
> df2 <- anova_result$QM$df2,
>
> p_value = anova_result$QMp
>
> )
>
>
>
> # Extract elements from 'Model Results' section for each level
> of the moderator
>
> mods_df <- data.frame(
>
> level = levels(mods_result$model)[[2]],
>
> estimate = coef(mods_result)$estimate,
>
> se = coef(mods_result)$se,
>
> ci_lb = confint(mods_result)[, "ci.lb"],
>
> ci_ub = confint(mods_result)[, "ci.ub"]
>
> )
>
>
>
> # Add the moderator name as a column to both data frames
>
> anova_df$moderator = moderator_name
>
> mods_df$moderator = moderator_name
>
>
>
> # Append the data frames to the extracted results data frame
>
> extracted_df <- rbind(extracted_df, anova_df, mods_df)
>
> }
>
>
>
> # Assign the extracted results for the current data frame to the output list
>
> extracted_results[[df_name]] <- extracted_df
>
> }
>
>
>
> return(extracted_results)
>
> }
>
>
>
>
>
> # Load the 'writexl' package if not already installed
>
> if (!requireNamespace("writexl", quietly = TRUE)) {
>
> install.packages("writexl")
>
> }
>
> library(writexl)
>
>
>
> # Extract the moderator results
>
> extracted_results <- extract_moderator_results(results.meta.moderator)
>
>
>
> # Save the extracted results to an Excel file
>
> output_file <- "moderator_results.xlsx"
>
> for (df_name in names(extracted_results)) {
>
> write_xlsx(extracted_results[[df_name]], path = output_file, sheet =
> df_name, append = TRUE)
>
> }
>
>
>
> Welcome your advice please. For context, I’m a novice when it comes to
> writing functions, so I suspect I’ve overlooked something simple or am
> likely out of my depth here!
>
>
>
> Cheers,
>
> Daniel
>
> _______________________________________________
> R-sig-meta-analysis mailing list @ R-sig-meta-analysis using r-project.org
> To manage your subscription to this mailing list, go to:
> https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis
More information about the R-sig-meta-analysis
mailing list