在上一篇博文中(【新提醒】科学网—从抽球到群落：基于蒙特卡洛零模型的构建解析（1）无资源限制 - 杜昕的博文)，我从一个极简的抽球过程出发，构建了一个无资源限制条件下的抽样—增长型零模型。在这一模型中，群落中个体数量的增加被视为一系列离散的增长事件；在每一轮中，系统仅发生一次“增加 1 个体”的事件，并按照当前物种的相对多度，将这一增长机会随机分配给某一物种。

这一设定隐含了一个关键假设：在缺乏额外生物学信息的前提下，不同物种被视为具有相同的内禀增长率。在零模型语境中，这可被理解为一种中性处理（neutral assumption / neutral theory），即不预设任何物种间增长策略差异，而将所有物种置于同一增长节奏之下。正因如此，“每一轮增加 1 个体”在不同物种之间具有统一的计数尺度，增长事件本身可被视为一个共享的时间基准；在此框架下，系统从初始群落状态演化到给定观测状态所经历的迭代轮数，便可直接由总个体数的变化倒推出，并与蒙特卡洛模拟结果进行对照，从而为后续引入更复杂设定提供一个最简且可比较的起点。

这一假设的意义在于零模型层面的可对照性，而非生物学上的同质性。通过将所有物种置于同一增长节奏之下，我们刻意剥离了增长策略差异、资源限制与种间作用等复杂因素，使得模型所刻画的仅是：在无约束增长背景下，仅由初始状态与随机分配机制所能产生的群落构成范围。

然而，这一极简处理也自然引出了下一个问题：如果不同物种在获得增长机会的能力上本就存在系统性差异，那么“增长事件”是否仍然可以被理解为具有统一尺度的时间刻度？

在经典生态学框架中，内禀增长率差异常被视为不同生活史策略的一种体现（例如通常所说的 r 型与 K 型策略）。需要注意的是，这类差异本身并不预设任何普适意义上的优势或劣势，其生态后果高度依赖于环境背景与时间尺度。在无资源限制的情形下，讨论增长率差异，并不等同于引入资源竞争或密度制约，而只是意味着：增长事件在物种之间的分配规则发生了改变。

从抽样—增长零模型的角度看，这一改变可以被理解为：在原有“按当前多度随机分配增长事件”的基础上，引入物种特异的权重，使得不同物种在同一轮增长中获得增长机会的概率不再完全相同。此时，模型仍然不涉及资源上限或显式种间作用，但增长事件已不再是一个对所有物种完全等价的时间单位。

因此，本文接下来的讨论将聚焦于这样一种扩展情形：在保持无资源限制与零模型结构不变的前提下，引入内禀增长率差异，考察增长策略差异本身能够在多大程度上影响群落构成的概率分布。这一分析并非试图判断不同策略的“优劣”，而是为后续逐步引入资源限制与密度依赖机制，提供一个清晰、可对照的基线。

———————————————————————————————————————————————————————————

内禀增长率的倍差与增长事件的权重化

物种增长模型的形式如下：

其中，r为内禀增长率。因此，在构建权重时，可将物种的特异性权重设置为exp(r_i)。其中，i是指示物种的下标。

———————————————————————————————————————————————————————————

#基于上述思想对上篇博文中的代码进行修改，修改部分加黑标注。

# 初始丰度（建议命名，输出更清晰）

N0 <- c(A = 2, B = 1, C = 1, D = 2, E = 1)  # 总数=7

# **物种差异参数：内禀增长率（可为相对值）**

r  <- c(A = 0.20, B = 0.05, C = 0.10, D = 0.30, E = 0.05)  # **新增**

# 循环次数：总共发生多少轮增长事件

n_steps <- 15          # 如果顺利，总数从 7 增到 22

# 蒙特卡洛重复次数

n_mc <- 99

