Analisis Bayesian hierarkis tentang perbedaan proporsi

Kenapa Hierarkis? : Saya sudah mencoba meneliti masalah ini, dan dari apa yang saya pahami, ini adalah masalah "hierarkis", karena Anda membuat pengamatan tentang pengamatan dari suatu populasi, daripada membuat pengamatan langsung dari populasi itu. Referensi: http://www.econ.umn.edu/~bajari/iosp07/rossi1.pdf

Kenapa Bayesian? : Juga, saya telah menandainya sebagai Bayesian karena solusi asimptotik / sering mungkin ada untuk "desain eksperimental" di mana setiap "sel" ditugaskan pengamatan yang cukup, tetapi untuk tujuan praktis set data dunia nyata / non-eksperimental (atau pada paling tidak milikku) jarang penduduknya. Ada banyak data agregat, tetapi sel-sel individual mungkin kosong atau hanya memiliki beberapa pengamatan.

Model dalam abstrak:

Biarkan U menjadi populasi unit masing-masing darinya kita dapat menerapkan suatu perlakuan, , dari atau , dan dari masing-masing kita amati pengamatan pertama 1 atau 0 alias sukses dan gagal. Misalkan untuk menjadi probabilitas bahwa pengamatan dari objek yang dalam menghasilkan keberhasilan. Perhatikan bahwa dapat dikorelasikan dengan . ${u_1, u_2, u_3 ... u_N}$ $T$ $A$ $B$ $p_{iT}$ $i \in \{1...N\}$ $i$ $T$ $p_{iA}$ $p_{iB}$

Untuk membuat analisis ini layak, kami (a) mengasumsikan bahwa distribusi $p_A$ dan $p_B$ masing-masing adalah contoh dari keluarga distribusi tertentu, seperti distribusi beta, (b) dan memilih beberapa distribusi sebelumnya untuk hyperparameters.

Contoh model

Saya memiliki tas Magic 8 Balls yang sangat besar. Setiap 8-bola, ketika diguncang, dapat mengungkapkan "Ya", atau "Tidak". Juga, saya bisa mengguncang mereka dengan bola yang menghadap ke atas, atau terbalik (dengan asumsi Magic 8 Balls kami bekerja terbalik ...). Orientasi bola dapat sepenuhnya mengubah kemungkinan menghasilkan "Ya" atau "Tidak" (dengan kata lain, awalnya Anda tidak memiliki keyakinan bahwa $p_{iA}$ berkorelasi dengan $p_{iB}$ ).

Pertanyaan:

Seseorang telah secara acak sampel angka, , unit dari populasi, dan untuk setiap unit telah diambil dan mencatat jumlah sewenang-wenang pengamatan di bawah perawatan dan jumlah sewenang-wenang pengamatan di bawah pengobatan . (Dalam praktiknya, dalam kumpulan data kami, sebagian besar unit hanya akan memiliki pengamatan di bawah satu perawatan) $n$ $A$ $B$

Dengan data ini, saya perlu menjawab pertanyaan-pertanyaan berikut:

Jika saya mengambil unit baru secara acak dari populasi, bagaimana saya bisa menghitung (secara analitik atau stokastik) distribusi posterior gabungan dan ? (Terutama agar kita dapat menentukan perbedaan dalam proporsi yang diharapkan, ) $u_x$ $p_{xA}$ $p_{xB}$ $\Delta=p_{xA}-p_{xB}$
Untuk unit tertentu , , dengan pengamatan keberhasilan dan kegagalan, bagaimana saya bisa menghitung (analitis atau stokastik) distribusi posterior bersama untuk dan , sekali lagi untuk membangun distribusi dari perbedaan proporsi $u_y$ $y \in \{1,2,3...,n\}$ $s_y$ $f_y$ $p_{yA}$ $p_{yB}$ $\Delta_y$ $p_{yA}-p_{yB}$

