贝叶斯核机器回归拓展R包：bkmrhat

`1.摘要`

bkmrhat包是用于扩展bkmr包的贝叶斯核机器回归（Bayesian Kernel Machine Regression, BKMR）分析工具，支持多链推断和诊断。该包利用future, rstan, 和coda包的功能，提供了在贝叶斯半参数广义线性模型下进行identity链接和 probit 链接的方法。

主要功能包括：多链合并、继续采样、诊断和预测等。包内包含多种函数，如kmbayes_parallel用于并行计算多个链，kmbayes_combine和kmbayes_combine_lowmem用于合并链，as.mcmc.bkmrfit将bkmrfit对象转换为MCMC对象以进行诊断，以及predict.bkmrfit用于生成预测。

2. 函数介绍

2.1 包介绍

提供了扩展bkmr包的贝叶斯核机器回归工具
支持多链推断和诊断
利用future, rstan, coda包

2.2 as.mcmc.bkmrfit

函数书写格式：as.mcmc()
将bkmrfit对象转换为coda包的MCMC对象
coda 包支持许多不同类型的单链 MCMC 诊断，包括 geweke.diag、traceplot 和 effectiveSize。还可以使用后总结，例如 HPDinterval 和summary.mcmc。
用于进行单链MCMC诊断和后验概括

示例代码1

# 加载bkmrhat包
library(bkmrhat)

# 例子
set.seed(111) #设置随机数种子
library(coda) #加载coda包
# 加载bkmr包
library(bkmr) 
# 生成模拟数据
dat <- bkmr::SimData(n = 50, M = 4)
# 提取数据
y <- dat$y
Z <- dat$Z
X <- dat$X

set.seed(111)
#   运行模型
fitkm <- kmbayes(y = y, Z = Z, X = X, iter = 500, verbose = FALSE,
                 varsel = FALSE)
# 应用as.mcmc函数
mcmcobj <- as.mcmc(fitkm, iterstart=251)    
# 从bkmr对象中提取MCMC链,模型参数的后验总结
summary(mcmcobj) 
# 与bkmr包中的默认值进行比较，该默认值省略了链的前1/2
summary(fitkm)

2.3 as.mcmc.list.bkmrfit.list

函数书写格式：as.mcmc.list()
转换多链bkmrfit对象为coda的mcmc.list对象，以进行 coda MCMC 诊断
coda 包支持许多不同类型的 MCMC 诊断，包括 geweke.diag、traceplot 和 effectiveSize。还可以使用后总结，例如 HPDinterval 和summary.mcmc。对于某些 MCMC 诊断，例如 gelman.diag 和 gelman.plot，需要使用多个链
适用于多链MCMC诊断

示例代码2

# 运行 2 个并行马尔可夫链（通常更好）
future::plan(strategy = future::multisession, workers=2)

# 使用kmbayes_parallel函数运行马尔可夫链
fitkm.list <- kmbayes_parallel(nchains=2, y = y, Z = Z, X = X, iter = 1000, verbose = FALSE, varsel = FALSE)

# 将结果转换为mcmc对象
mcmcobj = as.mcmc.list(fitkm.list)

# 打印马尔可夫链的摘要统计信息
summary(mcmcobj)

2.4 ExtractPIPs_parallel

函数书写格式：ExtractPIPs()
计算每个链的后验包含概率
“Posterior inclusion probabilities” 🔤后验包含概率🔤

示例代码3

# 设置并行计算策略为多会话（multisession），使用4个工作进程
future::plan(strategy = future::multisession, workers=2)

# 使用kmbayes_parallel函数运行并行马尔可夫链，生成马尔可夫链结果的列表fitkm.list
fitkm.list <- kmbayes_parallel(nchains=2, y = y, Z = Z, X = X, iter = 500, verbose = FALSE, varsel = TRUE)

# 将所有马尔可夫链的结果合并
bigkm = kmbayes_combine(fitkm.list, excludeburnin=FALSE)

# 从合并的马尔可夫链结果中提取参数估计
ests = ExtractEsts(bigkm)  

# 从合并的马尔可夫链结果中提取变量选择的后验概率
ExtractPIPs(bigkm)

2.5 kmbayes_combine

函数书写格式：kmbayes_combine()
合并多个bkmr链
组合包含 BKMR 的多个链适合不同的起始值
可设置自定义燃烧期和是否排除燃烧期

MCMC实施中的若干术语：

贝叶斯统计——6. 贝叶斯统计计算方法-CSDN博客https://blog.csdn.net/weixin_52505631/article/details/136207793

参数名	参数类型	中文解释
fitkm.list	output	bkmrfit 对象的列表，每个对象代表一个 MCMC 链的后验
burnin	numeric	自定义的燃烧数（每条链的燃烧迭代数）。如果为 NULL，则默认为每条链的一半
excludeburnin	logical	是否从最终链中排除燃烧迭代次数？注意，所有bkmr包函数会自动从计算中排除燃烧迭代次数
reorder	logical	确保合并链的前半部分只包含每个单独链的前半部分 - 这允许使用bkmr包中的标准函数，这些函数会自动修剪迭代的前半部分。这可用于后验总结，但某些诊断可能效果不好（自相关性，有效样本量），因此应在单独链上进行诊断检测

示例代码见示例代码3

2.6 kmbayes_combine_lowmem

类似于kmbayes_combine，但降低内存需求
在较低的内存设置中组合多个 BKMR 链
通过部分写入磁盘避免内存不足。此函数将一些结果写入磁盘，而不是尝试在内存中完全处理，这在某些情况下将避免 kmbayes_combine 可能发生的“内存不足【"out of memory"】”错误

