Memeriksa asumsi peluang proporsional berlaku dalam regresi logistik ordinal menggunakan fungsi polr

Saya telah menggunakan fungsi 'polr' dalam paket MASS untuk menjalankan regresi logistik ordinal untuk variabel respon kategoris ordinal dengan 15 variabel penjelas kontinu.

Saya telah menggunakan kode (ditunjukkan di bawah) untuk memeriksa bahwa model saya memenuhi asumsi peluang proporsional mengikuti saran yang diberikan dalam panduan UCLA . Namun, saya sedikit khawatir tentang output yang menyiratkan bahwa tidak hanya koefisien di berbagai titik potong serupa, tetapi mereka persis sama (lihat grafik di bawah).

FGV1b <- data.frame(FG1_val_cat=factor(FGV1b[,"FG1_val_cat"]), 
                    scale(FGV1[,c("X","Y","Slope","Ele","Aspect","Prox_to_for_FG", 
                          "Prox_to_for_mL", "Prox_to_nat_border", "Prox_to_village", 
                          "Prox_to_roads", "Prox_to_rivers", "Prox_to_waterFG", 
                          "Prox_to_watermL", "Prox_to_core", "Prox_to_NR", "PCA1", 
                          "PCA2", "PCA3")]))
b     <- polr(FG1_val_cat ~ X + Y + Slope + Ele + Aspect + Prox_to_for_FG + 
                            Prox_to_for_mL + Prox_to_nat_border + Prox_to_village + 
                            Prox_to_roads + Prox_to_rivers + Prox_to_waterFG + 
                            Prox_to_watermL + Prox_to_core + Prox_to_NR, 
              data=FGV1b, Hess=TRUE)

Lihat ringkasan model:

summary(b)
(ctableb <- coef(summary(b)))
q        <- pnorm(abs(ctableb[, "t value"]), lower.tail=FALSE) * 2
(ctableb <- cbind(ctableb, "p value"=q))

Dan sekarang kita bisa melihat interval kepercayaan untuk estimasi parameter:

(cib <- confint(b)) 
confint.default(b)

Tetapi hasil ini masih cukup sulit untuk ditafsirkan, jadi mari kita konversi koefisien menjadi rasio odds

exp(cbind(OR=coef(b), cib))

Memeriksa asumsi tersebut. Jadi kode berikut akan memperkirakan nilai yang akan digambarkan. Pertama itu menunjukkan kepada kita transformasi logit dari probabilitas menjadi lebih besar atau sama dengan setiap nilai variabel target

FG1_val_cat <- as.numeric(FG1_val_cat)
sf <- function(y) {
  c('VC>=1' = qlogis(mean(FG1_val_cat >= 1)),
    'VC>=2' = qlogis(mean(FG1_val_cat >= 2)),
    'VC>=3' = qlogis(mean(FG1_val_cat >= 3)),
    'VC>=4' = qlogis(mean(FG1_val_cat >= 4)),
    'VC>=5' = qlogis(mean(FG1_val_cat >= 5)),
    'VC>=6' = qlogis(mean(FG1_val_cat >= 6)),
    'VC>=7' = qlogis(mean(FG1_val_cat >= 7)),
    'VC>=8' = qlogis(mean(FG1_val_cat >= 8)))
}
(t <- with(FGV1b, summary(as.numeric(FG1_val_cat) ~ X + Y + Slope + Ele + Aspect + 
                             Prox_to_for_FG + Prox_to_for_mL + Prox_to_nat_border + 
                             Prox_to_village + Prox_to_roads + Prox_to_rivers + 
                             Prox_to_waterFG + Prox_to_watermL + Prox_to_core + 
                             Prox_to_NR, fun=sf)))

Tabel di atas menampilkan nilai prediksi (linier) yang akan kita dapatkan jika kita regresi variabel dependen kita pada variabel prediktor kita satu per satu, tanpa asumsi kemiringan paralel. Jadi sekarang, kita dapat menjalankan serangkaian regresi logistik biner dengan beragam cutpoint pada variabel dependen untuk memeriksa kesetaraan koefisien di seluruh cutpoint

par(mfrow=c(1,1))
plot(t, which=1:8, pch=1:8, xlab='logit', main=' ', xlim=range(s[,7:8]))

Mohon maaf bahwa saya bukan ahli statistik dan mungkin saya kehilangan sesuatu yang jelas di sini. Namun, saya telah menghabiskan waktu lama untuk mencari tahu apakah ada masalah dalam cara saya menguji asumsi model dan juga mencoba mencari cara lain untuk menjalankan model yang sama.

Sebagai contoh, saya membaca di banyak milis bantuan yang lain menggunakan fungsi vglm (dalam paket VGAM) dan fungsi lrm (dalam paket rms) (misalnya lihat di sini: Asumsi peluang proporsional dalam regresi logistik ordinal dalam R dengan paket-paket VGAM dan rms ). Saya telah mencoba menjalankan model yang sama tetapi terus-menerus menghadapi peringatan dan kesalahan.

