La combinación de datos de diferentes fuentes

Question

Quiero combinar datos de diferentes fuentes.

Digamos que yo quiero para estimar una propiedad química (por ejemplo, un coeficiente de partición):

Tengo algunos datos empíricos, variando debido al error de medición alrededor de la media.

Y, en segundo lugar, tengo un modelo de predicción de una estimación del resto de la información (el modelo tiene también algo de incertidumbre).

¿Cómo puedo combinar los dos conjuntos de datos? [La estimación conjunta será utilizado en otro modelo como predictor].

Meta-análisis bayesiano y métodos parecen ser adecuados. Sin embargo, no he encontrado muchas referencias e ideas de cómo implementar (estoy usando R, pero también familiarizados con python y C++).

Gracias.

Actualización

Ok, he aquí una más real, ejemplo:

Para estimar la toxicidad de un producto químico (normalmente expresado como $LC_{50}$ = concentración donde el 50% de los animales se mueren) laboratorio de experimentos se llevan a cabo. Felizmente los resultados de los experimentos están reunidos en una base de datos (EPA).

Aquí están algunos valores para el insecticida Lindano:

### 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)))

Sin embargo, también hay algunos modelos para predecir la toxicidad de las propiedades químicas (QSAR). Uno de estos modelos predice la toxicidad de la octanol/agua coeficiente de partición ($log~K_{OW}$):

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

The partitioning coefficient of Lindane is $log~K_{UJO} = 3.8$ and the predicted toxicity is $ log~LC_{50} [mol/L] = -4.902 $.

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

Is there a nice way to combine these two different informations (lab experimens and model predictions)?

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

The combined $LC_{50}$ será utilizado más adelante en un modelo como predictor. Por lo tanto, una sola (combinado) valor sería una solución simple.

Sin embargo, una distribución podría ser también práctico - si esto es posible en la modelización (¿cómo?).

Answer 1

Su modelo de estimación sería de gran utilidad antes.

He aplicado el siguiente enfoque en LeBauer et al 2013, y se han adaptado código de priors_demo.Rmd a continuación.

Para parametrizar este antes de utilizar la simulación, considere la posibilidad de su modelo

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

Assume $b_0 \sim N(0.94, 0.03)$ and $b_1 \sim N(1.33, 0.1)$; $\textrm{Lkow}$ is known (a fixed parameter; for example physical constants are often known very precisely relative to other parameters).

In addition, there is some model uncertainty, I'll make this $\epsilon \sim N(0,1)$, but should be an accurate representation of your information, for example the model's RMSE could be used to inform the scale of the standard deviation. I am intentionally making this an 'informative' prior.

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

Now imagine theprior is your prior and

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

is your data:

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

The easiest way to use the prior is going to be to parameterize a distribution that JAGS will recognize.

This can be done in many ways. Since the data don't have to be normal, you might consider finding a distribution using the package fitdistrplus. For simplicity, lets just assume that your prior is N(mean(theprior), sd(theprior)), or approximately $N(-4.9, 1.04)$. If you want to inflate the variance (to give the data more strength) you could use $N(-4.9, 2)$

A continuación, podemos ajustar un modelo de uso de PUNTAS

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)

Por último, una trama:

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")

Y se pueden considerar mu=5.08 a ser de su estimación de la media del valor del parámetro (rosa), y sd = 0.8 su desviación estándar; la posterior predicción de la estimación de la logLC_50 (donde es llegar sus muestras de) está en rojo.

Referencia

LeBauer, D. S., D. Wang, K. Richter, C. Davidson, & M. C. Dietze. (2013). Facilitar la retroalimentación entre las mediciones de campo y los modelos. Ecológica Monografías 83:133-154. doi:10.1890/12-0137.1