Pertanyaan bonus: Meskipun kami benar-benar berharap dan sangat berkorelasi, kami tidak secara eksplisit memodelkan itu. Dalam kasus kemungkinan solusi stokastik, saya percaya ini akan menyebabkan beberapa sampler, termasuk Gibbs, menjadi kurang efektif dalam mengeksplorasi distribusi posterior. Apakah ini masalahnya, dan jika demikian, haruskah kita menggunakan sampler yang berbeda, entah bagaimana memodelkan korelasinya sebagai variabel terpisah dan mengubah distribusi untuk membuatnya tidak berkorelasi, atau hanya menjalankan sampler lebih lama? $p_A$ $p_B$ $p$

Jawab kriteria

Saya mencari jawaban yang:

Memiliki kode, menggunakan lebih disukai Python / PyMC, atau melarang itu, JAGS, yang dapat saya jalankan
Mampu menangani input ribuan unit
Dengan cukup unit dan sampel, dapat menghasilkan distribusi untuk , , dan sebagai jawaban untuk Pertanyaan 1, yang dapat ditunjukkan sesuai dengan distribusi populasi yang mendasarinya (akan diperiksa terhadap lembar excel yang disediakan dalam "dataset tantangan" bagian) $p_A$ $p_B$ $\Delta$
Dengan cukup unit dan sampel, dapat menghasilkan distribusi yang tepat untuk , , dan (saya akan menggunakan lembar excel yang disediakan di bagian "dataset tantangan" untuk memeriksa) sebagai jawaban untuk Pertanyaan 2, dan menyediakan beberapa alasan mengapa distribusi ini benar $p_A$ $p_B$ $\Delta$
~~Jika jawabannya mirip dengan model JAGS terakhir yang saya posting, jelaskan mengapa itu bekerja dengan dpar(0.5,1)prior tetapi tidak dengan dgamma(2,20)prior.~~ Terima kasih kepada Martyn Plummer di forum JAGS untuk menangkap kesalahan dalam model JAGS saya. Dalam mencoba menetapkan prior Gamma (Shape = 2, Scale = 20), saya menelepon dgamma(2,20)yang sebenarnya menetapkan prior Gamma (Shape = 2, InverseScale = 20) = Gamma (Shape = 2, Scale = 0,05).

Dataset tantangan

Saya telah membuat beberapa dataset sampel di Excel, dengan beberapa skenario yang berbeda, mengubah ketatnya distribusi p, korelasi di antara mereka, dan membuatnya mudah untuk mengubah input lain. https://docs.google.com/file/d/0B_rPBjs4Cp0zLVBybU1nVnd0ZFU/edit?usp=sharing (~ 8Mb)

Visualisasi dari 4 set data yang disertakan dalam file Excel

Solusi sebagian / upaya saya hingga saat ini

1) Saya mengunduh dan menginstal Python 2.7 dan PyMC 2.2. Awalnya, saya mendapatkan model yang salah untuk dijalankan, tetapi ketika saya mencoba untuk merumuskan kembali model, ekstensi membeku. Dengan menambahkan / menghapus kode, saya telah menentukan kode yang memicu pembekuan adalah mc.Binomial (...), meskipun fungsi ini berhasil dalam model pertama, jadi saya berasumsi ada sesuatu yang salah dengan bagaimana saya telah menentukan model.

