使用R语言abc包进行近似贝叶斯计算时出现负估计值问题

问题背景 #

使用R语言abc包进行近似贝叶斯计算（ABC），参数化疾病暴发模型的不确定参数以重现观测模式。调用abc函数采用神经网络回归（method="neuralnet"）估计参数后验分布时，出现以下问题：

• 所有参数本应为正值，但输出结果中出现负估计值
• 收到警告：

Warning messages:
1: In density.default(x, weights = weights) :
  Selecting bandwidth *not* using 'weights'

先后尝试两种汇总统计量均出现该问题：

1. 真实技能统计量（TSS，范围-1至1）
2. 模拟TSS与观测TSS（值为1）的均方根误差（RMSE）

输出结果 #

> summary(abc_neuralnet)
Call: 
abc::abc(target = observed_summary_statistics, param = simulated_parameters, 
    sumstat = simulated_summary_statistics, tol = 1, method = "neuralnet")
Data:
 abc.out$adj.values (10000 posterior samples)
Weights:
 abc.out$weights

                               V1         V2         V3         V4         V5         V6         V7         V8         V9        V10        V11        V12        V13        V14        V15        V16
Min.:                     14.0275    -5.5947   475.6561 -2004.7401 -3410.4807    -0.0041    -0.2877    -0.4942    -0.0403    -0.6970    -0.0524    -0.3172    -0.1729    -0.2109     0.4302     0.0319
Weighted 2.5 % Perc.:     21.0431     1.5006   879.1652  2940.7780  -255.8276     0.2331    -0.1156    -0.1222    -0.0179     0.3437     0.0189    -0.1189     0.0248    -0.1586     0.4668     0.0605
Weighted Median:         130.9660    10.4621  5857.8215 18011.2848  2652.1901     0.4609     0.2799     0.5944     0.3044     0.7649     0.5835     0.5960     0.5098     0.5346     0.7495     0.1958
Weighted Mean:           138.2207    10.6428  5880.5932 18186.6218  2798.3743     0.5053     0.2919     0.5823     0.3040     0.7255     0.5612     0.5943     0.5130     0.5410     0.7519     0.1986
Weighted Mode:            41.2750    15.1685  2046.0809  7565.2865   473.4871     0.2850     0.0511     0.6869     0.3914     0.8491     0.9260     0.8935     0.2530     0.4488     0.5668     0.0947
Weighted 97.5 % Perc.:   283.2321    20.3251 10985.0459 34015.6282  6562.8585     0.9261     0.7165     1.2972     0.6245     0.9763     1.0293     1.2934     1.0046     1.2514     1.0480     0.3448
Max.:                    307.9988    20.4851 14114.8565 45242.1735  7852.7231     1.0047     1.1997     1.4287     0.6768     1.0072     1.1691     1.3752     1.0913     1.3801     1.0690     0.4321
                              V17        V18
Min.:                     -0.0444    -0.7655
Weighted 2.5 % Perc.:      0.1018    -0.2778
Weighted Median:           0.5202     0.4043
Weighted Mean:             0.5449     0.4108
Weighted Mode:             0.1790    -0.0694
Weighted 97.5 % Perc.:     1.0907     1.1251
Max.:                      1.2888     2.0197

代码 #

library(Metrics)
library(abc)

data <- read.csv("C:/Users/Test_abc.csv")

## Define the observed summary statistics
## observed_TSS <- 1
observed_summary_statistics <- c(RMSE_TSS_C1 = 0, RMSE_TSS_C2 = 0)

## Compute the root mean squared error (RMSE) between simulated and observed variables
data$RMSE_TSS_C1 <- sapply(data$TSS_C1, FUN = function(x){Metrics::rmse(1, x)})
data$RMSE_TSS_C2 <- sapply(data$TSS_C2, FUN = function(x){Metrics::rmse(1, x)})
## summary(data) 

## Retrieve the simulated summary statistics
simulated_summary_statistics <- data[, c("RMSE_TSS_C1", "RMSE_TSS_C2")]
## summary(simulated_summary_statistics)

## Retrieve the simulated parameters
simulated_parameters <- data[, c("V1", "V2", "V3", "V4", "V5", "V6", "V7", "V8", "V9", "V10", "V11", "V12", "V13", "V14", "V15", "V16", "V17", "V18")]
## summary(simulated_parameters)

## Run the "cv4abc" function
cv_rejection <- abc::cv4abc(param = simulated_parameters, sumstat = simulated_summary_statistics, nval = 5, tols = c(0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1), method = "rejection", transf = "none")

## Run the "abc" function
abc_neuralnet <- abc::abc(target = observed_summary_statistics, param = simulated_parameters, sumstat = simulated_summary_statistics, tol = 1, method = "neuralnet")

原因分析 #

1. 神经网络外推突破先验范围：
   观测汇总统计对应RMSE=0（即TSS=1），属于极端情况。如果模拟数据中没有生成TSS接近1的样本，神经网络预测时会进行外推，突破模拟参数的正值约束，输出负值。

2. 未施加参数非负约束：
   abc包的神经网络回归方法默认不对参数输出做约束，即使模拟参数全为正，模型也可能因拟合误差产生负值。

3. 汇总统计量信息不足/样本代表性差：
   RMSE=0的极端情况在模拟中可能极少出现，导致神经网络无法学习到参数与汇总统计的稳定映射关系，预测结果异常。同时tol=1使用了所有模拟样本，大量离群样本干扰模型拟合。

4. 带宽警告的潜在关联：
   警告提示密度估计不使用权重，反映后验权重分布可能存在异常，通常与样本中缺乏接近观测点的有效样本有关，间接导致模型拟合不稳定。

解决建议 #

• 检查模拟数据覆盖度：
  查看simulated_summary_statistics中RMSE的最小值，确认是否有样本接近0。如果没有，需调整模拟参数的先验范围，或优化模型使其能生成TSS=1的样本。

• 参数变换约束非负性：
  对参数做对数变换，将参数空间映射到实数域，神经网络预测后再指数变换回原空间，保证参数非负：

  # 对数变换参数
  simulated_parameters_log <- log(simulated_parameters)
  # 运行ABC
  abc_neuralnet_log <- abc::abc(target = observed_summary_statistics, param = simulated_parameters_log, sumstat = simulated_summary_statistics, tol = 1, method = "neuralnet")
  # 转换回原参数空间
  abc_neuralnet_log$adj.values <- exp(abc_neuralnet_log$adj.values)
  

• 优化tol参数：
  参考cv4abc的结果选择最优tol值（通常小于1），只保留接近观测汇总统计的样本，提升模型拟合精度。

• 尝试其他ABC方法：
  对比method="rejection"或method="loclinear"的结果，如果其他方法也出现负值，说明是数据或模型设计问题；若仅神经网络出现，可考虑更换方法或添加约束。

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