R: fungsi glm dengan spesifikasi family = "binomial" dan "weight"

14

Saya sangat bingung dengan bagaimana berat bekerja di glm dengan family = "binomial". Dalam pemahaman saya, kemungkinan glm dengan family = "binomial" ditentukan sebagai berikut: mana adalah "proporsi keberhasilan yang diamati" dan adalah jumlah percobaan yang diketahui.yn

f(y)=(nny)halny(1-hal)n(1-y)=exp(n[ycatatanhal1-hal-(-catatan(1-hal))]+catatan(nny))
yn

Dalam pemahaman saya, probabilitas keberhasilan adalah parametrized dengan beberapa koefisien linier sebagai dan fungsi glm dengan pencarian family = "binomial" untuk: Maka masalah optimasi ini dapat disederhanakan sebagai:β p = p ( β ) arg max β ihalβhal=hal(β)

argmaksβsayacatatanf(ysaya).

argmaksβsayacatatanf(ysaya)=argmaksβsayansaya[ysayacatatanhal(β)1-hal(β)-(-catatan(1-hal(β)))]+catatan(nsayansayaysaya)=argmaksβsayansaya[ysayacatatanhal(β)1-hal(β)-(-catatan(1-hal(β)))]

Oleh karena itu jika kita membiarkan nsaya=nsayac untuk semua saya=1,...,N untuk beberapa konstanta c , maka harus juga benar bahwa:
argmaksβsayacatatanf(ysaya)=argmaksβsayansaya[ysayacatatanhal(β)1-hal(β)-(-catatan(1-hal(β)))]
Dari sini, saya berpikir bahwa Scaling dari jumlah percobaan nsayadengan konstanta TIDAK mempengaruhi perkiraan kemungkinan maksimum β mengingat proporsi keberhasilan ysaya .

File bantuan glm mengatakan:

 "For a binomial GLM prior weights are used to give the number of trials 
  when the response is the proportion of successes" 

Oleh karena itu saya berharap bahwa penskalaan berat tidak akan mempengaruhi perkiraan β mengingat proporsi keberhasilan sebagai respons. Namun dua kode berikut mengembalikan nilai koefisien yang berbeda:

 Y <- c(1,0,0,0) ## proportion of observed success
 w <- 1:length(Y) ## weight= the number of trials
 glm(Y~1,weights=w,family=binomial)

Ini menghasilkan:

 Call:  glm(formula = Y ~ 1, family = "binomial", weights = w)

 Coefficients:
 (Intercept)  
      -2.197     

sementara jika saya kalikan semua bobot dengan 1000, koefisien yang diperkirakan berbeda:

 glm(Y~1,weights=w*1000,family=binomial)

 Call:  glm(formula = Y ~ 1, family = binomial, weights = w * 1000)

 Coefficients:
 (Intercept)  
    -3.153e+15  

Saya melihat banyak contoh lain seperti ini bahkan dengan beberapa skala dalam bobot. Apa yang terjadi disini?

FairyOnIce
sumber
3
Untuk apa nilainya, weightsargumen berakhir di dua tempat di dalam glm.fitfungsi (di glm.R ), yang merupakan pekerjaan di R: 1) dalam residu penyimpangan, dengan cara fungsi C binomial_dev_resids(dalam keluarga.c ) dan 2) pada langkah IWLS dengan cara Cdqrls(dalam lm.c ). Saya tidak tahu cukup C untuk bisa lebih membantu dalam melacak logika
shadowtalker
3
Periksa balasan di sini .
Stat
@ssdecontrol Saya membaca glm.fit di tautan yang Anda berikan kepada saya, tetapi saya tidak dapat menemukan di mana fungsi C "binomial_dev_resids" dipanggil di glm.fit. Apakah Anda keberatan jika Anda menunjukkannya?
FairyOnIce
@ssdecontrol Oh, maaf saya rasa saya mengerti. Setiap "keluarga" adalah daftar dan salah satu unsurnya adalah "dev.resids". Ketika saya mengetik binomial di konsol R, saya melihat definisi objek binomial dan memiliki garis: dev.resids <- function (y, mu, wt) .Call (C_binomial_dev_resids, y, mu, wt)
FairyOnIce

Jawaban:

4

Contoh Anda hanya menyebabkan kesalahan pembulatan di R. Bobot besar tidak berkinerja baik di glm. Benar bahwa penskalaan wdengan jumlah yang lebih kecil, seperti 100, mengarah ke perkiraan yang sama dengan yang tidak dihitung w.

Jika Anda ingin perilaku yang lebih andal dengan argumen bobot, coba gunakan svyglmfungsi dari surveypaket.

Lihat disini:

    > svyglm(Y~1, design=svydesign(ids=~1, weights=~w, data=data.frame(w=w*1000, Y=Y)), family=binomial)
Independent Sampling design (with replacement)
svydesign(ids = ~1, weights = ~w, data = data.frame(w = w * 1000, 
    Y = Y))

Call:  svyglm(formula = Y ~ 1, design = svydesign(ids = ~1, weights = ~w2, 
    data = data.frame(w2 = w * 1000, Y = Y)), family = binomial)

Coefficients:
(Intercept)  
     -2.197  

Degrees of Freedom: 3 Total (i.e. Null);  3 Residual
Null Deviance:      2.601 
Residual Deviance: 2.601    AIC: 2.843
AdamO
sumber
1

