Algoritma optimisasi mana yang digunakan dalam fungsi glm di R?

17

Seseorang dapat melakukan regresi logit di R menggunakan kode seperti:

> library(MASS)
> data(menarche)
> glm.out = glm(cbind(Menarche, Total-Menarche) ~ Age,
+                                              family=binomial(logit), data=menarche)
> coefficients(glm.out)
(Intercept)         Age 
 -21.226395    1.631968

Sepertinya algoritme pengoptimalan telah terkonvergensi - ada informasi tentang jumlah langkah algoritme penilaian skor:

Call:
glm(formula = cbind(Menarche, Total - Menarche) ~ Age, family = binomial(logit), 
    data = menarche)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.0363  -0.9953  -0.4900   0.7780   1.3675  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -21.22639    0.77068  -27.54   <2e-16 ***
Age           1.63197    0.05895   27.68   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 3693.884  on 24  degrees of freedom
Residual deviance:   26.703  on 23  degrees of freedom
AIC: 114.76

Number of Fisher Scoring iterations: 4

Saya ingin tahu tentang algoritma optimal apa itu? Apakah ini algoritma Newton-Raphson (keturunan urutan kedua)? Bisakah saya mengatur beberapa parameter untuk menggunakan algoritma Cauchy (keturunan gradien urutan pertama)?

Marcin Kosiński
sumber
5
Apakah Anda keberatan mengutip di mana algoritma Newton-Raphson disebut sebagai "gradient descent level II"?
Cliff AB
4
Skor FIsher itu sendiri terkait dengan, tetapi berbeda dari Newton-Raphson, pada dasarnya menggantikan Hessian dengan nilai yang diharapkan di bawah model.
Glen_b -Reinstate Monica
@CliffAB sory. Saya ment itu Newton's methodadalah metode keturunan urutan gradien kedua.
Marcin Kosiński
1
@ hxd1011, Anda tidak boleh mengedit untuk mengubah pertanyaan orang lain. Adalah satu hal untuk diedit ketika Anda berpikir Anda tahu apa artinya, tetapi pertanyaan mereka tidak jelas (mungkin bahasa Inggris bukan bahasa ibu mereka, misalnya), tetapi Anda tidak boleh membuat pertanyaan mereka berbeda (misalnya, lebih umum) daripada yang mereka miliki ingin. Sebaliknya, ajukan pertanyaan baru dengan apa yang Anda inginkan. Saya memutar kembali edit.
gung - Reinstate Monica
1
@ MarcinKosiński Metode Newton adalah Newton-Raphson, Raphson hanya membangun ide-ide Newton untuk kasus yang lebih umum.
AdamO

Jawaban:

19

Anda akan tertarik untuk mengetahui bahwa dokumentasi untuk glm, diakses melalui ?glmmenyediakan banyak wawasan yang bermanfaat: di bawah methodkami menemukan bahwa kuadrat terkecil yang diulang-ulang secara berulang adalah metode default untuk glm.fit, yang merupakan fungsi kuda kerja untuk glm. Selain itu, dokumentasi menyebutkan bahwa fungsi yang ditentukan pengguna dapat disediakan di sini, bukan sebagai standar.

Sycorax berkata Reinstate Monica
sumber
3
Anda juga bisa mengetikkan nama fungsi glmatau fit.glmsaat Rdiminta untuk mempelajari kode sumber.
Matthew Drury
@MatthewDrury Saya pikir Anda berarti pekerja keras glm.fityang tidak akan sepenuhnya direproduksi karena hal itu bergantung pada kode C C_Cdqrls.
AdamO
Yah, Anda benar, saya banyak mencampur pesanan dalam R. Apa maksud Anda tidak dapat direproduksi?
Matthew Drury
16

Metode yang digunakan disebutkan dalam output itu sendiri: itu adalah Fisher Scoring. Ini setara dengan Newton-Raphson dalam banyak kasus. Pengecualiannya adalah situasi di mana Anda menggunakan parameterisasi non-alami. Regresi risiko relatif adalah contoh dari skenario semacam itu. Di sana, informasi yang diharapkan dan diamati berbeda. Secara umum, Newton Raphson dan Fisher Scoring memberikan hasil yang hampir identik.

pp(1p)Skor Fisher. Selain itu, ini memberikan beberapa intuisi yang bagus untuk algoritma EM yang merupakan kerangka kerja yang lebih umum untuk memperkirakan kemungkinan yang rumit.

Pengoptimal umum default di R menggunakan metode numerik untuk memperkirakan momen kedua, pada dasarnya didasarkan pada linierisasi (berhati-hatilah terhadap kutukan dimensi). Jadi jika Anda tertarik untuk membandingkan efisiensi dan bias, Anda dapat menerapkan rutin kemungkinan maksimum logistik naif dengan sesuatu seperti

