[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,
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 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,
> >
> > 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
>
[[alternative HTML version deleted]]
More information about the R-sig-meta-analysis
mailing list