我目前在做一个项目，需要运行Leslie矩阵模拟来研究不同捕捞制度对种群的影响。我的思路是搭建一个模拟循环，针对3种不同的生物学场景，分别运行1000次模拟、每次模拟跨度50年。但我刚接触RStudio，数学基础也不算强，想请大家帮忙看看代码逻辑有没有大问题，能不能得到预期的结果？下面是我用的参数、3个基础矩阵和完整代码：

一、模拟参数定义 #

#### DEFINE PARAMETERS ####
harvest_levels <- c(31, 39, 77, 155) # number of females harvested annually 
                                    #out of 80, 100, 200, 400 total harvest
                                    # assuming females are 39% of harvest
                                    # 31 is current, assuming 80 animals harvested per year         
nSims <- 1000
nYears <- 50
Init_Min <- 2222 # min initial population
Init_Max <- 3965 # max initial population

SAD_Split <- c(0.67, 0.33) # assumed stable age distribution, we are assuming this from harvest 
harvest_weights <- c(0.67, 0.33) # weight of harvest, same as SAD (Juveniles, adults)

# Density dependence
K_total <- 4000 
K_females <- K_total * 0.39

二、3种生物学场景的基础Leslie矩阵 #

# DEFINE BASE MATRICES
# Butterflies and Rainbows Matrix 
bfar_mat <- matrix(c(0, # female offspring per juvenile
                     0.82, #female offspring per adult female
                     0.48, # juvenile survival Loforth et al. 2023
                     0.82), # Adult survival Loforth et al. 2023
                   nrow=2, byrow=T)
# Doom and Gloom Matrix
dag_mat <- matrix(c(0, # female offspring per juvenile
                    0.374, #female offspring per adult female
                    0.27, # juvenile survival Paragi et al. 1994
                    0.65), # Adult survival Paragi et al. 1994
                   nrow=2, byrow=T)
# Alternative Doom and Gloom
adag_mat <- matrix(c(0, 
                     0.374, # same as doom and gloom
                     0.83, # juvenile survival Kuntze et al. 2024
                     0.82), # adult survival Kuntze et al. 2024
                   nrow=2, byrow=T)

三、核心模拟循环函数 #

#### CREATE FUNCTION ####
run_leslie_loop <- function(base_mat,
                    harvest_levels,
                    nSims,
                    nYears,
                    Init_Min,
                    Init_Max,
                    K_females,
                    SAD_Split,
                    harvest_weights,
                    scenario_name) {
  
  # Pre-calculate growth rate
  ev <- eigen(base_mat)
  R_growth <- max(Re(ev$values))
  
  master_results <- list()
  
  # Loop over harvest levels
  for (h in harvest_levels) {
    
    sim_results <- vector("list", nSims)
    current_h_vec <- h * harvest_weights
    
    # Simulation loop
    for (s in 1:nSims) {
      
      allYears <- matrix(0, nrow = 2, ncol = nYears + 1)
      
      # Random starting population
      Rand_Start_Females <- runif(1, Init_Min, Init_Max) * 0.39
      allYears[, 1] <- Rand_Start_Females * SAD_Split
      
      for (t in 2:(nYears + 1)) {
        
        # Stochasticity
        stochastic_multiplier <- rlnorm(1, meanlog = 0, sdlog = 0.1)
        current_mat <- base_mat * stochastic_multiplier
        
        N_before <- current_mat %*% allYears[, t - 1]
        
        # Density dependence (Beverton-Holt)
        N_curr <- sum(allYears[, t - 1])
        dd <- 1 + ((R_growth - 1) * N_curr / K_females)
        
        N_after <- N_before / dd
        
        # Harvest (scaled)
        #allYears[, t] <- pmax(0, N_after - current_h_vec)
        u_eff <- min(1, h / sum(N_after))  # cap at 100%
        H_vec <- u_eff * N_after * harvest_weights
        allYears[, t] <- pmax(0, N_after - H_vec)
      }
      
      # Store total population
      sim_results[[s]] <- colSums(allYears) * 2.56
    }
    
    # Convert to tidy data frame
    temp_df <- as.data.frame(sim_results)
    names(temp_df) <- paste0("V", 1:nSims)
    
    temp_df$Year <- 0:nYears
    temp_df$HarvestLevel <- as.factor(h)
    temp_df$Scenario <- scenario_name
    
    master_results[[as.character(h)]] <- temp_df |>
      tidyr::pivot_longer(
        cols = starts_with("V"),
        names_to = "SimID",
        values_to = "Population"
      )
  }
}

我梳理的几个逻辑检查点&小建议 #

作为刚接触这类模拟的新手，你的代码逻辑整体已经很清晰了，下面几个点可以再确认下是否符合你的研究需求：

1. 密度依赖的作用时机：你用Beverton-Holt模型，把密度依赖作用在矩阵投影后的种群上（N_after = N_before / dd），这种整体缩放的方式是合理的简化。如果要更贴近实际，也可以考虑把密度依赖单独绑定在繁殖率上（比如只给矩阵里的繁殖元素乘密度依赖项），不过要看你的研究假设。
2. 捕捞策略的调整：你注释掉了固定数量捕捞，换成了按比例的有效捕捞率（避免捕捞量超过当前种群），这个调整非常贴合实际，点赞！不过要确认harvest_weights的年龄分配比例是否和实际捕捞的年龄偏好一致（你现在用了和稳定年龄分布一样的比例，这个设定如果有文献支撑就更稳妥）。
3. 随机性的设定：你给整个Leslie矩阵乘了同一个对数正态随机数，意味着繁殖率和存活率的波动是完全同步的。如果想更真实，可以考虑给矩阵的不同元素（比如繁殖率、幼体存活率、成体存活率）加独立的随机噪声，不过这个属于进阶优化，当前的简化版本已经能满足初步研究需求。
4. 初始种群的转换：你先随机生成总种群范围，再转成雌种群然后按年龄分配，这个逻辑没问题，只要确认Init_Min和Init_Max确实是总种群的范围即可。
5. 函数的返回值：目前函数最后没有return(master_results)语句，运行后不会返回结果，记得在函数末尾加上这一行哦！