如何基于ulam模型绘制wind_sd预测值与实际值对比图？

我是环境科学专业一年级研究生，正在R中开展项目，通过多个气象参数研究阵风的变异性。在为单个输入变量建立模型后，我构建了首个多水平模型，目前仅使用wind_spd（逐小时平均风速）、t_2m（2米高度气温）两个输入变量，后续可能添加更多变量，模型采用ulam函数构建。

导师要求我绘制输出参数wind_sd（逐小时风速标准差，用于衡量阵风变异性）的预测值（模型生成）与实际值（原始数据）对比图。经查找，现有绘图指南均基于lm()线性回归模型，且predict()函数不适用于ulam模型；我尝试使用extract.samples，但仅得到先验模拟值。

我尝试构建lm模型绘图，但该模型无法预测1.5标准差以上的阵风，而我的数据集（过去35年逐小时数据）关注的阈值为5标准差，且存在大量逐小时阵风远高于小时平均风速的情况。

我使用的是正态分布，是否需要改用二项分布？若无需更改，该如何提取输出参数wind_sd的模拟值或完成绘图？

以下是我构建ulam模型及lm绘图的代码：

dat_slim_MR <- list(
  wind_sd = new_combined_cc$wind_sd.x,
  t_2m = new_combined_cc$t_2m,
  wind_spd = new_combined_cc$wind_spd
)
str(dat_slim_MR)
dat_slim_MR_df <- as_tibble(dat_slim_MR)
dat_slim_MR_df <- dat_slim_MR_df %>% drop_na(wind_spd)

#1 chain to debug

mgustyMR <- ulam(
  alist(
    wind_sd ~ dnorm( mu , sigma ) ,
    mu <- a + b1 * wind_spd + b2 * t_2m,
    a ~ dnorm(0.5, 0.5),
    b1 ~ dnorm(0, 0.02),
    b2 ~ dnorm(0, 0.05),
    sigma ~ dexp(1)
  ) , data=dat_slim_MR_df, chains =1 )

#4 chains for accuracy

mgustyMR_4chains <- ulam(
  alist(
    wind_sd ~ dnorm( mu , sigma ) ,
    mu <- a + b1 * wind_spd + b2 * t_2m,
    a ~ dnorm(0.5, 0.5),
    b1 ~ dnorm(0, 0.02),
    b2 ~ dnorm(0, 0.05),
    sigma ~ dexp(1)
  ) , data=dat_slim_MR_df, chains =4, cores=4, iter=1000 )

#tried to extract samples but only gave values for priors...

postmgustyMR4chains <- extract.samples( mgustyMR_4chains ,  clean = FALSE)
str(postmgustyMR4chains)

#this isn't for the model I've made in ulam... it's A model but I want to use the model with the Bayesian priors

my_mod <- lm(wind_sd ~ wind_spd + t_2m, data = dat_slim_MR_df)
plot(x=predict(my_mod), y=dat_slim_MR_df$wind_sd,
     xlab='Predicted Values',
     ylab='Actual Values',
     xlim = c(0,5),
     ylim = c(0,5),
     main='Predicted vs. Actual Values')
abline(a = 0,
       b = 1,
       col = "red",
       lwd = 2)

我希望用ulam模型生成类似lm模型的对比图，但找不到适用函数，请问如何提取mgustyMR_4chains中的wind_sd模拟值并与原始数据集的wind_sd值绘图？

1. 分布选择：无需改用二项分布 #

二项分布适用于分类（0/1）或计数型离散响应变量，而wind_sd是连续型的风速标准差，当前使用的正态分布合理。若数据存在明显右偏（高值阵风占比高），可考虑改用对数正态分布或伽马分布拟合尾部数据，但二项分布完全不适用。

2. 提取ulam模型预测值并绘图 #

ulam属于rethinking包函数，需用link()计算线性预测器mu的后验样本，再通过sim()生成wind_sd的后验预测分布，最后用均值/中位数作为点预测值绘制对比图。

步骤1：获取后验预测样本 #

library(rethinking)

# 提取模型后验样本
post <- extract.samples(mgustyMR_4chains)

# 计算每个数据点的mu后验分布
mu <- link(mgustyMR_4chains, data = dat_slim_MR_df)

# 生成wind_sd的后验预测样本（每个数据点对应所有后验样本的模拟值）
wind_sd_pred <- sim(mgustyMR_4chains, data = dat_slim_MR_df)

步骤2：计算点预测值（以均值为例） #

# 对每个数据点的后验预测样本取均值，得到点预测值
wind_sd_pred_mean <- apply(wind_sd_pred, 2, mean)

步骤3：绘制预测值vs实际值对比图 #

# 绘制基础对比图
plot(x = wind_sd_pred_mean, y = dat_slim_MR_df$wind_sd,
     xlab = 'Predicted Values (ulam模型)',
     ylab = 'Actual Values',
     xlim = c(0, 5),
     ylim = c(0, 5),
     main = 'Predicted vs. Actual Values (Bayesian Model)')
abline(a = 0, b = 1, col = "red", lwd = 2)

# 可选：添加95%后验置信区间
wind_sd_pred_lower <- apply(wind_sd_pred, 2, quantile, 0.025)
wind_sd_pred_upper <- apply(wind_sd_pred, 2, quantile, 0.975)
segments(x0 = wind_sd_pred_lower, y0 = dat_slim_MR_df$wind_sd,
         x1 = wind_sd_pred_upper, y1 = dat_slim_MR_df$wind_sd,
         col = "blue", lwd = 0.5)

说明 #

• link()函数计算模型中mu <- a + b1*wind_spd + b2*t_2m的后验样本
• sim()函数基于模型似然（dnorm(mu, sigma)）生成wind_sd的后验预测样本，这就是你需要的模型模拟值
• 若关注5标准差的阈值表现，可调整xlim和ylim到对应范围，查看模型对极端值的拟合效果

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