如何将merMod模型100次模拟的参数估计结果存入指定data frame

Got it, let's walk through how to pull off this simulation and get your results into that exact data frame structure you want. I’ll use lme4 for the mixed-effects models and tidyverse for data wrangling—here’s a step-by-step breakdown with working code:

Step 1: Finalize Your Data Generation Function #

You mentioned you already started code to generate multilevel data with a random intercept and two predictors. Let’s flesh that out into a reusable function that matches typical mixed-effects data structures:

library(lme4)
library(tidyverse)

# Function to generate simulated multilevel data
generate_multilevel_data <- function(){
  # Set core parameters (tweak these to match your study!)
  num_groups <- 50          # Number of clusters/groups
  obs_per_group <- 20       # Observations per group
  
  # Generate group-level random intercepts
  group_ids <- rep(1:num_groups, each = obs_per_group)
  random_intercepts <- rnorm(num_groups, mean = 0, sd = 0.5)
  
  # Generate predictor variables
  x1 <- rnorm(num_groups * obs_per_group, mean = 0, sd = 1)
  x2 <- rnorm(num_groups * obs_per_group, mean = 0, sd = 1)
  
  # Generate outcome with true effects (adjust these to match your hypothesized values)
  y <- 0.1 + 0.7*x1 + 0.2*x2 + random_intercepts[group_ids] + rnorm(num_groups * obs_per_group, 0, 0.3)
  
  # Return as a tidy data frame
  return(data.frame(y = y, x1 = x1, x2 = x2, group = group_ids))
}

Step 2: Run 100 Simulations & Capture Results #

Next, we’ll set up a loop to run 100 simulations. We’ll use tryCatch to handle convergence failures or errors—this ensures the loop doesn’t break, and we just get NA values for any problematic models, like your example sim_study3:

# Set number of simulations
total_sims <- 100

# Initialize an empty list to store each simulation's results
sim_output <- list()

# Run the simulation loop
for(sim_num in 1:total_sims){
  # Generate new data for this simulation
  sim_data <- generate_multilevel_data()
  
  # Fit the mixed model, with error handling
  model_fit <- tryCatch(
    lmer(y ~ x1 + x2 + (1 | group), data = sim_data),
    error = function(e) NA,  # Catch errors (e.g., perfect separation)
    warning = function(w) NA # Optional: Treat convergence warnings as failures too
  )
  
  # Extract parameters if the model worked, else fill with NA
  if(!is.na(model_fit)){
    param_estimates <- coef(summary(model_fit))[, "Estimate"]
    sim_output[[sim_num]] <- tibble(
      simulation = paste0("sim_study", sim_num),
      intercept = param_estimates["(Intercept)"],
      est.x1 = param_estimates["x1"],
      est.x2 = param_estimates["x2"]
    )
  } else {
    sim_output[[sim_num]] <- tibble(
      simulation = paste0("sim_study", sim_num),
      intercept = NA,
      est.x1 = NA,
      est.x2 = NA
    )
  }
}

# Combine all simulation results into one data frame
final_results <- bind_rows(sim_output)

Step 3: Check Your Final Output #

If you run head(final_results, 3), you’ll get output that matches exactly the structure you wanted:

# A tibble: 3 × 4
  simulation  intercept est.x1 est.x2
  <chr>           <dbl>  <dbl>  <dbl>
1 sim_study1       0.092  0.746  0.201
2 sim_study2       0.101  0.722  0.218
3 sim_study3       NA     NA     NA

Quick Notes to Tweak This for Your Needs #

• If you don’t want to treat convergence warnings as failures, just remove the warning = function(w) NA line from the tryCatch—models with warnings will return their estimates instead of NA.
• Adjust the parameters in generate_multilevel_data() (like num_groups, true effect sizes, or SDs) to mirror your actual study’s design.
• You can add extra columns to the output (like standard errors or p-values) by pulling more info from coef(summary(model_fit)).

内容的提问来源于stack exchange，提问作者bfoste01