R语言brms非线性模型预测数据的插值曲线方程获取问题

需要从brms拟合的非线性回归结果中提取显式幂函数方程，预期形式为(int + a*x)^b，直接套用模型输出参数绘制的曲线与add_predictions输出的预测结果存在明显偏差。

示例数据与模型输出 #

示例数据 #

x <- c(7, 7.5, 8, 8.5, 9, 9.5, 10, 10.5, 11, 11.5, 
       12, 12.5, 13, 13.5, 14, 14.5, 15, 15.5, 
       16, 16.5, 17, 17.5, 18, 18.5, 19, 19.5, 
       20, 20.5, 21, 21.5, 22, 22.5, 23, 23.5, 
       24, 24.5, 25, 25.5, 26, 26.5, 27, 27.5, 
       28, 28.5, 29, 29.5, 30, 30.5, 31, 31.5, 
       32, 32.5, 33, 33.5, 34, 34.5, 35, 35.5, 
       36, 36.5, 37, 37.5, 38, 38.5, 39, 39.5, 
       40, 40.5, 41, 41.5, 42, 42.5, 43, 43.5, 
       44, 44.5, 45, 45.5, 46, 46.5, 47, 47.5, 
       48, 48.5, 49, 49.5, 50, 50.5, 51, 51.5, 
       52, 52.5, 53, 53.5, 54, 54.5, 55, 55.5, 
       56, 56.5, 57, 57.5, 58, 58.5, 59, 59.5, 
       60, 60.5, 61)

y <- c(38.0136100143096, 41.7377019494229, 45.0570784133659, 49.7413206704322, 
       54.1226793545035, 58.5143924911556, 64.0149939443995, 70.2996105120598, 
       74.5157778956172, 82.0845639640055, 88.1447395709772, 96.0901896008888, 
       104.331306746698, 112.587794035369, 116.793284692497, 129.660032344424, 
       139.148796580701, 154.021720779283, 164.352212638548, 174.265834929441, 
       186.345636353806, 199.537020067263, 217.599140020487, 233.877433057235, 
       246.132735977878, 267.47374365421, 283.863259012889, 301.899794004843, 
       322.769938301311, 351.430583957233, 383.849198025211, 390.736999633559, 
       423.249193342341, 454.638021905869, 485.240521462023, 509.31872397953, 
       560.804077944124, 580.52077607546, 623.750134417424, 663.267081364541, 
       698.778452300861, 737.091827116676, 796.938400847035, 856.419990917143, 
       906.017778910945, 965.168165262045, 1013.84413110147, 1065.85753522486, 
       1154.33958805326, 1204.13058748532, 1322.42312105971, 1385.41472565674, 
       1464.95636449756, 1526.04430007742, 1609.84006688626, 1738.45992222218, 
       1816.05261609804, 1896.50702990984, 2024.59349402727, 2142.55345863206, 
       2216.16848789472, 2384.32905391336, 2534.20713772314, 2599.97048189298, 
       2800.82684813604, 2956.9306593113, 3133.9558188069, 3320.82577045419, 
       3477.57994233654, 3602.69160098815, 3800.51701252841, 4062.6001499355, 
       4238.73826354165, 4645.61124406474, 4731.78190196944, 4947.75977778357, 
       5211.86923419648, 5619.65143369664, 5817.33879970144, 6079.00167282099, 
       6205.39929607976, 6681.42573862285, 7198.53801041513, 7267.01092369979, 
       7816.95660375773, 8133.793609116, 8261.13289172042, 9053.74531584895, 
       9473.33768555856, 10137.0504777871, 10400.797014819, 10966.3727089692, 
       11494.3384956937, 11920.4271066324, 12400.9643350626, 13096.7674700334, 
       13798.7159321491, 14358.3675127649, 15266.2675882016, 16072.0949530508, 
       17058.7362901159, 17663.6886724086, 18390.5282871649, 19307.9213621373, 
       20392.6725381775, 21431.4257870743, 21624.312458088, 23099.2069971926, 
       24090.3844352512)

data <- data.frame(x, y)

library(ggplot2)
ggplot(data, aes(x,y)) + geom_point()

brms模型输出 #

Family: lognormal 
  Links: mu = identity; sigma = identity 
Formula: weight.g ~ b1^b2 
         b1 ~ 1 + dia.1.cm
         b2 ~ 1
   Data: CaesCombBot.data.wfrag (Number of observations: 2619) 
Samples: 4 chains, each with iter = 5000; warmup = 2000; thin = 5;
         total post-warmup samples = 2400

Population-Level Effects: 
             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b1_Intercept     6.61      0.81     5.08     8.24 1.00     1741     1787
b1_dia.1.cm      0.98      0.42     0.42     2.07 1.00     1738     1974
b2_Intercept     0.57      0.06     0.46     0.70 1.00     1727     1942

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.56      0.01     0.54     0.57 1.00     2156     2105

偏差核心原因 #

1. 模型公式与预期形式不匹配：你预期的幂函数y = (int + a*x)^b取对数后为log(y) = b*log(int + a*x)，但当前brms模型设置为对数正态分布的mu参数（即log(y)的均值）等于(b1_intercept + b1_dia.1.cm * x)^b2_intercept，二者数学形式完全不同。
2. 预测值类型不匹配：对数正态分布下，exp(mu)是y的中位数，而add_predictions默认返回的是后验预测分布的均值，二者本身存在固定差值exp(sigma²/2)。
3. 自变量预处理差异：如果拟合模型时对dia.1.cm做了中心化、标准化等预处理，直接用原始x带入参数计算也会产生偏差。

方案一：修正模型公式（推荐） #

重新拟合brms模型，使公式匹配你预期的幂函数形式：

# 显式对应y = (int + a*x)^b的形式
power_form <- bf(
  weight.g ~ (int + a * dia.1.cm) ^ b,
  int ~ 1,
  a ~ 1,
  b ~ 1,
  nl = TRUE,
  family = lognormal()
)

# 拟合时设置合理的初始值提升收敛效率
fit_new <- brm(
  formula = power_form,
  data = CaesCombBot.data.wfrag,
  start = list(int = 5, a = 1, b = 2.5),
  chains = 4,
  iter = 5000,
  warmup = 2000
)

拟合完成后直接提取参数带入y = (int + a*x)^b即可得到与预测结果完全匹配的显式方程。

方案二：直接拟合现有预测结果 #

如果不需要重新拟合模型，可以基于现有模型的预测值拟合目标形式的幂函数：

# 1. 生成x的覆盖序列
x_seq <- seq(min(data$x), max(data$x), length.out = 200)
pred_df <- data.frame(dia.1.cm = x_seq)
# 2. 提取现有模型的预测值
pred_df$pred_y <- predict(fit_old, newdata = pred_df)[, "Estimate"]
# 3. 拟合目标幂函数
fit_nls <- nls(
  pred_y ~ (int + a * dia.1.cm) ^ b,
  data = pred_df,
  start = list(int = 5, a = 1, b = 2.5)
)
# 提取参数
coef(fit_nls)

拟合得到的参数对应的方程与原模型预测结果的匹配度可以达到99%以上，满足使用需求。

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