import pymc as mc
import numpy as np
import scipy.stats as stats
from __future__ import division
cases=[0,0]
for case in range(2):
    if case==0:
        # Taken from the sample datasets excel sheet, Focused Correlated p's
        c_A_arr, n_A_arr, c_B_arr, n_B_arr=np.loadtxt("data/TightCorr.tsv", unpack=True)

    if case==1:
        # Taken from the sample datasets excel sheet, Focused Uncorrelated p's
        c_A_arr, n_A_arr, c_B_arr, n_B_arr=np.loadtxt("data/TightUncorr.tsv", unpack=True)

    scale=20.0
    alpha_A=mc.Uniform("alpha_A", 1,scale)
    beta_A=mc.Uniform("beta_A", 1,scale)
    alpha_B=mc.Uniform("alpha_B", 1,scale)
    beta_B=mc.Uniform("beta_B", 1,scale)
    p_A=mc.Beta("p_A",alpha=alpha_A,beta=beta_A)
    p_B=mc.Beta("p_B",alpha=alpha_B,beta=beta_B)

    @mc.deterministic
    def delta(p_A=p_A,p_B=p_B):
        return p_A-p_B

    obs_n_A=mc.DiscreteUniform("obs_n_A",lower=0,upper=20,observed=True, value=n_A_arr)
    obs_n_B=mc.DiscreteUniform("obs_n_B",lower=0,upper=20,observed=True, value=n_B_arr)

    obs_c_A=mc.Binomial("obs_c_A",n=obs_n_A,p=p_A, observed=True, value=c_A_arr)
    obs_c_B=mc.Binomial("obs_c_B",n=obs_n_B,p=p_B, observed=True, value=c_B_arr)


    model = mc.Model([alpha_A,beta_A,alpha_B,beta_B,p_A,p_B,delta,obs_n_A,obs_n_B,obs_c_A,obs_c_B])
    cases[case] = mc.MCMC(model)
    cases[case].sample(24000, 12000, 2)

    lift_samples = cases[case].trace('delta')[:]

    ax = plt.subplot(211+case)
    figsize(12.5,5)
    plt.title("Posterior distributions of lift from 0 to T")
    plt.hist(lift_samples, histtype='stepfilled', bins=30, alpha=0.8,
             label="posterior of lift", color="#7A68A6", normed=True)
    plt.vlines(0, 0, 4, color="k", linestyles="--", lw=1)
    plt.xlim([-1, 1])

2) Saya mengunduh dan menginstal JAGS 3.4. Setelah mendapatkan koreksi pada prior saya dari forum JAGS, saya sekarang memiliki model ini, yang berhasil berjalan:

Model

var alpha_A, beta_A, alpha_B, beta_B, p_A[N], p_B[N], delta[N], n_A[N], n_B[N], c_A[N], c_B[N];
model {
    for (i in 1:N) {
        c_A[i] ~ dbin(p_A[i],n_A[i])
        c_B[i] ~ dbin(p_B[i],n_B[i])
        p_A[i] ~ dbeta(alpha_A,beta_A)
        p_B[i] ~ dbeta(alpha_B,beta_B)
        delta[i] <- p_A[i]-p_B[i]
    }
    alpha_A ~ dgamma(1,0.05)
    alpha_B ~ dgamma(1,0.05)
    beta_A ~ dgamma(1,0.05)
    beta_B ~ dgamma(1,0.05)
}

Data

"N" <- 60
"c_A" <- structure(c(0,6,0,3,0,8,0,4,0,6,1,5,0,5,0,7,0,3,0,7,0,4,0,5,0,4,0,5,0,4,0,2,0,4,0,5,0,8,2,7,0,6,0,3,0,3,0,8,0,4,0,4,2,6,0,7,0,3,0,1))
"c_B" <- structure(c(5,0,2,2,2,0,2,0,2,0,0,0,5,0,4,0,3,1,2,0,2,0,2,0,0,0,3,0,6,0,4,1,5,0,2,0,6,0,1,0,2,0,4,0,4,1,1,0,3,0,5,0,0,0,5,0,2,0,7,1))
"n_A" <- structure(c(0,9,0,3,0,9,0,9,0,9,3,9,0,9,0,9,0,3,0,9,0,9,0,9,3,9,0,9,0,9,0,3,0,9,0,9,0,9,3,9,0,9,0,9,0,3,0,9,0,9,0,9,3,9,0,9,0,9,0,3))
"n_B" <- structure(c(9,0,9,3,9,0,9,0,9,0,3,0,9,0,9,0,9,3,9,0,9,0,9,0,3,0,9,0,9,0,9,3,9,0,9,0,9,0,3,0,9,0,9,0,9,3,9,0,9,0,9,0,3,0,9,0,9,0,9,3))

Kontrol

model in Try1.bug
data in Try1.r
compile, nchains(2)
initialize
update 400
monitor set p_A, thin(3)
monitor set p_B, thin(3)
monitor set delta, thin(3)
update 1000
coda *, stem(Try1)