set.seed(1234)
x <- rnorm(1000)
y <- rbinom(1000, 1, exp(-2.3 + 0.1*x)/(1+exp(-2.3 + 0.1*x)))
f <- function(b) {
  p <- exp(b[1] + b[2]*x)/(1+exp(b[1] + b[2]*x))
  -sum(dbinom(y, 1, p, log=TRUE))
}
ans <- nlm(f, p=0:1, hessian=TRUE)

memberi saya

> ans$estimate
[1] -2.2261225  0.1651472
> coef(glm(y~x, family=binomial))
(Intercept)           x 
 -2.2261215   0.1651474 
AdamO
sumber
apa perbedaannya dibandingkan dengan gradient decent / newton's method / BFGS? Saya pikir Anda menjelaskan, tetapi saya tidak cukup mengikuti.
Haitao Du
@ hxd1011 lihat "Metode Numerik untuk Optimasi Tanpa Batas dan Persamaan Nonlinier" Dennis, JE dan Schnabel, RB
AdamO
1
@ hxd1011 sejauh yang saya tahu, Newton Raphson tidak memerlukan atau memperkirakan Goni dalam langkah-langkahnya. Metode Newton-Type dalam nlmmemperkirakan gradien secara numerik kemudian menerapkan Newton Raphson. Dalam BFGS saya pikir gradien diperlukan seperti halnya dengan Newton Raphson, tetapi langkah-langkah berturut-turut dievaluasi menggunakan pendekatan orde kedua yang membutuhkan perkiraan Hessian. BFGS baik untuk optimasi yang sangat nonlinier. Tetapi untuk GLM, mereka biasanya berperilaku sangat baik.
AdamO
2
Jawaban yang ada menyebutkan "kuadrat terkecil yang berulang secara berulang", bagaimana cara memasukkan gambar, dibandingkan dengan algoritma yang Anda sebutkan di sini?
Amuba kata Reinstate Monica
@AdamO dapatkah Anda menjawab komentar amuba? Terima kasih
Haitao Du
1

Sebagai catatan, implementasi R sederhana murni dari algoritma R's glm, berdasarkan skoring Fisher (berulang kuadrat terkecil berulang), seperti yang dijelaskan dalam jawaban lain diberikan oleh:

glm_irls = function(X, y, weights=rep(1,nrow(X)), family=poisson(log), maxit=25, tol=1e-16) {
    if (!is(family, "family")) family = family()
    variance = family$variance
    linkinv = family$linkinv
    mu.eta = family$mu.eta
    etastart = NULL

    nobs = nrow(X)    # needed by the initialize expression below
    nvars = ncol(X)   # needed by the initialize expression below
    eval(family$initialize) # initializes n and fitted values mustart
    eta = family$linkfun(mustart) # we then initialize eta with this
    dev.resids = family$dev.resids
    dev = sum(dev.resids(y, linkinv(eta), weights))
    devold = 0
    beta_old = rep(1, nvars)

    for(j in 1:maxit)
    {
      mu = linkinv(eta) 
      varg = variance(mu)
      gprime = mu.eta(eta)
      z = eta + (y - mu) / gprime # potentially -offset if you would have an offset argument as well
      W = weights * as.vector(gprime^2 / varg)
      beta = solve(crossprod(X,W*X), crossprod(X,W*z), tol=2*.Machine$double.eps)
      eta = X %*% beta # potentially +offset if you would have an offset argument as well
      dev = sum(dev.resids(y, mu, weights))
      if (abs(dev - devold) / (0.1 + abs(dev)) < tol) break
      devold = dev
      beta_old = beta
    }
    list(coefficients=t(beta), iterations=j)
}

Contoh:

## Dobson (1990) Page 93: Randomized Controlled Trial :
y <- counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
X <- model.matrix(counts ~ outcome + treatment)

coef(glm.fit(x=X, y=y, family = poisson(log))) 
  (Intercept)      outcome2      outcome3    treatment2    treatment3 
 3.044522e+00 -4.542553e-01 -2.929871e-01 -7.635479e-16 -9.532452e-16

coef(glm_irls(X=X, y=y, family=poisson(log)))
     (Intercept)   outcome2   outcome3    treatment2   treatment3
[1,]    3.044522 -0.4542553 -0.2929871 -3.151689e-16 -8.24099e-16

Diskusi yang baik tentang algoritma pemasangan GLM, termasuk perbandingan dengan Newton-Raphson (yang menggunakan Hessian yang diamati sebagai kebalikan dari Hessian yang diharapkan dalam algoritma IRLS) dan algoritma hybrid (yang dimulai dengan IRLS, karena ini lebih mudah untuk diinisialisasi, tetapi kemudian selesaikan dengan optimasi lebih lanjut menggunakan Newton-Raphson) dapat ditemukan dalam buku "Generalized Linear Models and Extensions" oleh James W. Hardin & Joseph M. Hilbe .

Tom Wenseleers
sumber