#Import SAS dataset of simulated raw datasets for meta-analyses #to determine bias and coverage # Update the metafor package #update.packages(ask = FALSE, repos = "https://cran.r-project.org") #install.packages("remotes") #remotes::install_github("wviechtb/metafor") library(tidyverse) library(haven) library(metafor) #import the data fullsasdata <- read_sas("simrawpb90.sas7bdat", NULL) #sasdata <- subset(fullsasdata,Sim==1) sasdata <- subset(fullsasdata,Sim<2501) #sasdata <- subset(fullsasdata,Sim>0 & Sim<201) #sasdata$tValue <- NULL; sasdata$Sex <- as.factor(sasdata$Sex) #glimpse(sasdata) # Get unique values of Sim present in the data frame uniquesim <- unique(sasdata$Sim) # Get true values from the data frame TrueHeteroSD <- sasdata$TrueHeteroSD[1] TrueFemaleMean <- sasdata$TrueFemale[1] TrueMaleMean <- sasdata$TrueMale[1] # Convert the dataset to a list # not used in this program, which now uses for-loop processing #dat_list <- split(sasdata, sasdata$Sim) # Lists to store meta-analysis results meta_results <- list() adjusted_meta <- list() adjusted_confint <- list() # Function to perform meta-analysis using rma.uni perform_meta_analysis <- function(xxx) { # Perform meta-analysis using rma.uni meta_analysis <- rma.uni(yi = xxx$Ydelta, vi = xxx$YdeltaSEsq, mods=~xxx$Sex -1, control=list(maxiter=400), method = c("REML","DL"), level=90) # Return results return(meta_analysis) } if(FALSE) { # Function to perform selection model analysis using selmodel # This was for list processing; not used currently perform_selmodel <- function(xxx) { adjusted_meta <- selmodel(xxx, type="step", steps=(0.025)) # adjusted_meta <- selmodel(xxx, type="beta") return(adjusted_meta) } } # Define a function to apply metafor::selmodel() apply_selmodel <- function(x) { # Try processing the element tryCatch({ # This code may produce an error #result <- metafor::selmodel(x, type="beta") #result <- metafor::selmodel(x, type="step", steps=(0.025), pval=meta_data$pValue) result <- metafor::selmodel(x, type="step", steps=(0.025)) return(result) }, error = function(e) { # Handle the error (e.g., print a message) cat("selmodel failed with this:\n", conditionMessage(e), "\n") # Return nothing (skip the element) when an error occurs return(invisible()) }) } # Function to eventually get Tau2 UCL using confint on results from selmodel perform_confint <- function(xxx) { adjusted_confint <- confint(xxx, level=90, fixed=FALSE, tau2=TRUE) # adjusted_confint <- confint(xxx, level=90, fixed=FALSE, tau2=TRUE) return(adjusted_confint) } # Perform meta-analysis for each dataset for (i in uniquesim) { # Subset data for the current meta-analysis meta_data <- subset(sasdata, Sim == i) # Perform meta-analysis using rma.uni resultmeta <- perform_meta_analysis(meta_data) # Store result meta_results[[i]] <- resultmeta # Run selmodel and store result resultselmodel <- apply_selmodel(resultmeta) if (!is.null(resultselmodel)) { adjusted_meta[[i]] <- resultselmodel } #resultselmodel <- perform_selmodel(resultmeta) #adjusted_meta[[i]] <- resultselmodel # Run confint and store result if (!is.null(resultselmodel)) { resultconfint <- perform_confint(resultselmodel) adjusted_confint[[i]] <- resultconfint } } # Remove any remaining NULLs meta_results <- meta_results[sapply(meta_results, function(x) !is.null(x))] adjusted_meta <- adjusted_meta[sapply(adjusted_meta, function(x) !is.null(x))] adjusted_confint <- adjusted_confint[sapply(adjusted_confint, function(x) !is.null(x))] # The following strips out a blank line at the start of each object in adjusted_confint # Define a function to extract df from a given model extract_tau <- function(model) { df <- model$random return(df) } # Define a function to extract all tau2 from a list extract_all_df <- function(xxxx) { all_df <- lapply(xxxx, extract_tau) return(all_df) } all_tau2 <- extract_all_df(adjusted_confint) #########omit from here (processing lists separately) if(FALSE) { # Define a function to apply metafor::rma() apply_rma <- function(xxx) { rma.uni(yi=Ydelta, vi=YdeltaSEsq, mods=~Sex -1, data=xxx, control=list(maxiter=400), method=c("REML","DL"), level=90) } # Apply the function to each element of dat_list using lapply meta_results <- lapply(dat_list, apply_rma) # Define a function to apply metafor::selmodel() apply_selmodel <- function(x) { # Try processing the element tryCatch({ # This code may produce an error #result <- metafor::selmodel(x, type="beta") result <- metafor::selmodel(x, type="step", steps=(0.025)) return(result) }, error = function(e) { # Handle the error (e.g., print a message) cat("selmodel failed with this:\n", conditionMessage(e), "\n") # Return nothing (skip the element) when an error occurs return(invisible()) }) } # Clear the adjusted_meta list, to avoid any confusion adjusted_meta <- list() # Apply the function to each element of meta_results using lapply adjusted_meta <- lapply(meta_results, apply_selmodel) # Remove the NULLs adjusted_meta <- adjusted_meta[sapply(adjusted_meta, function(x) !is.null(x))] # Define a function to calculate the confidence interval for tau^2 calculate_tau2_ci <- function(model) { confint(model, tau2 = TRUE) } # Apply the function to each element of sel using lapply tau2_ci <- lapply(adjusted_meta, calculate_tau2_ci) # Define a function to extract df from a given model extract_tau <- function(model) { df <- model$random return(df) } # Define a function to extract all tau2 from a list extract_all_df <- function(xxxx) { all_df <- lapply(xxxx, extract_tau) return(all_df) } all_tau2 <- extract_all_df(tau2_ci) } ######to here ###################################################### # To get a results.txt file, un-comment sink("results.txt") sink("results.txt") # Get unique values of Sim present in the data frame uniquesim <- unique(sasdata$Sim) simsdf <- data.frame(Sims=uniquesim) cat("Number of simulations available for meta-analysis:\n") colSums(!is.na(simsdf)) cat("\n") # Extract estimates from meta-analysis results into a data matrix (metafixed) and vectors (tau2, tau2se) metafixed <- sapply(meta_results, function(x) x$beta) metafixedlb <- sapply(meta_results, function(x) x$ci.lb) metafixedub <- sapply(meta_results, function(x) x$ci.ub) metatau2 <- sapply(meta_results, function(x) x$tau2) metatau2se <- sapply(meta_results, function(x) x$se.tau2) # Split the metafixed into two parts: one for SexFemale and one for SexMale metafemale <- data.frame(SexFemale=metafixed[1,]) # Get CLs and coverage for females metafemale$lb <- metafixedlb[1,] metafemale$ub <- metafixedub[1,] metafemale$CLpm <- (metafemale$ub-metafemale$lb)/2 metafemale$TrueMean <- TrueFemaleMean metafemale$Coverage <- 0 metafemale$Coverage[metafemale$lb