Combinando datos de diferentes fuentes

Quiero combinar datos de diferentes fuentes.

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

Tengo algunos datos empíricos, que varían debido al error de medición en torno a la media.

Y, en segundo lugar, tengo un modelo que predice una estimación a partir de otra información (el modelo también tiene cierta incertidumbre).

¿Cómo puedo combinar esos dos conjuntos de datos? [La estimación combinada se utilizará en otro modelo como predictor].

El metanálisis y los métodos bayesianos parecen ser adecuados. Sin embargo, no he encontrado muchas referencias e ideas sobre cómo implementarlo (estoy usando R, pero también estoy familiarizado con Python y C ++).

Gracias.

Actualizar

Ok, aquí hay un ejemplo más real:

Para estimar la toxicidad de un químico (típicamente expresado como$LC_{50}$= concentración donde muere el 50% de los animales) se realizan experimentos de laboratorio. Afortunadamente, los resultados de los experimentos se recopilan en una base de datos (EPA) .

Estos son 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 disponibles para predecir la toxicidad de las propiedades químicas ( QSAR ). Uno de estos modelos predice la toxicidad del coeficiente de partición octanol / agua ($log~K_{OW}$):

El coeficiente de partición de lindano es $log~K_{OW} = 3.8$ y la toxicidad prevista es $log~LC_{50} [mol/L] = -4.902$.

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

¿Hay una buena manera de combinar estas dos informaciones diferentes (experimentos de laboratorio y predicciones de modelos)?

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

El combinado $LC_{50}$se usará más adelante en un modelo como predictor. Por lo tanto, un solo valor (combinado) sería una solución simple.

Sin embargo, una distribución también puede ser útil, si esto es posible en el modelado (¿cómo?).

Su modelo estimado sería un útil previo.

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

Para parametrizar esto antes de usar la simulación, considere su modelo

Asumir $b_0 \sim N(0.94, 0.03)$ y $b_1 \sim N(1.33, 0.1)$; $\textrm{Lkow}$ es conocido (un parámetro fijo; por ejemplo, las constantes físicas a menudo se conocen con mucha precisión en relación con otros parámetros).

Además, hay cierta incertidumbre en el modelo, haré esto $\epsilon \sim N(0,1)$, pero debe ser una representación precisa de su información, por ejemplo, el RMSE del modelo podría usarse para informar la escala de la desviación estándar. Intencionalmente estoy haciendo esto un "informativo" previo.

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

Ahora imagina que thepriores tu anterior y

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

son sus datos:

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

La forma más fácil de usar el anterior será parametrizar una distribución que JAGS reconocerá.

Esto puede hacerse de muchas maneras. Como los datos no tienen que ser normales, puede considerar encontrar una distribución usando el paquete fitdistrplus. Para simplificar, supongamos que su anterior es N(mean(theprior), sd(theprior)), o aproximadamente$N(-4.9, 1.04)$. Si desea inflar la varianza (para dar más fuerza a los datos), puede usar$N(-4.9, 2)$

Entonces podemos ajustar un modelo usando 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)

Finalmente, 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 puede considerar mu=5.08su estimación del valor medio del parámetro (rosa) y sd = 0.8su desviación estándar; la estimación predictiva posterior del logLC_50 (de donde obtiene las muestras) está en rojo.

Referencia

LeBauer, DS, D. Wang, K. Richter, C. Davidson y MC Dietze. (2013) Facilitar retroalimentaciones entre mediciones de campo y modelos de ecosistemas. Monografías ecológicas 83: 133-154. doi: 10.1890 / 12-0137.1