Format input untuk respons dalam binomial glm di R

Di R, ada tiga metode untuk memformat data input untuk regresi logistik menggunakan glmfungsi:

Data dapat dalam format "biner" untuk setiap pengamatan (misalnya, y = 0 atau 1 untuk setiap pengamatan);
Data dapat dalam format "Wilkinson-Rogers" (misalnya, y = cbind(success, failure)) dengan setiap baris mewakili satu perlakuan; atau
Data dapat dalam format tertimbang untuk setiap pengamatan (misalnya, y = 0,3, bobot = 10).

Ketiga pendekatan menghasilkan estimasi koefisien yang sama, tetapi berbeda dalam derajat kebebasan dan menghasilkan nilai penyimpangan dan skor AIC. Dua metode terakhir memiliki lebih sedikit pengamatan (dan karenanya derajat kebebasan) karena mereka menggunakan setiap perlakuan untuk jumlah pengamatan sedangkan yang pertama menggunakan setiap pengamatan untuk jumlah pengamatan.

Pertanyaan saya: Apakah ada kelebihan numerik atau statistik untuk menggunakan satu format input daripada yang lain? Satu-satunya keuntungan yang saya lihat adalah tidak harus memformat ulang data seseorang Runtuk digunakan dengan model.

Saya telah melihat dokumentasi glm , mencari di web, dan situs ini dan menemukan satu posting yang berhubungan secara tangensial , tetapi tidak ada panduan tentang topik ini.

Berikut adalah contoh simulasi yang menunjukkan perilaku ini:

# Write function to help simulate data
drc4 <- function(x, b =1.0, c = 0, d = 1, e = 0){
    (d - c)/ (1 + exp(-b * (log(x)  - log(e))))
}
# simulate long form of dataset
nReps = 20
dfLong <- data.frame(dose = rep(seq(0, 10, by = 2), each = nReps))
dfLong$mortality <-rbinom(n = dim(dfLong)[1], size = 1,
                              prob = drc4(dfLong$dose, b = 2, e = 5))

# aggregate to create short form of dataset
dfShort <- aggregate(dfLong$mortality, by = list(dfLong$dose), 
                     FUN = sum)
colnames(dfShort) <- c("dose", "mortality")
dfShort$survival <- nReps - dfShort$mortality 
dfShort$nReps <- nReps
dfShort$mortalityP <- dfShort$mortality / dfShort$nReps

fitShort <- glm( cbind(mortality, survival) ~ dose, 
                 data = dfShort, 
                 family = "binomial")
summary(fitShort)

fitShortP <- glm( mortalityP ~ dose, data = dfShort, 
                  weights = nReps,     
                  family = "binomial")
summary(fitShortP)

fitLong <- glm( mortality ~ dose, data = dfLong, 
                family = "binomial")
summary(fitLong)

r logistic generalized-linear-model Richard Erickson
sumber

Dalam contoh Anda, perbedaan antara penyimpangan nol dan residual sama untuk ketiga model. Jika Anda menambah atau menghapus parameter, perubahan AIC juga sama untuk ketiganya.

Jonny Lomond

Anda menjawab sendiri: jika menggunakan data yang sama, untuk hasil yang sama, lalu bagaimana mereka bisa berbeda? Selain itu, jika metode ini akan memberikan hasil yang berbeda hanya karena menyediakan data dalam format berbeda, sesuatu akan sangat salah dengan itu, atau dengan implementasinya.

Tim

Format WR pada akhirnya kemungkinan tertimbang. Masalah dengan bobot adalah R tidak dapat menentukan apakah itu bobot frekuensi, bobot probabilitas, atau lainnya. Dengan bobot survei, misalnya, Anda mungkin hanya memiliki n pengamatan tetapi mereka mewakili segmen populasi / kerangka sampling. Jadi derajat kebebasannya memang 100. svyglmdari paket survei memberi Anda metode yang lebih baik untuk menangani argumen bobot.

AdamO

Tetapi jika Anda akan cocok dengan model quasibinomial menggunakan ketiga cara pengkodean mereka akan menghasilkan hasil yang berbeda, kan, karena seseorang bisa memiliki penyebaran berlebihan positif ketika dikodekan sebagai data binomial tetapi tidak ketika dikodekan sebagai data logistik / biner. Atau apakah saya salah dalam hal ini?

Tom Wenseleers

Jawaban:

Tidak ada alasan statistik untuk memilih satu dari yang lain, selain kejelasan konseptual. Meskipun nilai penyimpangan yang dilaporkan berbeda, perbedaan ini sepenuhnya karena model jenuh. Jadi perbandingan apa pun yang menggunakan penyimpangan relatif antara model tidak terpengaruh, karena model jenuh log-kemungkinan dibatalkan.

Saya pikir ini berguna untuk pergi melalui perhitungan penyimpangan eksplisit.

