Asimtot regresi binomial

Regresi logistik binomial memiliki asimtot atas dan bawah masing-masing 1 dan 0. Namun, data akurasi (hanya sebagai contoh) mungkin memiliki asimtot atas dan bawah yang sangat berbeda dengan 1 dan / atau 0. Saya dapat melihat tiga solusi potensial untuk ini:

1. Jangan khawatir tentang hal itu jika Anda mendapatkan pasangan yang cocok dalam bidang yang diminati. Jika Anda tidak mendapatkan kecocokan yang baik maka:
2. Transformasikan data sehingga jumlah minimum dan maksimum dari respons yang benar dalam sampel memberikan proporsi 0 dan 1 (bukannya katakan 0 dan 0,15).
   atau
3. Gunakan regresi non-linier sehingga Anda dapat menentukan asymptotes atau meminta penggantinya untuk Anda.

Tampaknya bagi saya bahwa opsi 1 & 2 akan lebih disukai daripada opsi 3 sebagian besar untuk alasan kesederhanaan, dalam hal mana opsi 3 mungkin merupakan opsi yang lebih baik karena dapat menghasilkan lebih banyak informasi?

sunting
Ini contohnya. Total yang mungkin benar untuk akurasi adalah 100, tetapi akurasi maksimum dalam hal ini adalah ~ 15.

accuracy <- c(0,0,0,0,0,1,3,5,9,13,14,15,14,15,16,15,14,14,15)
x<-1:length(accuracy)
glmx<-glm(cbind(accuracy, 100-accuracy) ~ x, family=binomial)
ndf<- data.frame(x=x)
ndf$fit<-predict(glmx, newdata=ndf, type="response")
plot(accuracy/100 ~ x)
with(ndf, lines(fit ~ x))

Opsi 2 (sesuai komentar dan untuk memperjelas makna saya) kemudian akan menjadi model

glmx2<-glm(cbind(accuracy, 16-accuracy) ~ x, family=binomial)

Opsi 3 (untuk kelengkapan) akan menjadi sesuatu yang mirip dengan:

fitnls<-nls(accuracy ~ upAsym + (y0 - upAsym)/(1 + (x/midPoint)^slope), 
  start = list("upAsym" = max(accuracy), "y0" = 0, "midPoint" = 10, "slope" = 5), 
  lower = list("upAsym" = 0, "y0" = 0, "midPoint" = 1, "slope" = 0), 
  upper = list("upAsym" = 100, "y0" = 0, "midPoint" = 19, hillslope = Inf), 
  control = nls.control(warnOnly = TRUE, maxiter=1000),
  algorithm = "port")

Pertanyaan menarik. Kemungkinan yang muncul di benak saya adalah memasukkan parameter tambahan untuk mengontrol batas atas fungsi 'tautan'.$p\in[0,1]$

Biarkan , menjadi pengamatan independen, di mana , , adalah vektor dari variabel penjelas, adalah vektor koefisien regresi dan adalah fungsi tautan. Kemudian fungsi kemungkinan diberikan oleh$\{{\bf x}_j,y_j,n_j\}$$j=1,...,n$$y_j\sim \text{Binomial}\{n_i,pF({\bf x}_j^T\beta)\}$$p\in[0,1]$${\bf x}_j=(1,x_{j1},...,x_{jk})^T$$\beta=(\beta_0,...,\beta_k)$$F^{-1}$

Langkah selanjutnya adalah memilih tautan, ucapkan distribusi logistik, dan temukan MLE yang sesuai .$(\beta,p)$

Pertimbangkan contoh mainan simulasi berikut menggunakan model dosis-respons dengan dan$(\beta_0,\beta_1,p)=(0.5,0.5,0.25)$$n=31$

dose = seq(-15,15,1)
a = 0.5
b = 0.5
n=length(dose)
sim = rep(0,n)
for(i in 1:n) sim[i] = rbinom(1,100,0.25*plogis(a+b*dose[i]))

plot(dose,sim/100)

lp = function(par){
if(par[3]>0& par[3]<1) return(-(n*mean(sim)*log(par[3]) +  sum(sim*log(plogis(par[1]+par[2]*dose)))  + sum((100-sim)*log(1-par[3]*plogis(par[1]+par[2]*dose))) ))
else return(-Inf)
}

optim(c(0.5,0.5,0.25),lp)

Salah satu hasil yang saya dapatkan adalah . Karena itu tampaknya akurat. Tentu saja, eksplorasi yang lebih rinci dari model ini akan diperlukan karena memasukkan parameter dalam model regresi biner bisa rumit dan masalah pengidentifikasian atau keberadaan MLE dapat melompat pada tahap 1 2 .$(\hat\beta_0,\hat\beta_1,\hat p)=( 0.4526650, 0.4589112, 0.2395564)$

Edit

Diberikan hasil edit (yang mengubah masalah secara signifikan), metode yang saya usulkan sebelumnya dapat dimodifikasi agar sesuai dengan data yang Anda berikan. Pertimbangkan modelnya

di mana adalah CDF logistik, adalah parameter lokasi, adalah parameter skala, dan parameter mengontrol ketinggian kurva sama seperti pada model sebelumnya. Model ini dapat dipasang menggunakan Nonlinear Least Squares . Kode R berikut menunjukkan cara melakukan ini untuk data Anda.$F$$\mu$$\sigma$$p$

rm(list=ls())
y = c(0,0,0,0,0,1,3,5,9,13,14,15,14,15,16,15,14,14,15)/100
x = 1:length(y)
N = length(y)

plot(y ~ x)

Data = data.frame(x,y)

nls_fit = nls(y ~ p*plogis(x,m,s), Data, start = list(m = 10, s = 1,  p = 0.2) )

lines(Data$x, predict(nls_fit), col = "red")

Saya akan menggunakan maksimum vektor X sebagai jumlah total kemungkinan keberhasilan. (Ini adalah perkiraan bias dari jumlah keberhasilan sebenarnya yang sebenarnya, tetapi ini akan bekerja dengan cukup baik jika Anda memiliki cukup data).

accuracy <- c(0,0,0,0,0,1,3,5,9,13,14,15,14,15,16,15,14,14,15)
x<-1:length(accuracy)
glmx<-glm(cbind(accuracy, max(accuracy)-accuracy) ~ x, family=binomial)
ndf<- data.frame(x=x)
ndf$fit<-predict(glmx, newdata=ndf, type="response")
plot(accuracy/max(accuracy) ~ x)
with(ndf, lines(fit ~ x))

Ini menciptakan plot yang terlihat seperti:

Perhatikan bahwa regresi binomial didasarkan pada memiliki respons biner untuk setiap kasus individu. setiap respons individu harus mampu mengambil satu dari dua nilai. Jika ada beberapa batasan proporsi maka pasti ada beberapa kasus yang hanya bisa mengambil satu nilai.

Sepertinya Anda tidak berurusan dengan data biner tetapi dengan data pada rentang yang terbatas. jika ini masalahnya, maka regresi beta terdengar lebih tepat. Kami dapat menulis distribusi beta sebagai:

Anda kemudian mengatur sama dengan fungsi tautan apa pun yang memetakan interval ke dalam real. Ada paket R yang dapat digunakan agar sesuai dengan model ini, meskipun saya pikir Anda perlu tahu batas-batasnya. Jika ya, maka definisikan ulang variabel baru .$g(\mu_i)=x_i^T\beta$$[L,U]$$y_i=\frac{d_i-L}{U-L}$