Apakah ada cara untuk memaksa hubungan antara koefisien dalam regresi logistik?

Saya ingin menentukan model regresi logistik di mana saya memiliki hubungan berikut:

$E[Y_i|X_i] = f(\beta x_{i1} + \beta^2x_{i2})$ dimana $f$ adalah fungsi invit logit.

Apakah ada cara "cepat" untuk melakukan ini dengan fungsi R yang sudah ada atau apakah ada nama untuk model seperti ini? Saya menyadari bahwa saya dapat memodifikasi algoritma Newton-Raphson yang digunakan untuk regresi logistik, tetapi ini adalah banyak teori dan pekerjaan pengkodean dan saya sedang mencari jalan pintas.

EDIT: dapatkan estimasi poin untuk $\beta$ cukup mudah menggunakan optim () atau pengoptimal lain di R untuk memaksimalkan kemungkinan. Tapi saya butuh kesalahan standar pada orang-orang ini.

r logistic generalized-linear-model nonlinear-regression TrynnaDoStat
sumber

Bisakah Anda menjelaskan situasi di mana Anda ingin melakukan itu selamanya? Saya hanya penasaran. Saya pribadi belum melihat ini dan mungkin harus kode itu menggunakan optimasi terbatas.

wolfsatthedoor

Saya tidak bisa masuk ke spesifik tetapi alasan saya ingin melakukan ini pada dasarnya adalah untuk kekuatan statistik. Jika saya percaya hubungan ini ada maka memaksa model untuk memiliki hubungan ini akan membuat saya lebih dekat dengan yang sebenarnya

β

$\beta$ nilai. Adapun untuk mendapatkan estimasi titik pada

β

$\beta$ , itu cukup mudah menggunakan optim () atau pengoptimal lain di R untuk memaksimalkan kemungkinan. Tapi saya butuh kesalahan standar pada orang-orang ini.

TrynnaDoStat

Yang ini tidak sesulit kelihatannya: cukup mudah untuk mendapatkan MLE dari prinsip pertama, karena yang Anda miliki hanyalah satu parameter. Tuliskan kemungkinan log dan bedakan sehubungan dengan

β

$\beta$ . Temukan nol. Itu akan dilakukan secara numerik. Jika Anda memiliki masalah memulai pencarian, pas dengan model dua parameter yang biasa dan gunakan (katakanlah) koefisien

x_{2}

$x_2$ sebagai nilai awal.

whuber

Prosedur NLIN dalam SAS dapat ditampung dengan rumus regresi nonlinier seperti ini dan akan menghasilkan koefisien kesalahan dan kesalahan. Apakah Anda menikah dengan R atau ada fleksibilitas?

RobertF

Mendapatkan SES tidak lebih sulit: ini adalah kasus mudah dari teori standar, di mana hanya ada satu parameter. Tetapi mengingat sifat nonlinear ketergantungannya pada parameter, saya akan enggan bergantung pada teori perkiraan atau pada optimasi numerik brute-force: plot kemungkinan fungsi itu sendiri, setidaknya dalam beberapa kasus, sehingga Anda dapat memahami kualitatifnya perilaku dekat puncak.

whuber

Jawaban:

Ini cukup mudah dilakukan dengan fungsi optimal di R. Pemahaman saya adalah bahwa Anda ingin menjalankan regresi logistik di mana y adalah biner. Anda cukup menulis fungsi dan kemudian memasukkannya ke dalam optim. Di bawah ini adalah beberapa kode yang tidak saya jalankan (kode semu).

#d is your data frame and y is normalized to 0,1
your.fun=function(b)
{

    EXP=exp(d$x1*b +d$x2*b^2)

    VALS=( EXP/(1+EXP) )^(d$y)*( 1/(1+EXP) )^(1-d$y) 
    return(-sum(log(VALS)))
}

result=optim(0,your.fun,method="BFGS",hessian=TRUE)
# estimates 
 result$par
    #standard errors
    sqrt(diag(inv(result$hessian)))
# maximum log likelihood
-result$value

