Combiner des données de différentes sources

Je souhaite combiner des données provenant de différentes sources.

Disons que je veux estimer une propriété chimique (par exemple un coefficient de partage ):

J'ai quelques données empiriques, variant en raison d'une erreur de mesure autour de la moyenne.

Et, deuxièmement, j'ai un modèle prédisant une estimation à partir d'autres informations (le modèle présente également une certaine incertitude).

Comment puis-je combiner ces deux jeux de données? [L'estimation combinée sera utilisée dans un autre modèle comme prédicteur].

La méta-analyse et les méthodes bayésiennes semblent convenir. Cependant, je n'ai pas trouvé beaucoup de références et d'idées pour l'implémenter (j'utilise R, mais je suis également familier avec python et C ++).

Merci.

Mise à jour

Ok, voici un exemple plus réel:

Pour estimer la toxicité d'un produit chimique (généralement exprimée en $LC_{50}$= concentration où 50% des animaux meurent), des expériences en laboratoire sont menées. Heureusement, les résultats des expériences sont rassemblés dans une base de données (EPA) .

Voici quelques valeurs pour l'insecticide Lindane :

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

Cependant, il existe également certains modèles pour prédire la toxicité des propriétés chimiques ( QSAR ). L'un de ces modèles prédit la toxicité à partir du coefficient de partage octanol / eau ($log~K_{OW}$):

Le coefficient de partage du lindane est $log~K_{OW} = 3.8$ et la toxicité prévue est $log~LC_{50} [mol/L] = -4.902$.

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

Existe-t-il une belle façon de combiner ces deux informations différentes (expériences de laboratoire et prédictions de modèles)?

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

Le combiné $LC_{50}$sera utilisé plus tard dans un modèle comme prédicteur. Par conséquent, une seule valeur (combinée) serait une solution simple.

Cependant, une distribution peut également être utile - si cela est possible dans la modélisation (comment?).

Votre estimation de modèle serait un préalable utile.

J'ai appliqué l'approche suivante dans LeBauer et al 2013 , et j'ai adapté le code de priors_demo.Rmd ci-dessous.

Pour paramétrer cela avant d'utiliser la simulation, considérez votre modèle

Présumer $b_0 \sim N(0.94, 0.03)$ et $b_1 \sim N(1.33, 0.1)$; $\textrm{Lkow}$ est connu (un paramètre fixe; par exemple, les constantes physiques sont souvent connues très précisément par rapport à d'autres paramètres).

De plus, il y a une certaine incertitude du modèle, je vais en faire $\epsilon \sim N(0,1)$, mais doit être une représentation précise de vos informations, par exemple, le RMSE du modèle peut être utilisé pour informer l'échelle de l'écart-type. J'en fais intentionnellement un avant «informatif».

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

Imaginez maintenant thepriorest votre avant et

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

sont vos données:

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

Le moyen le plus simple d'utiliser l'a priori sera de paramétrer une distribution que JAGS reconnaîtra.

Cela peut se faire de plusieurs manières. Étant donné que les données ne doivent pas être normales, vous pouvez envisager de trouver une distribution à l'aide du package fitdistrplus. Pour simplifier, supposons simplement que votre a priori est N(mean(theprior), sd(theprior)), ou approximativement$N(-4.9, 1.04)$. Si vous voulez gonfler la variance (pour donner plus de force aux données), vous pouvez utiliser$N(-4.9, 2)$

Ensuite, nous pouvons adapter un modèle à l'aide de 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)

Enfin, un complot:

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

Et vous pouvez considérer mu=5.08comme votre estimation de la valeur moyenne du paramètre (rose) et de sd = 0.8son écart-type; l'estimation prédictive postérieure du logLC_50 (d'où vous obtenez vos échantillons) est en rouge.

Référence

LeBauer, DS, D. Wang, K. Richter, C. Davidson et MC Dietze. (2013). Faciliter les rétroactions entre les mesures sur le terrain et les modèles d'écosystème. Monographies écologiques 83: 133-154. doi: 10.1890 / 12-0137.1