Skip to content

[R-meta] Extracting meta regression elements automatically

3 messages · Reza Norouzian, Daniel Gucciardi

#
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
#
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 at r-project.org> wrote:
1 day later
#
Many thanks for the prompt response and information Reza - your tips were
exactly what I needed to achieve my goal! For those interested, here is the
final function:

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

        # Extract the coefficients table from the summary
        coef_table <- coef(summary(mods_result))

        # Create the moderator results data frame for the current moderator
        moderator_df <- data.frame(
          Data_Frame = df_name,
          Moderator = moderator,
          Estimate = anova_result$QM,
          Df1 = anova_result$QMdf[1],
          Df2 = anova_result$QMdf[2],
          p_value = anova_result$QMp,
          estimate = coef_table[, "estimate"],
          se = coef_table[, "se"],
          ci.lb = coef_table[, "ci.lb"],
          ci.ub = coef_table[, "ci.ub"],
          stringsAsFactors = FALSE
        )

        # Store the moderator results for the current data frame
        moderator_results[[moderator]] <- moderator_df

      }, 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
  }

  # Merge results for all data frames into a single data frame
  output_df <- do.call(rbind, unlist(output_list, recursive = FALSE))

  return(output_df)
}

# 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)

print(results.meta.moderator)

# Convert the data frame to Excel format
write_xlsx(results.meta.moderator, path = "meta.regression.xlsx")
write.csv(results.meta.moderator, file = "meta.regression.csv", row.names =
TRUE)

Cheers,
Daniel
On Mon, Jun 26, 2023 at 6:41?AM Reza Norouzian <rnorouzian at gmail.com> wrote: