为含区域随机效应的Ricker模型添加拟合线与均值线的实现需求

以下是几种给Ricker模型添加区域（Clean变量）随机效应的实现方案，附带可视化代码，适配你已有的无随机效应模型基础：

1. 非线性混合效应模型（nlme包） #

适合直接拟合原始Ricker非线性形式，支持灵活的随机效应结构：

library(nlme)
library(ggplot2)

# 定义Ricker模型函数
ricker_fun <- function(Nt_minus1, lambda, alpha) {
  lambda * Nt_minus1 * exp(-alpha * Nt_minus1)
}

# 拟合模型：全局固定效应 + 区域水平的lambda、alpha随机效应
ricker_mixed <- nlme(
  formula = Nt ~ ricker_fun(Nt_minus1, lambda, alpha),
  data = dat,
  fixed = lambda + alpha ~ 1,  # 全局均值参数
  random = lambda + alpha ~ 1 | Clean,  # 区域随机偏移
  start = c(lambda = 2, alpha = 0.01)  # 初始值根据你的数据调整
)

# 计算各区域拟合值（固定+随机效应）和全局均值拟合值
ranef_vals <- ranef(ricker_mixed)
fixed_vals <- fixef(ricker_mixed)
dat$lambda_i <- fixed_vals["lambda"] + ranef_vals[dat$Clean, "lambda"]
dat$alpha_i <- fixed_vals["alpha"] + ranef_vals[dat$Clean, "alpha"]
dat$fit_i <- ricker_fun(dat$Nt_minus1, dat$lambda_i, dat$alpha_i)
dat$fit_global <- ricker_fun(dat$Nt_minus1, fixed_vals["lambda"], fixed_vals["alpha"])

# 可视化
ggplot(dat, aes(x = Nt_minus1, y = Nt)) +
  geom_point(aes(color = Clean), alpha = 0.6) +
  geom_line(aes(y = fit_i, color = Clean), linetype = "dashed") +  # 各区域拟合线
  geom_line(aes(y = fit_global), color = "black", size = 1.2) +  # 全局均值线
  labs(x = "N(t-1)", y = "N(t)", color = "区域") +
  theme_bw()

2. 对数线性化混合效应模型（lme4包） #

将Ricker模型取对数转化为线性形式，拟合效率更高，适合大数据集：

library(lme4)

# 对数线性化模型（假设响应无零值，误差正态分布）
ricker_lmer <- lmer(
  formula = log(Nt) ~ log(Nt_minus1) - Nt_minus1 + (log(Nt_minus1) - Nt_minus1 | Clean),
  data = dat,
  start = c(`(Intercept)` = log(2), `-Nt_minus1` = 0.01)
)

# 若存在过度离散，改用负二项混合模型
# ricker_glmer <- glmer.nb(
#   log(Nt) ~ log(Nt_minus1) - Nt_minus1 + (log(Nt_minus1) - Nt_minus1 | Clean),
#   data = dat,
#   start = c(`(Intercept)` = log(2), `-Nt_minus1` = 0.01)
# )

# 生成拟合值
dat$fit_i <- exp(predict(ricker_lmer, type = "response"))
dat$fit_global <- exp(predict(ricker_lmer, re.form = ~0, type = "response"))

# 可视化代码同nlme方案

3. 带随机效应的最大似然拟合（bbmle包） #

如果要延续你之前的mle2思路，可通过积分随机效应的边际似然实现：

library(bbmle)
library(mvtnorm)

# 定义负对数似然函数（整合随机效应）
ricker_mle2_ll <- function(lambda0, alpha0, sd_lambda, sd_alpha, rho, sd_resid) {
  cov_mat <- matrix(
    c(sd_lambda^2, rho*sd_lambda*sd_alpha, rho*sd_lambda*sd_alpha, sd_alpha^2),
    nrow = 2
  )
  total_ll <- 0
  for (region in unique(dat$Clean)) {
    sub_dat <- dat[dat$Clean == region, ]
    # 积分掉随机效应u1（lambda的偏移）和u2（alpha的偏移）
    region_ll <- integrate(function(u) {
      lambda <- lambda0 + u[1]
      alpha <- alpha0 + u[2]
      sum(dnorm(sub_dat$Nt, mean = lambda*sub_dat$Nt_minus1*exp(-alpha*sub_dat$Nt_minus1),
                sd = sd_resid, log = TRUE)) +
        dmvnorm(u, mean = c(0,0), sigma = cov_mat, log = TRUE)
    }, lower = c(-Inf, -Inf), upper = c(Inf, Inf))$value
    total_ll <- total_ll + region_ll
  }
  -total_ll
}

# 拟合模型（计算量较大，适合小数据集）
ricker_mle2_random <- mle2(
  ricker_mle2_ll,
  start = list(lambda0 = 2, alpha0 = 0.01, sd_lambda = 0.1,
               sd_alpha = 0.001, rho = 0, sd_resid = 0.5)
)

4. 贝叶斯非线性混合模型（brms包） #

适合需要不确定性估计的场景，灵活性强：

library(brms)

# 拟合贝叶斯Ricker混合模型
ricker_brms <- brm(
  bf(Nt ~ lambda * Nt_minus1 * exp(-alpha * Nt_minus1),
     lambda ~ 1 + (1 | Clean),
     alpha ~ 1 + (1 | Clean),
     nl = TRUE),
  data = dat,
  family = gaussian(),
  prior = c(prior(normal(2, 1), nlpar = "lambda"),
            prior(normal(0.01, 0.005), nlpar = "alpha"),
            prior(exponential(1), class = "sd")),
  chains = 4, iter = 2000
)

# 提取拟合值与置信区间
dat$fit_i <- fitted(ricker_brms, re_formula = NULL)[,1]
dat$fit_global <- fitted(ricker_brms, re_formula = NA)[,1]
global_ci <- fitted(ricker_brms, re_formula = NA)[,c(2,3)]
dat$global_lower <- global_ci[,1]
dat$global_upper <- global_ci[,2]

# 带置信区间的可视化
ggplot(dat, aes(x = Nt_minus1, y = Nt)) +
  geom_point(aes(color = Clean), alpha = 0.6) +
  geom_line(aes(y = fit_i, color = Clean), linetype = "dashed") +
  geom_line(aes(y = fit_global), color = "black", size = 1.2) +
  geom_ribbon(aes(ymin = global_lower, ymax = global_upper), alpha = 0.2, fill = "black") +
  labs(x = "N(t-1)", y = "N(t)", color = "区域") +
  theme_bw()

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