8
$\begingroup$

I want to combine data from different sources.

Let's say I want to estimate a chemical property (e.g. a partitioning coefficient):

I have some empirical data, varying due to measurement error around the mean.

And, secondly, I have a model predicting an estimate from other information (the model has also some uncertainty).

How can I combine those two datasets? [The combined estimate will be used in another model as predictor].

Meta-analysis and bayesian methods seem to be suitable. However, haven't found many references and ideas how to implement it (I am using R, but also familiar with python and C++).

Thanks.

Update

Ok, here's a more real example:

To estimate the toxicity of a chemical (typically expressed as $LC_{50}$ = concentration where 50% of animals die) lab-experiments are conducted. Happily the results of the experiments are gathered in a database (EPA).

Here are some values for the 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)))

However, there are also some models available to predict toxicity from chemical properties (QSAR). One of these models predicts toxicity from the octanol/water partition coefficient ($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_{OW} = 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}$ will be used later on in a model as predictor. Therefore, a single (combined) value would be a simple solution.

However, a distribution might be also handy - if this is possible in modelling (how?).

$\endgroup$
  • 2
    $\begingroup$ Although others might find enough here to respond to, I do not yet see that there is enough information to support a well-reasoned answer. Would it be possible to be a little more specific about the data you plan to combine? $\endgroup$ – whuber Jul 8 '14 at 19:01
  • $\begingroup$ @whuber: Thanks for the comment. I added a more specific example and hop this clarifies what I am looking for. $\endgroup$ – EDi Jul 9 '14 at 8:59
  • $\begingroup$ The clarification is helpful--thank you. But could you add a few words about what the result of a "combination" of these results would be? Would it be a single $LC_{50}$? A range of them? A confidence interval for them? An assessment of how well the prediction seems to work? Something else? And, regardless of how they would be combined, ultimately interest will focus on using the $LC_{50}$ information for making decisions, such as regulating the manufacture, use, or disposal of chemicals. How these decisions are made usually has a (strong) bearing on the proper method of combination to use. $\endgroup$ – whuber Jul 9 '14 at 13:23
  • $\begingroup$ Sounds like you could apply one of the prior estimation approaches I developed here, with examples in this priors_demo.Rmd. $\endgroup$ – David LeBauer Jul 9 '14 at 18:12
  • $\begingroup$ @David. Thanks for the paper - I'll take a look at. $\endgroup$ – EDi Jul 9 '14 at 18:51
5
$\begingroup$

Your model estimate would be a useful prior.

I have applied the following approach in LeBauer et al 2013, and have adapted code from priors_demo.Rmd below.

To parameterize this prior using simulation, consider your model

$$ \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)$

Then we can fit a model using 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)

Finally, a plot:

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

And you can consider mu=5.08 to be your estimate of the mean parameter value (pink), and sd = 0.8 its standard deviation; the posterior predictive estimate of the logLC_50 (where you are getting your samples from) is in red.

enter image description here

Reference

LeBauer, D.S., D. Wang, K. Richter, C. Davidson, & M.C. Dietze. (2013). Facilitating feedbacks between field measurements and ecosystem models. Ecological Monographs 83:133–154. doi:10.1890/12-0137.1

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ I should have replaced -1.33 with b1 in the prior calculation, but I don't have time to fix it right now. It won't make much difference. $\endgroup$ – David LeBauer Jul 10 '14 at 6:29
  • $\begingroup$ @EDi thanks - please cite the included reference if you use it! $\endgroup$ – David LeBauer Jul 10 '14 at 15:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.