R语言非线性模型拟合点置信区间高效生成方法咨询

• 核心原理：你使用的y = a * exp(b * conc)指数模型属于可解析求导的形式，直接用delta方法计算预测值的方差即可，不需要调用效率极低的predictNLS()。
  方差计算公式为：Var(y_pred) = (∂y/∂a)²*Var(a) + (∂y/∂b)²*Var(b) + 2*(∂y/∂a)*(∂y/∂b)*Cov(a,b)
  对应偏导可直接手动推导：∂y/∂a = exp(b*conc)，∂y/∂b = a*conc*exp(b*conc)
  方差开方得到预测值标准误后，乘t分布临界值即可得到置信区间，需要预测区间的话再加残差方差项即可。
• 代码实现步骤
  第一步：提取每个分组模型的关键汇总信息，不需要拆分单个参数

  library(dplyr)
  library(tidyr)
  library(purrr)
  
  # 提取每个Run对应模型的参数、协方差矩阵、残差信息
  model_summary <- model_ex %>%
    mutate(
      # 提取参数估计值
      a_est = map_dbl(model, ~coef(.x)[["a"]]),
      b_est = map_dbl(model, ~coef(.x)[["b"]]),
      # 提取参数协方差矩阵
      vcov_mat = map(model, ~vcov(.x)),
      # 提取残差自由度
      df_resid = map_dbl(model, ~df.residual(.x)),
      # 提取残差标准差（计算预测区间需要）
      sigma_resid = map_dbl(model, ~sigma(.x))
    )
  

  第二步：将待预测的新数据集和模型汇总表按分组合并，全向量化计算所有预测值和区间，处理百万级数据无压力

  # 示例待预测数据集，替换为你自己的实际数据即可
  new_data <- expand.grid(Run = unique(DNAse$Run), conc = seq(0, 10, 0.1))
  
  # 批量计算所有预测值和上下限
  pred_result <- new_data %>%
    left_join(model_summary, by = "Run") %>%
    mutate(
      # 点预测值
      y_pred = a_est * exp(b_est * conc),
      # 计算偏导数
      dyda = exp(b_est * conc),
      dydb = a_est * conc * exp(b_est * conc),
      # 提取协方差矩阵元素
      var_a = map_dbl(vcov_mat, ~.x["a", "a"]),
      var_b = map_dbl(vcov_mat, ~.x["b", "b"]),
      cov_ab = map_dbl(vcov_mat, ~.x["a", "b"]),
      # delta方法计算预测值方差、标准误
      var_pred = dyda^2 * var_a + dydb^2 * var_b + 2*dyda*dydb*cov_ab,
      se_pred = sqrt(var_pred),
      # 95%置信度对应的t临界值
      t_crit = qt(0.975, df = df_resid),
      # 置信区间上下限
      ci_lwr = y_pred - t_crit * se_pred,
      ci_upr = y_pred + t_crit * se_pred,
      # 可选：计算预测区间上下限
      pi_lwr = y_pred - t_crit * sqrt(var_pred + sigma_resid^2),
      pi_upr = y_pred + t_crit * sqrt(var_pred + sigma_resid^2)
    ) %>%
    # 按需删除中间计算列
    select(-c(a_est, b_est, vcov_mat, df_resid, sigma_resid, dyda, dydb, var_a, var_b, cov_ab, var_pred, se_pred, t_crit))
  

• 扩展兼容：如果后续使用更复杂的非线性模型，无法手动推导偏导，可以调用numDeriv::grad()函数数值计算偏导数，依然是向量化处理，效率远高于predictNLS()。

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