Misalnya, ketika saya mencoba menyesuaikan model vglm dengan argumen 'parallel = FALSE' (karena tautan sebelumnya menyebutkan penting untuk menguji asumsi odds proporsional), saya menemukan kesalahan berikut:

Kesalahan dalam lm.fit (X.vlm, y = z.vlm, ...): NA / NaN / Inf di 'y'
Selain itu: Pesan peringatan:
In Deviance.categorical.data.vgam (mu = mu, y = y, w = w, residual = residual,: nilai yang dipasang mendekati 0 atau 1

Saya ingin bertanya tolong jika ada orang yang mungkin mengerti dan dapat menjelaskan kepada saya mengapa grafik yang saya buat di atas terlihat seperti itu. Jika memang itu berarti ada sesuatu yang tidak beres, bisa tolong bantu saya menemukan cara untuk menguji asumsi peluang proporsional ketika hanya menggunakan fungsi polr. Atau jika itu tidak mungkin, maka saya akan mencoba menggunakan fungsi vglm, tetapi kemudian akan membutuhkan bantuan untuk menjelaskan mengapa saya terus mendapatkan kesalahan yang diberikan di atas.

CATATAN: Sebagai latar belakang, ada 1000 titik data di sini, yang sebenarnya merupakan titik lokasi di seluruh area studi. Saya mencari untuk melihat apakah ada hubungan antara variabel respon kategoris dan 15 variabel penjelas ini. Semua 15 variabel penjelas tersebut adalah karakteristik spasial (misalnya, ketinggian, koordinat xy, kedekatan dengan hutan, dll.). 1000 titik data dialokasikan secara acak menggunakan GIS, tetapi saya mengambil pendekatan pengambilan sampel bertingkat. Saya memastikan bahwa 125 poin dipilih secara acak dalam masing-masing dari 8 tingkat respons kategoris yang berbeda. Saya harap informasi ini juga bermanfaat.

Variabel dependen memiliki 8 level yang dipesan sehingga dalam grafik untuk menguji asumsi peluang proporsional, Anda akan melihat 8 simbol yang berbeda untuk setiap variabel independen. Anda hanya melihat 2 simbol untuk setiap variabel independen mungkin karena Anda telah memilih interval yang terlalu pendek untuk nilai sumbu x. Jika dugaan saya benar, Anda hanya perlu menggunakan interval yang lebih luas untuk nilai sumbu x. Coba kode ini:

par(mfrow=c(1,1))
plot(t, which=1:8, pch=1:8, xlab='logit', main=' ', xlim=range(s[,3:9]))

Jadi, saya menemukan ini melalui googling sekitar dan saya pikir jawaban mungkin masih berguna karena alasan itu. Saya pikir ada kesalahan

sf <- function(y) {
  c('VC>=1' = qlogis(mean(FG1_val_cat >= 1)),
    'VC>=2' = qlogis(mean(FG1_val_cat >= 2)),
    'VC>=3' = qlogis(mean(FG1_val_cat >= 3)),
    'VC>=4' = qlogis(mean(FG1_val_cat >= 4)),
    'VC>=5' = qlogis(mean(FG1_val_cat >= 5)),
    'VC>=6' = qlogis(mean(FG1_val_cat >= 6)),
    'VC>=7' = qlogis(mean(FG1_val_cat >= 7)),
    'VC>=8' = qlogis(mean(FG1_val_cat >= 8)))
}

di mana Anda menggunakan FG1_val_catdaripada y. Menggunakan contoh dari Strategi Pemodelan Regresi Harrell:

library(Hmisc)
getHdata(support)
support <- support[complete.cases(support[, c("sfdm2", "adlsc", "sex", "age", "meanbp")]), ]
sfdm <- as.integer (support$sfdm2 ) - 1

sf1 <- function (y) {
  c(' Y ≥ 1 ' = qlogis (mean(sfdm >= 1)), 
    ' Y ≥ 2 ' = qlogis (mean(sfdm >= 2)),
    ' Y ≥ 3 ' = qlogis (mean(sfdm >= 3))
  )
}

sf2 <- function (y) {
  c(' Y ≥ 1 ' = qlogis (mean(y >= 1)), 
    ' Y ≥ 2 ' = qlogis (mean(y >= 2)),
    ' Y ≥ 3 ' = qlogis (mean(y >= 3))
  )
}

s1 <- summary(sfdm ~ adlsc + sex + age + meanbp, fun=sf1,
              data = support)
s2 <- summary(sfdm ~ adlsc + sex + age + meanbp, fun=sf2,
              data = support)  

plot(s1, which =1:3, pch =1:3, xlab = ' logit ', main = ' ', width.factor = 1.4, cex.lab = 0.75)

plot(s2, which =1:3, pch =1:3, xlab = ' logit ', main = ' ', width.factor = 1.4, cex.lab = 0.75)

vs.