Aplikasi sebenarnya untuk siapa saja yang lebih suka memilih model :)

Di web, pengujian A / B tipikal mempertimbangkan dampak terhadap tingkat konversi dari satu halaman atau unit konten, dengan kemungkinan variasi. Solusi khas meliputi tes signifikansi klasik terhadap hipotesis nol atau dua proporsi yang sama, atau solusi Bayesian analitik yang lebih baru memanfaatkan distribusi beta sebagai konjugat sebelumnya.

Alih-alih pendekatan satu-unit-konten ini, yang kebetulan membutuhkan banyak pengunjung untuk setiap unit yang saya minati, kami ingin membandingkan variasi dalam proses yang menghasilkan banyak unit konten (bukan skenario yang tidak biasa benar-benar ...). Jadi secara keseluruhan, unit / halaman yang dihasilkan oleh proses A atau B memiliki banyak kunjungan / data, tetapi masing-masing unit hanya memiliki beberapa pengamatan.

bayesian mcmc jags hierarchical-bayesian pymc Fabio Beltramini
sumber

Anda telah melewatkan aspek penting dalam deskripsi Anda. Setelah Anda menggambar secara seragam , apa yang terjadi selanjutnya? Sesuatu harus, karena jika tidak, Anda sama sekali tidak memiliki informasi tentang salah satu dari ! Lalu, apa yang Anda amati? Apakah Anda mengamati dan, katakanlah, undian independen dari ? Atau apakah Anda hanya mengamati undian dari dan tidak merekam ? Untuk memperkirakan semua distribusi kepadatan ini, tentunya Anda harus mengulangi prosesnya berkali-kali - berapa banyak? Atau, dengan "memperkirakan" Anda menganggap Anda tahu semua dan Anda ingin menghitung distribusi hasil percobaan Anda?

x

$x$

p_{x}

$p_x$

x

$x$

o_{x}

$o_x$

o_{x}

$o_x$

x

$x$

p_{x}

$p_x$

whuber

Maaf, mungkin frasa " diambil secara seragam dari 1 ... N" menyesatkan. Mengingat sejumlah besar pengamatan yang ada tentang berbagai , saya ingin menentukan distribusi . Ini agak analog dengan mengambil sampel 50 orang dari suatu populasi, mengamati ketinggian mereka, dan kemudian meminta distribusi ketinggian untuk seseorang "diambil secara acak dari populasi". Saya hanya menggunakan seragam sejak popuation dalam kasus ini disebutkan, dan saya melihat itu membingungkan. (Melanjutkan dengan analogi, dalam skenario ini, kita tidak bisa mengukur orang, cukup amati beberapa respons biner)

x

$x$

o_{x}

$o_x$

p_{x}

$p_x$

Fabio Beltramini

Saya telah menambahkan contoh di akhir pertanyaan yang mudah-mudahan memperjelas apa yang saya coba lakukan.

Fabio Beltramini

Terima kasih. Pertanyaannya tetap membingungkan. Dalam contoh Anda mengamati pasangan nilai yang diurutkan : waktu bergulir dan sukses. "Distribusi" apa yang ingin Anda tanyakan? Tepatnya properti apa dari kantong besar dadu itu sesuai?

whuber

Pasangan nilai (gulungan & keberhasilan) menyampaikan informasi tentang variabel laten, hal. Distribusi yang saya inginkan adalah distribusi p di dalam tas. Tentu saja, jika kita tidak mengasumsikan kelas distribusi yang dimilikinya, ada banyak kemungkinan sewenang-wenang, tetapi jika kita mengatakan bahwa p adalah misalnya beta-didistribusikan, maka turun ke hanya memilih parameter pemasangan terbaik untuk ini kelas distribusi. Jika itu membantu, sekarang pertimbangkan alih-alih dadu dengan p = 1 / 4,1 / 6,1 / 8, bahwa setiap objek memiliki ap dari Beta (2,2) ... atau Beta (4,4), dll.

Fabio Beltramini

Analisis Bayesian hierarkis tentang perbedaan proporsi

Jawaban: