# [R-meta] Extracting meta regression elements automatically

Daniel Gucciardi d@n|e|@|@gucc|@rd| @end|ng |rom gm@||@com
Tue Jun 27 00:45:57 CEST 2023

```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,
df5 = df.use.protective,
df6 = df.use.risk)
df_names <- c("df.psychosocial.protective",
"df.psychosocial.risk",
"df.use.protective",
"df.use.risk")

"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 using gmail.com> wrote:

> 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,
> >
> >
> >
> >                 df5 = df.use.protective,
> >
> >                 df6 = df.use.risk)
> >
> > df_names <- c("df.psychosocial.protective",
> >
> >               "df.psychosocial.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)
> >
> > }
> >
> >
> >
> >
> >
> >
> > 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)
> >
> > }
> >
> >
> >
> > 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
>

[[alternative HTML version deleted]]

```