simulate_community_mc <- function(N0, r, n_steps, n_mc) {  # **修改：加入 r**

  N0 <- as.integer(N0)

  r  <- as.numeric(r)                                     # **新增**

  stopifnot(length(N0) == length(r))                      # **新增**

  stopifnot(sum(N0) > 0, n_steps > 0, n_mc > 0) 

  S <- length(N0)

  sp_names <- names(N0)

  if (is.null(sp_names)) sp_names <- paste0("sp", seq_len(S))

  # **将 r 映射为正权重（增长机会权重）**

  w <- exp(r)                                             # **新增**

  names(w) <- sp_names                                    # **新增（对齐物种名）**

  out <- matrix(NA_integer_, nrow = n_mc, ncol = S,

                dimnames = list(NULL, sp_names)) 

  for (m in seq_len(n_mc)) {  

    N <- N0 

    for (t in seq_len(n_steps)) { 

      # **按（当前多度 × 物种权重）抽一个物种**

      weights <- N * w                                    # **修改**

      i <- sample.int(S, size = 1, prob = weights)         # **修改** 

      N[i] <- N[i] + 1

    }  

    out[m, ] <- N

  }

  out

}

#输出结果

result <- simulate_community_mc(N0, r, n_steps, n_mc)      # **修改：传入 r**

head(result)

————————————————————————————————————————————————

#绘图

# ============================================================

# Monte Carlo null model: sampling–growth process

# Output: mean trajectory + 95% envelope (no spaghetti lines)

# ============================================================

# -----------------------------

# 0. Settings

# -----------------------------

set.seed(123)

# 初始群落（可自行改）

N0 <- c(A = 2, B = 1, C = 1, D = 2, E = 1)  # total = 7

# 增长轮数（事件数）

n_steps <- 15   # total becomes 22

# Monte Carlo 次数

n_mc <- 99

# -----------------------------

# 1. Simulation function

# -----------------------------

simulate_paths_mc <- function(N0, n_steps, n_mc, w = NULL) {

  N0 <- as.integer(N0)

  stopifnot(sum(N0) > 0, n_steps > 0, n_mc > 0)

  S <- length(N0)

  sp_names <- names(N0)

  if (is.null(sp_names)) sp_names <- paste0("sp", seq_len(S))

  # 相同内禀增长率：w = 1

  if (is.null(w)) {

    w <- rep(1, S)

  } else {

    w <- as.numeric(w)

    stopifnot(length(w) == S)

  }

  names(w) <- sp_names

  # 保存每一步

  paths <- array(

    NA_integer_,

    dim = c(n_steps + 1, S, n_mc),

    dimnames = list(step = 0:n_steps, species = sp_names, sim = NULL)

  )

  for (m in seq_len(n_mc)) {

    N <- N0

    paths[1, , m] <- N

    for (t in seq_len(n_steps)) {

      probs <- N * w

      i <- sample.int(S, size = 1, prob = probs)

      N[i] <- N[i] + 1

      paths[t + 1, , m] <- N

    }

  }

  paths

}

# -----------------------------

# 2. Summarize: mean + 95% envelope

# -----------------------------

summarize_paths <- function(paths) {

  steps <- as.integer(dimnames(paths)$step)

  sp_names <- dimnames(paths)$species

  out <- lapply(sp_names, function(sp) {

    mat <- paths[, sp, ]   # [step, sim]

    data.frame(

      step    = steps,

      species = sp,

      mean    = apply(mat, 1, mean),

      q025    = apply(mat, 1, quantile, probs = 0.025),

      q975    = apply(mat, 1, quantile, probs = 0.975)

    )

  })

  do.call(rbind, out)

}

# -----------------------------

# 3. Run simulation

# -----------------------------

paths <- simulate_paths_mc(N0, n_steps, n_mc)

# 若考虑物种差异（下一篇再用）：

# r <- c(A = 0.20, B = 0.05, C = 0.10, D = 0.30, E = 0.05)

# paths <- simulate_paths_mc(N0, n_steps, n_mc, w = exp(r))

df_stats <- summarize_paths(paths)

# -----------------------------

# 4. Plot (white board)

# -----------------------------

library(ggplot2)

ggplot(df_stats, aes(x = step)) +

  # 95% 包络（虚线）

  geom_line(aes(y = q025), linetype = "dashed", linewidth = 0.6) +

  geom_line(aes(y = q975), linetype = "dashed", linewidth = 0.6) +

  # 均值轨迹（实线）

  geom_line(aes(y = mean), linewidth = 1.1) +

  facet_wrap(~ species, scales = "free_y") +

  labs(

    x = "Growth step (event count)",

    y = "Abundance",

    title = "Monte Carlo mean trajectory with 95% envelope",

    subtitle = "Solid line: mean; Dashed lines: 2.5%–97.5% quantiles"

  ) +

  theme_classic() +

  theme(

    panel.grid = element_blank(),

    strip.background = element_blank()

  )