Penyusutan yang sangat tinggi dalam regresi logistik efek acak

8

Perhatikan contoh sederhana berikut ini:

library( rms )
library( lme4 )
params <- structure(list(Ns = c(181L, 191L, 147L, 190L, 243L, 164L, 83L, 
                            383L, 134L, 238L, 528L, 288L, 214L, 502L, 307L, 302L, 199L, 156L, 
                            183L), means = c(0.09, 0.05, 0.03, 0.06, 0.07, 0.07, 0.1, 0.1, 
                                             0.06, 0.11, 0.1, 0.11, 0.07, 0.11, 0.1, 0.09, 0.1, 0.09, 0.08
                            )), .Names = c("Ns", "means"), row.names = c(NA, -19L), class = "data.frame")
SimData <- data.frame( ID = as.factor( rep( 1:nrow( params ), params$Ns ) ),
                   Res = do.call( c, apply( params, 1, function( x ) c( rep( 0, x[ 1 ]-round( x[ 1 ]*x[ 2 ] ) ),
                                                                        rep( 1, round( x[ 1 ]*x[ 2 ] ) ) ) ) ) )
tapply( SimData$Res, SimData$ID, mean )
dd <- datadist( SimData )
options( datadist = "dd" )
fitFE <- lrm( Res ~  ID, data = SimData )
fitRE <- glmer( Res ~ ( 1|ID ), data = SimData, family = binomial( link = logit ), nAGQ = 50 )

Yaitu kami memberikan efek tetap dan model efek acak untuk masalah yang sama, sangat sederhana (regresi logistik, intersep saja).

Ini adalah bagaimana model efek tetap terlihat seperti:

plot( summary( fitFE ) )

Model efek tetap

Dan ini adalah bagaimana efek acak:

dotplot( ranef( fitRE, condVar = TRUE ) )

Model efek acak

Penyusutan itu sendiri tidak mengejutkan, tetapi luasnya. Berikut ini adalah perbandingan yang lebih langsung:

xyplot( plogis(fe)~plogis(re),
    data = data.frame( re = coef( fitRE )$ID[ , 1 ],
                       fe = c( 0, coef( fitFE )[ -1  ] )+coef( fitFE )[ 1 ] ),
    abline = c( 0, 1 ) )

Probabilitas yang diprediksi dari model efek tetap dan acak

Perkiraan efek tetap berkisar antara kurang dari 3% hingga lebih dari 11, namun efek acaknya adalah antara 7,5 dan 9,5%. (Dimasukkannya kovariat membuat ini lebih ekstrim.)

Saya bukan ahli dalam efek acak dalam regresi logistik, tetapi dari regresi linier, saya mendapat kesan bahwa penyusutan substansial hanya dapat terjadi dari ukuran kelompok yang sangat-sangat kecil. Namun, di sini, bahkan kelompok terkecil memiliki hampir seratus pengamatan, dan ukuran sampel mencapai di atas 500.

Apa alasannya? Atau apakah saya mengabaikan sesuatu ...?

EDIT (28 Jul 2017). Sesuai saran @Ben Bolker, saya mencoba apa yang terjadi jika responsnya berkelanjutan (sehingga kami menghapus masalah tentang ukuran sampel efektif , yang khusus untuk data binomial).

Baru SimDataoleh karena itu

SimData <- data.frame( ID = as.factor( rep( 1:nrow( params ), params$Ns ) ),
                   Res = do.call( c, apply( params, 1, function( x ) c( rep( 0, x[ 1 ]-round( x[ 1 ]*x[ 2 ] ) ),
                                                                        rep( 1, round( x[ 1 ]*x[ 2 ] ) ) ) ) ),
                   Res2 = do.call( c, apply( params, 1, function( x ) rnorm( x[1], x[2], 0.1 ) ) ) )
data.frame( params, Res = tapply( SimData$Res, SimData$ID, mean ), Res2 = tapply( SimData$Res2, SimData$ID, mean ) )

dan model-model baru

fitFE2 <- ols( Res2 ~ ID, data = SimData )
fitRE2 <- lmer( Res2 ~ ( 1|ID ), data = SimData )

Hasilnya dengan

xyplot( fe~re, data = data.frame( re = coef( fitRE2 )$ID[ , 1 ],
                       fe = c( 0, coef( fitFE2 )[ -1  ] )+coef( fitFE2 )[ 1 ] ),
    abline = c( 0, 1 ) )

adalah

masukkan deskripsi gambar di sini

Sejauh ini bagus!

Namun, saya memutuskan untuk melakukan pemeriksaan lagi untuk memverifikasi ide Ben, tetapi hasilnya ternyata sangat aneh. Saya memutuskan untuk memeriksa teorinya dengan cara lain: Saya kembali ke hasil biner, tetapi meningkatkan sarana sehingga ukuran sampel yang efektif menjadi lebih besar. Saya hanya berlari params$means <- params$means + 0.5dan mencoba lagi contoh aslinya, inilah hasilnya:

masukkan deskripsi gambar di sini

Meskipun ukuran sampel minimum (efektif) memang meningkat drastis ...

> summary(with(SimData,tapply(Res,list(ID),
+                             function(x) min(sum(x==0),sum(x==1)))))
Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
33.0    72.5    86.0   100.3   117.5   211.0 

... penyusutan benar-benar meningkat ! (Menjadi total, dengan estimasi nol varians.)

Tamas Ferenci
sumber
Anda merencanakan rasio odds di plot pertama dan rasio odds log di plot kedua.
Douglas Bates
Ya, tetapi plot ketiga, yang benar-benar membandingkan mereka, dan menunjukkan masalah pertanyaan ini, menggunakan skala yang sama untuk keduanya! Sama seperti perintah verbal saya di bawah plot.
Tamas Ferenci

Jawaban:

5

Saya menduga bahwa jawabannya di sini berkaitan dengan definisi "ukuran sampel efektif". Aturan praktis (dari buku Harrell's Regression Modeling Strategies ) adalah bahwa ukuran sampel yang efektif untuk variabel Bernoulli adalah minimum dari jumlah keberhasilan dan kegagalan; mis. sampel ukuran 10.000 dengan hanya 4 keberhasilan lebih seperti memiliki daripada . Ukuran sampel yang efektif di sini tidak kecil, tetapi mereka jauh lebih kecil dari jumlah pengamatan.n=4n=104

Ukuran sampel efektif per kelompok:

summary(with(SimData,tapply(Res,list(ID),
                      function(x) min(sum(x==0),sum(x==1)))))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   4.00   11.00   16.00   21.63   29.00   55.00 

Ukuran sampel per grup:

summary(c(table(SimData$ID)))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   83.0   172.5   199.0   243.8   295.0   528.0 

Salah satu cara untuk menguji penjelasan ini adalah dengan melakukan contoh analog dengan respon yang terus menerus bervariasi (Gamma atau Gaussian).

Ben Bolker
sumber
2
wow, ukuran sampel yang efektif, saya tidak pernah memikirkan hal ini. Terima kasih! Eksperimen saya dengan respons Gaussian mengkonfirmasi gagasan Anda, tetapi meningkatkan ukuran sampel minimum efektif tidak; lihat hasil edit saya ...
Tamas Ferenci