Perhatikan bahwa your.fun adalah negatif dari fungsi kemungkinan log. Jadi optim memaksimalkan log kemungkinan (secara default optim meminimalkan semua yang mengapa saya membuat fungsi negatif). Jika Y bukan biner, buka http://fisher.osu.edu/~schroeder.9/AMIS900/ch5.pdf untuk formulir fungsi multinomial dan kondisional dalam model logit.

Zachary Blumenfeld
sumber

Hargai tanggapan Zachary! Namun, ini tidak akan berfungsi karena saya memerlukan kesalahan standar pada perkiraan saya. Saya berpikir untuk menggabungkan bootstrap dan optim () tetapi akan lebih suka metode non-bootstrap jika memungkinkan. Memodifikasi Newton-Raphson akan jauh lebih memuaskan tetapi jauh lebih sulit untuk diterapkan.

TrynnaDoStat

Mungkin saya tidak mengerti, kesalahan standar estimasi berasal dari fungsi kemungkinan maksimum goni dievaluasi pada estimasi. Cara Anda menulis fungsi Anda, Anda hanya memiliki satu parameter. Anda tentu dapat mem-bootstrap kode di atas untuk mendapatkan kesalahan standar juga.

Zachary Blumenfeld

@ ZacharyBlumenfeld Saya mengerti apa yang Anda katakan sekarang! Saya bingung karena pemahaman saya tentang teori asimptotik MLE adalah bahwa pengamatan kami harus benar (rata-rata tentu bervariasi di sini sehingga pengamatan saya tidak benar). Namun, saya melihat sekarang bahwa pengamatan tidak harus iid dalam kondisi keteraturan tertentu ( en.wikipedia.org/wiki/Maximum_likelihood#Asymptotic_normality ). Sekarang saya hanya perlu mengecek situasi saya memenuhi persyaratan keteraturan. Terima kasih lagi!

TrynnaDoStat

Catatan: Jika

Y = X β + ϵ

$Y = X\beta + \epsilon$ maka itu

ϵ

$\epsilon$ yang mungkin dianggap iid, belum tentu

Y

$Y$ .

conjugateprior

Jawaban di atas benar. Untuk referensi, berikut adalah beberapa kode R yang berfungsi untuk menghitungnya. Saya telah mengambil kebebasan untuk menambahkan intersep, karena Anda mungkin memang menginginkan salah satunya.

## make some data
set.seed(1234)
N <- 2000
x1 <- rnorm(N)
x2 <- rnorm(N)
## create linear predictor
lpred <- 0.5 + 0.5 * x1 + 0.25 * x2
## apply inverse link function
ey <- 1/(1 + exp(-lpred))
## sample some dependent variable
y <- rbinom(N, prob=ey, size=rep(1,N))

dat <- matrix(c(x1, x2, y), nrow=N, ncol=3)
colnames(dat) <- c('x1', 'x2', 'y')

Sekarang buat fungsi kemungkinan log untuk memaksimalkan, menggunakan di sini dbinomkarena ada, dan menjumlahkan hasilnya

## the log likelihood function
log.like <- function(beta, dat){
  lpred <- beta[1] + dat[,'x1'] * beta[2] + dat[,'x2'] * beta[2]**2
  ey <- 1/(1 + exp(-lpred))
  sum(dbinom(dat[,'y'], prob=ey, size=rep(1,nrow(dat)), log=TRUE))
}

dan sesuaikan model dengan kemungkinan maksimum. Saya belum repot-repot menawarkan gradien atau memilih metode optimasi, tetapi Anda mungkin ingin melakukan keduanya.

## fit
res <- optim(par=c(1,1), ## starting values 
             fn=log.like,
             control=list(fnscale=-1), ## maximise not minimise
             hessian=TRUE, ## for SEs
             dat=dat)

Sekarang lihat hasilnya. Estimasi parameter ML dan gejala asimtotik adalah:

## results
data.frame(coef=res$par,
           SE=sqrt(diag(solve(-res$hessian))))

yang seharusnya

##        coef         SE
## 1 0.4731680 0.04828779
## 2 0.5799311 0.03363505

atau ada bug (yang selalu memungkinkan).

Peringatan biasa tentang kesalahan standar turunan Hessian berlaku.

conjugateprior
sumber