glm.fitfamily$initializeglm.fitWXXW

$intializeKode yang relevan adalah:

if (NCOL(y) == 1) {
    if (is.factor(y)) 
        y <- y != levels(y)[1L]
    n <- rep.int(1, nobs)
    y[weights == 0] <- 0
    if (any(y < 0 | y > 1)) 
        stop("y values must be 0 <= y <= 1")
    mustart <- (weights * y + 0.5)/(weights + 1)
    m <- weights * y
    if (any(abs(m - round(m)) > 0.001)) 
        warning("non-integer #successes in a binomial glm!")
}

Ini adalah versi sederhana glm.fityang menunjukkan poin saya

> #####
> # setup
> y <- matrix(c(1,0,0,0), ncol = 1)
> weights <- 1:nrow(y) * 1000
> nobs <- length(y)
> family <- binomial()
> X <- matrix(rep(1, nobs), ncol = 1) # design matrix used later
> 
> # set mu start as with family$initialize
> if (NCOL(y) == 1) {
+   n <- rep.int(1, nobs)
+   y[weights == 0] <- 0
+   mustart <- (weights * y + 0.5)/(weights + 1)
+   m <- weights * y
+   if (any(abs(m - round(m)) > 0.001)) 
+     warning("non-integer #successes in a binomial glm!")
+ }
> 
> mustart # starting value
             [,1]
[1,] 0.9995004995
[2,] 0.0002498751
[3,] 0.0001666111
[4,] 0.0001249688
> (eta <- family$linkfun(mustart))
          [,1]
[1,]  7.601402
[2,] -8.294300
[3,] -8.699681
[4,] -8.987322
> 
> #####
> # Start loop to fit
> mu <- family$linkinv(eta)
> mu_eta <- family$mu.eta(eta)
> z <- drop(eta + (y - mu) / mu_eta)
> w <- drop(sqrt(weights * mu_eta^2 / family$variance(mu = mu)))
> 
> # code is simpler here as (X^T W X) is a scalar
> X_w <- X * w
> (.coef <- drop(crossprod(X_w)^-1 * ((w * z) %*% X_w)))
[1] -5.098297
> (eta <- .coef * X)
          [,1]
[1,] -5.098297
[2,] -5.098297
[3,] -5.098297
[4,] -5.098297
> 
> # repeat a few times from "start loop to fit"

Kita dapat mengulangi bagian terakhir dua kali lagi untuk melihat bahwa metode Newton-Raphson berbeda:

> #####
> # Start loop to fit
> mu <- family$linkinv(eta)
> mu_eta <- family$mu.eta(eta)
> z <- drop(eta + (y - mu) / mu_eta)
> w <- drop(sqrt(weights * mu_eta^2 / family$variance(mu = mu)))
> 
> # code is simpler here as (X^T W X) is a scalar
> X_w <- X * w
> (.coef <- drop(crossprod(X_w)^-1 * ((w * z) %*% X_w)))
[1] 10.47049
> (eta <- .coef * X)
         [,1]
[1,] 10.47049
[2,] 10.47049
[3,] 10.47049
[4,] 10.47049
> 
> 
> #####
> # Start loop to fit
> mu <- family$linkinv(eta)
> mu_eta <- family$mu.eta(eta)
> z <- drop(eta + (y - mu) / mu_eta)
> w <- drop(sqrt(weights * mu_eta^2 / family$variance(mu = mu)))
> 
> # code is simpler here as (X^T W X) is a scalar
> X_w <- X * w
> (.coef <- drop(crossprod(X_w)^-1 * ((w * z) %*% X_w)))
[1] -31723.76
> (eta <- .coef * X)
          [,1]
[1,] -31723.76
[2,] -31723.76
[3,] -31723.76
[4,] -31723.76

Ini tidak terjadi jika Anda memulai weights <- 1:nrow(y)atau mengatakan weights <- 1:nrow(y) * 100.

Perhatikan bahwa Anda dapat menghindari perbedaan dengan menetapkan mustartargumen. Misalnya lakukan

> glm(Y ~ 1,weights = w * 1000, family = binomial, mustart = rep(0.5, 4))

Call:  glm(formula = Y ~ 1, family = binomial, weights = w * 1000, mustart = rep(0.5, 
    4))

Coefficients:
(Intercept)  
     -2.197  

Degrees of Freedom: 3 Total (i.e. Null);  3 Residual
Null Deviance:      6502 
Residual Deviance: 6502     AIC: 6504
Benjamin Christoffersen
sumber
Saya pikir bobot mempengaruhi lebih dari argumen untuk diinisialisasi. Dengan regresi logistik, Newton Raphson memperkirakan kemungkinan maksimum yang ada dan unik ketika data tidak dipisahkan. Memasok nilai awal yang berbeda ke pengoptimal tidak akan sampai pada nilai yang berbeda, tetapi mungkin akan memakan waktu lebih lama untuk sampai ke sana.
AdamO
"Memasok nilai awal yang berbeda ke pengoptimal tidak akan sampai pada nilai yang berbeda ..." . Yah metode Newton tidak menyimpang dan menemukan maksimum unik dalam contoh terakhir di mana saya menetapkan nilai awal (lihat contoh di mana saya memberikan mustart argumen). Sepertinya masalah yang terkait dengan perkiraan awal yang buruk .
Benjamin Christoffersen