Daten aus verschiedenen Quellen kombinieren

8

Ich möchte Daten aus verschiedenen Quellen kombinieren.

Angenommen, ich möchte eine chemische Eigenschaft (z. B. einen Verteilungskoeffizienten ) abschätzen :

Ich habe einige empirische Daten, die aufgrund von Messfehlern um den Mittelwert variieren.

Und zweitens habe ich ein Modell, das eine Schätzung aus anderen Informationen vorhersagt (das Modell weist auch einige Unsicherheiten auf).

Wie kann ich diese beiden Datensätze kombinieren? [Die kombinierte Schätzung wird in einem anderen Modell als Prädiktor verwendet].

Metaanalyse und Bayes'sche Methoden scheinen geeignet zu sein. Ich habe jedoch nicht viele Referenzen und Ideen zur Implementierung gefunden (ich verwende R, bin aber auch mit Python und C ++ vertraut).

Vielen Dank.

Aktualisieren

Ok, hier ist ein realeres Beispiel:

Um die Toxizität einer Chemikalie (typischerweise ausgedrückt als = Konzentration, bei der 50% der Tiere sterben) Laborexperimente durchgeführt. Glücklicherweise werden die Ergebnisse der Experimente in einer Datenbank (EPA) gesammelt .LC50

Hier einige Werte für das Insektizid Lindan :

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

Es gibt jedoch auch einige Modelle zur Vorhersage der Toxizität anhand chemischer Eigenschaften ( QSAR ). Eines dieser Modelle sagt die Toxizität anhand des Octanol / Wasser-Verteilungskoeffizienten ( ) voraus :log KOW

log LC50[mol/L]=0.94 (±0.03) log KOW  1.33(± 0.1)

Der Verteilungskoeffizient von Lindan beträgt und die vorhergesagte Toxizität ist .log KOW=3.8log LC50[mol/L]=4.902

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

Gibt es eine gute Möglichkeit, diese beiden unterschiedlichen Informationen (Laborexperimente und Modellvorhersagen) zu kombinieren?

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

Das kombinierte wird später in einem Modell als Prädiktor verwendet. Daher wäre ein einzelner (kombinierter) Wert eine einfache Lösung.LC50

Eine Distribution kann jedoch auch nützlich sein - wenn dies bei der Modellierung möglich ist (wie?).

EDi
quelle
2
Obwohl andere hier vielleicht genug finden, um darauf zu antworten, sehe ich noch nicht genug Informationen, um eine gut begründete Antwort zu unterstützen. Wäre es möglich, die Daten, die Sie kombinieren möchten, etwas genauer zu beschreiben?
whuber
@whuber: Danke für den Kommentar. Ich habe ein spezifischeres Beispiel hinzugefügt und hoffe, dies verdeutlicht, wonach ich suche.
EDi
Die Klarstellung ist hilfreich - danke. Aber könnten Sie ein paar Worte darüber hinzufügen, was das Ergebnis einer "Kombination" dieser Ergebnisse sein würde? Wäre es ein einzelner ? Eine Reihe von ihnen? Ein Konfidenzintervall für sie? Eine Einschätzung, wie gut die Vorhersage zu funktionieren scheint? Etwas anderes? Unabhängig davon, wie sie kombiniert werden sollen, wird sich das Interesse letztendlich darauf konzentrieren, die -Informationen für Entscheidungen zu verwenden, z. B. für die Regulierung der Herstellung, Verwendung oder Entsorgung von Chemikalien. Wie diese Entscheidungen getroffen werden, hat normalerweise einen (starken) Einfluss auf die richtige Kombinationsmethode. L C 50LC50LC50
whuber
Klingt so, als könnten Sie einen der vorherigen Schätzungsansätze anwenden, die ich hier entwickelt habe , mit Beispielen in diesem priors_demo.Rmd .
David LeBauer
@David. Danke für das Papier - ich werde es mir ansehen.
EDi

Antworten:

5

Ihre Modellschätzung wäre ein nützlicher Vorgänger.

Ich habe den folgenden Ansatz in LeBauer et al. 2013 angewendet und den Code von priors_demo.Rmd unten angepasst .

Berücksichtigen Sie Ihr Modell, um dies vor der Verwendung der Simulation zu parametrisieren

logLC50=b0X+b1

Angenommen, und ; ist bekannt (ein fester Parameter; beispielsweise sind physikalische Konstanten im Vergleich zu anderen Parametern oft sehr genau bekannt).B 1 ~ N ( 1,33 , 0,1 ) Lkowb0N(0.94,0.03)b1N(1.33,0.1)Lkow

Darüber hinaus gibt es einige Modellunsicherheiten. Ich mache dieses , sollte aber eine genaue Darstellung Ihrer Informationen sein, z. B. könnte der RMSE des Modells verwendet werden, um den Maßstab des Standards zu bestimmen Abweichung. Ich mache dies absichtlich zu einem "informativen" Prior.ϵ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

Stellen thepriorSie sich jetzt vor, Ihr Prior und

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

sind Ihre Daten:

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

Der einfachste Weg, den Prior zu verwenden, besteht darin, eine Verteilung zu parametrisieren, die JAGS erkennt.

Dies kann auf viele Arten erfolgen. Da die Daten nicht normal sein müssen, können Sie eine Verteilung mithilfe des Pakets suchen fitdistrplus. wir der Einfachheit halber einfach an, dass Ihr Prior N(mean(theprior), sd(theprior))oder ungefähr . Wenn Sie die Varianz aufblasen möchten (um den Daten mehr Stärke zu verleihen), können SieN ( - 4,9 , 2 )N(4.9,1.04)N(4.9,2)

Dann können wir ein Modell mit JAGS anpassen

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)

Zum Schluss noch eine Handlung:

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

Und Sie können davon ausgehen mu=5.08, dass Sie den mittleren Parameterwert (pink) und sd = 0.8seine Standardabweichung schätzen . Die hintere prädiktive Schätzung des logLC_50 (von dem Sie Ihre Proben beziehen) ist rot.

Geben Sie hier die Bildbeschreibung ein

Referenz

LeBauer, DS, D. Wang, K. Richter, C. Davidson und MC Dietze. (2013). Erleichterung von Rückkopplungen zwischen Feldmessungen und Ökosystemmodellen. Ecological Monographs 83: 133–154. doi: 10.1890 / 12-0137.1

David LeBauer
quelle
Ich hätte in der vorherigen Berechnung -1,33 durch b1 ersetzen sollen, aber ich habe momentan keine Zeit, dies zu beheben. Es wird keinen großen Unterschied machen.
David LeBauer
@EDi danke - bitte zitieren Sie die beiliegende Referenz, wenn Sie es verwenden!
David LeBauer