Penyimpangan model adalah 2 * (LL (Model Jenuh) - LL (Model)). Misalkan Anda memiliki sel yang berbeda, di mana adalah jumlah pengamatan dalam sel , adalah model prediksi untuk semua pengamatan di sel , dan adalah nilai yang diamati (0 atau 1) untuk pengamatan ke- di sel . $i$ $n_i$ $i$ $p_i$ $i$ $y_{ij}$ $j$ $i$

Bentuk panjang

Kemungkinan log dari model (diusulkan atau null) adalah

\sum_{i} \sum_{j} (\log (p_{i}) y_{i j} + \log (1 - p_{i}) (1 - y_{i j}))

$\sum_i\sum_j\left(\log(p_i)y_{ij} + \log(1 - p_i)(1 - y_{ij})\right)$

dan kemungkinan log dari model jenuh adalahIni sama dengan 0, karena adalah 0 atau 1. Note tidak terdefinisi, tetapi untuk kenyamanan silakan baca sebagai singkatan untuk , yaitu 0.

\sum_{i} \sum_{j} (\log (y_{i j}) y_{i j} + \log (1 - y_{i j}) (1 - y_{i j})) .

$\sum_i\sum_j \left(\log(y_{ij})y_{ij} + \log(1 - y_{ij})(1-y_{ij})\right).$

y_{i j}

$y_{ij}$

\log (0)

$\log(0)$

0 \log (0)

$0\log(0)$

lim_{x \to 0^{+}} x \log (x)

$\lim_{x \to 0^+}x\log(x)$

Bentuk pendek (berbobot)

Perhatikan bahwa distribusi binomial sebenarnya tidak dapat mengambil nilai-nilai non-integer, tetapi kami tetap dapat menghitung "log likelihood" dengan menggunakan sebagian kecil dari keberhasilan yang diamati dalam setiap sel sebagai respons, dan menimbang setiap musim panas dalam perhitungan log-likelihood dengan jumlah pengamatan dalam sel itu.

\sum_{i} n_{i} (\log (p_{i}) \sum_{j} y_{i j} / n_{i} + \log (1 - p_{i}) (1 - \sum_{j} (y_{i j} / n_{i}))

$\sum_in_i \left(\log(p_i)\sum_jy_{ij}/n_i + \log(1 - p_i)(1 - \sum_j(y_{ij}/n_i)\right)$

Ini persis sama dengan penyimpangan model yang kami hitung di atas, yang dapat Anda lihat dengan menarik jumlah di atas dalam persamaan bentuk panjang sejauh mungkin. $j$

Sementara penyimpangan jenuh berbeda. Karena kita tidak lagi memiliki 0-1 tanggapan, bahkan dengan satu parameter per pengamatan kita tidak bisa mendapatkan tepat 0. Sebaliknya model log-kemungkinan jenuh adalah

\sum_{i} n_{i} (\log (\sum_{j} y_{i j} / n_{i}) \sum_{j} y_{i j} / n_{i} + \log (1 - \sum_{j} y_{i j} / n_{i}) (1 - \sum_{j} y_{i j} / n_{i})) .

$\sum_i n_i\left(\log(\sum_jy_{ij}/n_i)\sum_jy_{ij}/n_i + \log(1 - \sum_jy_{ij}/n_i)(1-\sum_jy_{ij}/n_i)\right).$

Dalam contoh Anda, Anda dapat memverifikasi bahwa dua kali jumlah ini adalah perbedaan antara nilai null dan nilai penyimpangan yang dilaporkan untuk kedua model.

ni = dfShort$nReps
yavg = dfShort$mortalityP
sum.terms <-ni*(log(yavg)*yavg + log(1 - yavg)*(1 - yavg))
# Need to handle NaN when yavg is exactly 0
sum.terms[1] <- log(1 - yavg[1])*(1 - yavg[1])

2*sum(sum.terms)
fitShortP$deviance - fitLong$deviance

Jonny Lomond
sumber

Saya pikir Anda harus mengklarifikasi ekspresi untuk penyimpangan model jenuh. Log 0 tidak beroperasi dengan baik.

AdamO

Terima kasih, saya seharusnya menjelaskan apa yang saya maksud. Saya menambahkan edit untuk menjelaskan bahwa dengan 0log (0) maksud saya 0 dalam konteks ini.

Jonny Lomond

OK tapi saya benar-benar bingung (maafkan saya, saya tidak pernah membahas penyimpangan dengan sangat rinci): jika Anda memiliki log (y) y - log (1-y) (1-y) sebagai penyimpangan model jenuh, tidak setiap observasi hanya 0?

AdamO

"Model jenuh" adalah model yang dibayangkan dengan satu parameter per observasi. Jadi probabilitas yang diprediksi untuk setiap pengamatan adalah 1 atau 0, tergantung pada nilai sebenarnya yang diamati. Jadi dalam hal ini kemungkinan log dari model jenuh adalah 0, data adalah satu-satunya data yang dapat dihasilkan oleh model jenuh.

Jonny Lomond

Tom Wenseleers