机器算法验证 - 结合不同来源的数据 - 吾爱随笔录

结合不同来源的数据

机器算法验证 r 贝叶斯荟萃分析锯齿

2022-03-18 07:41:19

我想合并来自不同来源的数据。

假设我想估计一种化学性质（例如分配系数）：

我有一些经验数据，由于平均值周围的测量误差而有所不同。

其次，我有一个模型可以根据其他信息预测估计值（该模型也有一些不确定性）。

如何结合这两个数据集？[组合估计将在另一个模型中用作预测器]。

Meta 分析和贝叶斯方法似乎是合适的。但是，还没有找到很多参考和想法如何实现它（我使用的是R，但也熟悉python和C++）。

谢谢。

更新

好的，这是一个更真实的例子：

为了估计化学物质的毒性（通常表示为 = 50% 动物死亡的浓度），进行了实验室实验。令人高兴的是，实验结果收集在数据库(EPA)中。 $LC_{50}$

以下是杀虫剂林丹的一些值：

### Toxicity of Lindane in ug/L
epa <- c(850 ,6300 ,6500 ,8000, 1990 ,516, 6442 ,1870, 1870, 2000 ,250 ,62000,
         2600,1000,485,1190,1790,390,1790,750000,1000,800
)
hist(log10(epa))

# or in mol / L
# molecular weight of Lindane
mw = 290.83 # [g/mol]
hist(log10(epa/ (mw * 1000000)))

然而，也有一些模型可用于预测化学性质的毒性 ( QSAR )。其中一种模型根据辛醇/水分配系数 ( ) 预测毒性： $log~K_{OW}$

l o g L C_{50} [m o l / L] = 0.94 (\pm 0.03) l o g K_{O W} - 1.33 (\pm 0.1)

$log~LC_{50} [mol/L] = 0.94~(\pm 0.03)~log~K_{OW}~-~1.33 (\pm~0.1)$

林丹的分配系数为，预测毒性为。 $log~K_{OW} = 3.8$ $log~LC_{50} [mol/L] = -4.902$

lkow = 3.8
mod1 <- -0.94 * lkow - 1.33
mod1

有没有一种很好的方法来结合这两种不同的信息（实验室实验和模型预测）？

hist(log10(epa/ (mw * 1000000)))
abline(v = mod1, col = 'steelblue')

组合的稍后将在模型中用作预测器。因此，单个（组合）值将是一个简单的解决方案。 $LC_{50}$

但是，分布也可能很方便 - 如果这在建模中是可能的（如何？）。

1个回答

您的模型估计将是一个有用的先验。

我在LeBauer et al 2013中应用了以下方法，并从下面的priors_demo.Rmd中改编了代码。

要使用模拟参数化此先验，请考虑您的模型

{logLC}_{50} = b_{0} X + b_{1}

$\textrm{logLC}_{50} = b_0 X+b_1$

假设和；是已知的（一个固定的参数；例如，物理常数通常相对于其他参数非常精确地已知）。 $b_0 \sim N(0.94, 0.03)$ $b_1 \sim N(1.33, 0.1)$ $\textrm{Lkow}$

此外，还有一些模型不确定性，我会做这个，但应该是您信息的准确表示，例如模型的 RMSE 可用于告知标准的规模偏差。我故意将其作为“信息”之前。 $\epsilon \sim N(0,1)$

b0 <- rnorm(1000, -0.94, 0.03)
b1 <- rnorm(1000, -1.33, 0.1)
e <- rnorm(1000, 0, 1)
lkow <- 3.8
theprior <- b0 * lkow + b1 + e

现在想象theprior是你的先前和

thedata <- log10(epa/ (mw * 1000000))

是你的数据：

library(ggplot2)
ggplot() + geom_density(aes(theprior)) + theme_bw() + geom_rug(aes(thedata))

使用先验的最简单方法是参数化 JAGS 将识别的分布。

这可以通过多种方式完成。由于数据不必是正常的，您可以考虑使用 package 查找分布fitdistrplus。为简单起见，让我们假设您的先验是N(mean(theprior), sd(theprior))或大约。如果您想扩大方差（以使数据更具强度），您可以使用 $N(-4.9, 1.04)$ $N(-4.9, 2)$

然后我们可以使用 JAGS 拟合模型

writeLines(con = "mymodel.bug",
           text = "
           model{
             for(k in 1:length(Y)) {
               Y[k] ~ dnorm(mu, tau)
             }

             # informative prior on mu
             mu ~ dnorm(-4.9, 0.25) # precision tau = 1/variance
             # weak prior 
             tau ~ dgamma(0.01, 0.01)
             sd <- 1 / sqrt(tau)
           }")

require(rjags)
j.model  <- jags.model(file = "mymodel.bug", 
                                  data = data.frame(Y = thedata), 
                                  n.adapt = 500, 
                                  n.chains = 4)
mcmc.object <- coda.samples(model = j.model, variable.names = c('mu', 'tau'),
                            n.iter = 10000)
library(ggmcmc)

## look at diagnostics
ggmcmc(ggs(mcmc.object), file = NULL)

## good convergence, but can start half-way through the simulation
mcmc.o     <- window(mcmc.object, start = 10000/2)
summary(mcmc.o)

最后，一个情节：

ggplot() + theme_bw() + xlab("mu") + 
     geom_density(aes(theprior), color = "grey") + 
     geom_rug(aes(thedata)) + 
     geom_density(aes(unlist(mcmc.o[,"mu"])), color = "pink") +
     geom_density(aes(unlist(mcmc.o[,"pred"])), color = "red")

并且您可以考虑mu=5.08是您对平均参数值（粉红色）sd = 0.8及其标准差的估计；logLC_50（您从中获取样本的位置）的后验预测估计为红色。

在此处输入图像描述

参考

LeBauer、DS、D. Wang、K. Richter、C. Davidson 和 MC Dietze。（2013）。促进现场测量和生态系统模型之间的反馈。生态专着八十三：133—154。doi:10.1890/12-0137.1

其它你可能感兴趣的问题

上一篇R混淆中的两样本Kolmogorov-Smirnov检验p值下一篇计算困惑度