Menggabungkan data dari berbagai sumber

8

Saya ingin menggabungkan data dari berbagai sumber.

Katakanlah saya ingin memperkirakan properti kimia (misalnya koefisien partisi ):

Saya memiliki beberapa data empiris, bervariasi karena kesalahan pengukuran di sekitar rata-rata.

Dan, kedua, saya memiliki model yang memperkirakan perkiraan dari informasi lain (model juga memiliki beberapa ketidakpastian).

Bagaimana saya bisa menggabungkan kedua set data itu? [Estimasi gabungan akan digunakan dalam model lain sebagai prediktor].

Meta-analisis dan metode bayesian tampaknya cocok. Namun, belum menemukan banyak referensi dan ide bagaimana mengimplementasikannya (saya menggunakan R, tetapi juga akrab dengan python dan C ++).

Terima kasih.

Memperbarui

Oke, ini contoh yang lebih nyata:

Untuk memperkirakan toksisitas bahan kimia (biasanya dinyatakan sebagai L.C50= konsentrasi ketika 50% hewan mati) percobaan laboratorium dilakukan. Untungnya hasil percobaan dikumpulkan dalam database (EPA) .

Berikut adalah beberapa nilai untuk insektisida 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)))

Namun, ada juga beberapa model yang tersedia untuk memprediksi toksisitas dari sifat kimia ( QSAR ). Salah satu model ini memprediksi toksisitas dari koefisien partisi oktanol / air (lHaig KHAIW):

lHaig L.C50[mHail/L.]=0,94 (±0,03) lHaig KHAIW - 1.33(± 0,1)

Koefisien partisi Lindane adalah lHaig KHAIW=3.8 dan toksisitas yang diprediksi adalah lHaig L.C50[mHail/L.]=-4.902.

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

Apakah ada cara yang baik untuk menggabungkan dua informasi yang berbeda ini (percobaan laboratorium dan prediksi model)?

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

Gabungan L.C50akan digunakan nanti dalam model sebagai prediktor. Oleh karena itu, nilai tunggal (gabungan) akan menjadi solusi sederhana.

Namun, distribusi mungkin juga berguna - jika ini dimungkinkan dalam pemodelan (bagaimana?).

EDI
sumber
2
Walaupun orang lain mungkin menemukan cukup di sini untuk menanggapi, saya belum melihat bahwa ada cukup informasi untuk mendukung jawaban yang beralasan. Apakah mungkin untuk sedikit lebih spesifik tentang data yang Anda rencanakan untuk digabungkan?
whuber
@whuber: Terima kasih atas komentarnya. Saya menambahkan contoh yang lebih spesifik dan melompat ini menjelaskan apa yang saya cari.
EDi
Klarifikasi ini membantu - terima kasih. Tetapi bisakah Anda menambahkan beberapa kata tentang apa hasil "kombinasi" dari hasil ini? Apakah itu tunggalL.C50? Berbagai dari mereka? Interval kepercayaan untuk mereka? Penilaian seberapa bagus prediksi tersebut berfungsi? Sesuatu yang lain Dan, terlepas dari bagaimana mereka akan digabungkan, pada akhirnya minat akan fokus pada penggunaanL.C50informasi untuk membuat keputusan, seperti mengatur pembuatan, penggunaan, atau pembuangan bahan kimia. Bagaimana keputusan ini dibuat biasanya memiliki pengaruh (kuat) pada metode kombinasi yang tepat untuk digunakan.
whuber
Sepertinya Anda bisa menerapkan salah satu pendekatan estimasi sebelumnya yang saya kembangkan di sini , dengan contoh-contoh dalam priors_demo.Rmd ini .
David LeBauer
@ David. Terima kasih untuk kertasnya - Saya akan melihat.
EDi

Jawaban:

5

Perkiraan model Anda akan menjadi prior yang bermanfaat.

Saya telah menerapkan pendekatan berikut dalam LeBauer et al 2013 , dan telah mengadaptasi kode dari priors_demo.Rmd di bawah ini.

Untuk parameter ini sebelum menggunakan simulasi, pertimbangkan model Anda

logLC50=b0X+b1

Menganggap b0N(0,94,0,03) dan b1N(1.33,0,1); Lkow diketahui (parameter tetap; misalnya konstanta fisik sering dikenal sangat relatif terhadap parameter lainnya).

Selain itu, ada beberapa ketidakpastian model, saya akan membuatnya ϵN(0,1), tetapi harus merupakan representasi akurat dari informasi Anda, misalnya RMSE model dapat digunakan untuk menginformasikan skala deviasi standar. Saya sengaja menjadikan ini sebagai 'informatif' sebelumnya.

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

Sekarang bayangkan theprioradalah prioritas Anda

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

adalah data Anda:

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

Cara termudah untuk menggunakan prior adalah dengan parameterisasi distribusi yang akan dikenali JAGS.

Ini bisa dilakukan dengan banyak cara. Karena data tidak harus normal, Anda mungkin mempertimbangkan untuk menemukan distribusi menggunakan paket fitdistrplus. Untuk mempermudah, anggap saja bahwa prior Anda adalah N(mean(theprior), sd(theprior)), atau kurang-lebihN(-4.9,1.04). Jika Anda ingin mengembang varians (untuk memberikan data lebih banyak kekuatan) yang bisa Anda gunakanN(-4.9,2)

Lalu kita bisa muat model menggunakan 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)

Akhirnya, sebuah 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")

Dan Anda dapat mempertimbangkan mu=5.08untuk menjadi estimasi nilai parameter rata-rata (merah muda), dan sd = 0.8standar deviasinya; estimasi prediktif posterior dari logLC_50 (di mana Anda mendapatkan sampel Anda) berwarna merah.

masukkan deskripsi gambar di sini

Referensi

LeBauer, DS, D. Wang, K. Richter, C. Davidson, & MC Dietze. (2013). Memfasilitasi umpan balik antara pengukuran lapangan dan model ekosistem. Monografi Ekologis 83: 133–154. doi: 10.1890 / 12-0137.1

David LeBauer
sumber
Saya seharusnya mengganti -1,33 dengan b1 pada perhitungan sebelumnya, tapi saya tidak punya waktu untuk memperbaikinya sekarang. Itu tidak akan membuat banyak perbedaan.
David LeBauer
@EDi terima kasih - sebutkan referensi yang disertakan jika Anda menggunakannya!
David LeBauer