示例代码4

# 设置并行计算策略为多会话（multisession），使用2个工作进程
future::plan(strategy = future::multisession, workers=2)

# 使用kmbayes_parallel函数运行并行马尔可夫链，生成马尔可夫链结果的列表fitkm.list
fitkm.list <- kmbayes_parallel(nchains=2, y = y, Z = Z, X = X, iter = 500, verbose = FALSE, varsel = TRUE)

# 将所有马尔可夫链的结果合并，使用kmbayes_combine_lowmem函数，保留全部样本
bigkm = kmbayes_combine_lowmem(fitkm.list, excludeburnin=FALSE)

# 从合并的马尔可夫链结果中提取参数估计
ests = ExtractEsts(bigkm)  # 默认保留样本后半部分

# 从合并的马尔可夫链结果中提取变量选择的后验概率
ExtractPIPs(bigkm)

2.7 kmbayes_continue

继续现有bkmr拟合的采样
不完全从先验开始，而是从最后的参数值开始
使用场景：当您使用 kmbayes 函数进行 MCMC 采样，但您没有获取足够的样本且不想重新开始时，请使用此选项

示例代码5

# 调用 bkmr::kmbayes 函数进行贝叶斯分析，生成初始模型 fitty1
fitty1 = bkmr::kmbayes(y = y, Z = Z, X = X, est.h = TRUE, iter = 100)

# 进行一些诊断分析，以判断100次迭代是否足够（默认设置的迭代次数）
# 添加100个额外的迭代（仅作为示例，仍然不足够）
fitty2 = kmbayes_continue(fitty1, iter = 100)

# 将 fitty2 转换为 mcmc 对象
cobj = as.mcmc(fitty2)

# 输出 mcmc 对象中的变量名称
varnames(cobj)

2.8 kmbayes_diagnose

kmbayes_diag
使用rstan包进行MCMC诊断
报告R-hat等指标：使用 Rhat、ess_bulk 和 ess_tail 函数从 rstan 包中为 MCMC 提供诊断。请注意，仅针对 kmbayes_parallel 中的 bkmrfit.list 对象报告 r-hat

示例代码6

# 使用 kmbayes_parallel 函数进行贝叶斯分析，并创建 fitkm.list 列表对象
# nchains=2 表示使用2个链，y、Z、X 是输入的数据，iter=1000 表示迭代次数为1000，verbose=FALSE 表示关闭冗长的输出，varsel=TRUE 表示进行变量选择
fitkm.list <- kmbayes_parallel(nchains = 2, y = y, Z = Z, X = X, iter = 1000, verbose = FALSE, varsel = TRUE)

# 运行 kmbayes_diag 函数对 fitkm.list 进行诊断分析
kmbayes_diag(fitkm.list)

# 运行 kmbayes_diag 函数对第一个链（fitkm.list[[1]]）进行诊断分析
kmbayes_diag(fitkm.list[[1]])

# 关闭所有连接
closeAllConnections()

2.9 kmbayes_parallel

并行运行多个bkmr链
利用future包加速计算：从 kmbayes 函数拟合平行链。这些链利用future包中的并行处理，可以加速拟合并实现依赖于分散初始值的多个马尔可夫链的诊断。

示例代码见示例代码7

2.10 kmbayes_parallel_continue

继续kmbayes_parallel拟合的采样
使用场景：当您使用 kmbayes_parallel 函数进行 MCMC 采样，但您没有获取足够的样本且不想重新开始时，请使用此选项。
返回多链bkmrfit对象

示例代码7

# 并行计算策略，同时指定2个工作进程
future::plan(strategy = future::multisession, workers = 2)

# 使用 kmbayes_parallel 函数进行贝叶斯分析，创建 fitty1p 对象
fitty1p = kmbayes_parallel(nchains = 2, y = y, Z = Z, X = X)

# 使用 kmbayes_parallel_continue 函数继续在 fitty1p 上进行贝叶斯分析，迭代次数设置为3000，创建 fitty2p 对象
fitty2p = kmbayes_parallel_continue(fitty1p, iter = 3000)

# 将 fitty2p 转换为 mcmc.list 格式，创建 cobj 对象
cobj = as.mcmc.list(fitty2p)

# 绘制 MCMC 对象 cobj 的图形
plot(cobj)

2.11 predict.bkmrfit

函数书写格式：predict()
生成基于后验均值或标准差的预测值
适合与SuperLearner等集成：提供基于后验均值的观察水平预测，或者生成观察预测的后验标准差。此函数对于与仅使用点估计的 SuperLearner 等集成机器学习包进行交互时非常有用。

示例代码8

# 加载bkmr库
library(bkmr)

# 设置种子以确保结果的可复现性
set.seed(111)

# 生成模拟数据
dat <- bkmr::SimData(n = 50, M = 4)
y <- dat$y # 响应变量
Z <- dat$Z # 块变量
X <- dat$X # 协变量

# 再次设置种子以确保结果的可复现性
set.seed(111)

# 拟合贝叶斯知识迁移回归模型
fitkm <- kmbayes(y = y, Z = Z, X = X, iter = 200, verbose = FALSE, varsel = TRUE)

# 预测后验均值
postmean = predict(fitkm)

# 使用Z的一半值进行预测，得到预测结果的均值差异
postmean2 = predict(fitkm, Znew = Z/2)

# 计算后验均值的均值差异
mean(postmean - postmean2)