Nilai awal standar pas dengan regresi logistik dengan glm

Saya bertanya-tanya bagaimana nilai awal default ditentukan dalam glm.

Posting ini menunjukkan bahwa nilai default ditetapkan sebagai nol. Ini salah satu mengatakan bahwa ada sebuah algoritma di balik itu, namun link yang relevan rusak.

Saya mencoba menyesuaikan model regresi logistik sederhana dengan penelusuran algoritme:

set.seed(123)

x <- rnorm(100)
p <- 1/(1 + exp(-x))
y <- rbinom(100, size = 1, prob = p)

# to see parameter estimates in each step
trace(glm.fit, quote(print(coefold)), at = list(c(22, 4, 8, 4, 19, 3)))

Pertama, tanpa spesifikasi nilai awal:

glm(y ~ x, family = "binomial")

Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  .... step 22,4,8,4,19,3 
NULL
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  .... step 22,4,8,4,19,3 
[1] 0.386379 1.106234
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  .... step 22,4,8,4,19,3 
[1] 0.3991135 1.1653971
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  .... step 22,4,8,4,19,3 
[1] 0.3995188 1.1669508

Pada langkah pertama, nilai awal adalah NULL.

Kedua, saya menetapkan nilai awal menjadi nol:

glm(y ~ x, family = "binomial", start = c(0, 0))

Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  .... step 22,4,8,4,19,3 
[1] 0 0
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  .... step 22,4,8,4,19,3 
[1] 0.3177530 0.9097521
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  .... step 22,4,8,4,19,3 
[1] 0.3909975 1.1397163
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  .... step 22,4,8,4,19,3 
[1] 0.3994147 1.1666173
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  .... step 22,4,8,4,19,3 
[1] 0.3995191 1.1669518

Dan kita dapat melihat bahwa iterasi antara pendekatan pertama dan kedua berbeda.

Untuk melihat nilai awal yang ditentukan oleh glmsaya mencoba mencocokkan model dengan hanya satu iterasi:

glm(y ~ x, family = "binomial", control = list(maxit = 1))

Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  .... step 22,4,8,4,19,3 
NULL

Call:  glm(formula = y ~ x, family = "binomial", control = list(maxit = 1))

Coefficients:
(Intercept)            x  
     0.3864       1.1062  

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:      134.6 
Residual Deviance: 115  AIC: 119

Estimasi parameter (tidak mengherankan) sesuai dengan estimasi pendekatan pertama dalam iterasi kedua yaitu, [1] 0.386379 1.106234 Menetapkan nilai-nilai ini sebagai nilai awal mengarah ke urutan iterasi yang sama seperti pada pendekatan pertama:

glm(y ~ x, family = "binomial", start = c(0.386379, 1.106234))

Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  .... step 22,4,8,4,19,3 
[1] 0.386379 1.106234
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  .... step 22,4,8,4,19,3 
[1] 0.3991135 1.1653971
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  .... step 22,4,8,4,19,3 
[1] 0.3995188 1.1669508

Jadi pertanyaannya adalah, bagaimana nilai-nilai ini dihitung?

TL; DR

• start=c(b0,b1)menginisialisasi eta menjadi b0+x*b1(mu ke 1 / (1 + exp (-eta)))
• start=c(0,0) menginisialisasi eta ke 0 (mu ke 0,5) terlepas dari nilai y atau x.
• start=NULL menginisialisasi eta = 1.098612 (mu = 0.75) jika y = 1, terlepas dari nilai x.

• start=NULL menginisialisasi eta = -1.098612 (mu = 0.25) jika y = 0, terlepas dari nilai x.

• Setelah eta (dan akibatnya mu dan var (mu)) telah dihitung, wdan zdihitung dan dikirim ke pemecah QR, dalam semangat qr.solve(cbind(1,x) * w, z*w).

Bentuk panjang

Membangun off komentar Roland: Saya membuat glm.fit.truncated(), di mana aku mengambil glm.fitturun ke C_Cdqrlspanggilan, dan kemudian berkomentar itu. glm.fit.truncatedmenampilkan nilai zdan w(serta nilai jumlah yang digunakan untuk menghitung zdan w) yang kemudian akan diteruskan ke C_Cdqrlspanggilan:

## call Fortran code via C wrapper
fit <- .Call(C_Cdqrls, x[good, , drop = FALSE] * w, z * w,
             min(1e-7, control$epsilon/1000), check=FALSE) 

Lebih banyak dapat dibaca di C_Cdqrls sini . Untungnya, fungsi qr.solvepada basis R mengetuk langsung ke versi LINPACK yang dipanggil glm.fit().

Jadi kami menjalankan glm.fit.truncateduntuk spesifikasi nilai awal yang berbeda, dan kemudian melakukan panggilan ke qr.solvedengan nilai w dan z, dan kami melihat bagaimana "nilai awal" (atau nilai iterasi yang ditampilkan pertama) dihitung. Seperti yang ditunjukkan Roland, menentukan start=NULLatau start=c(0,0)dalam glm () memengaruhi perhitungan untuk w dan z, bukan untuk start.

Untuk awal = NULL: zadalah vektor di mana elemen memiliki nilai 2.431946 atau -2.431946 dan wmerupakan vektor di mana semua elemen adalah 0.4330127:

start.is.null <- glm.fit.truncated(x,y,family=binomial(), start=NULL)
start.is.null
w <- start.is.null$w
z <- start.is.null$z
## if start is NULL, the first displayed values are:
qr.solve(cbind(1,x) * w, z*w)  
# > qr.solve(cbind(1,x) * w, z*w)  
#                 x 
# 0.386379 1.106234 

Untuk awal = c (0,0): zadalah vektor di mana elemen memiliki nilai 2 atau -2 dan wmerupakan vektor di mana semua elemen 0,5:

## if start is c(0,0)    
start.is.00 <- glm.fit.truncated(x,y,family=binomial(), start=0)
start.is.00
w <- start.is.00$w
z <- start.is.00$z
## if start is c(0,0), the first displayed values are:    
qr.solve(cbind(1,x) * w, z*w)  
# > qr.solve(cbind(1,x) * w, z*w)  
#                   x 
# 0.3177530 0.9097521 

Jadi itu semua baik dan bagus, tetapi bagaimana kita menghitung wdan z? Di dekat bagian bawah glm.fit.truncated()kita lihat

z <- (eta - offset)[good] + (y - mu)[good]/mu.eta.val[good]
w <- sqrt((weights[good] * mu.eta.val[good]^2)/variance(mu)[good])

Lihatlah perbandingan berikut antara nilai yang dihasilkan dari jumlah yang digunakan untuk menghitung zdan w:

cbind(y, start.is.null$mu, start.is.00$mu)
cbind(y, start.is.null$eta, start.is.00$eta)
cbind(start.is.null$var_mu, start.is.00$var_mu)
cbind(start.is.null$mu.eta.val, start.is.00$mu.eta.val)

Perhatikan bahwa start.is.00akan memiliki vektor mudengan hanya nilai 0,5 karena eta diatur ke 0 dan mu (eta) = 1 / (1 + exp (-0)) = 0,5. start.is.nullmenetapkan yang dengan y = 1 menjadi mu = 0,75 (yang sesuai dengan eta = 1,098612) dan yang dengan y = 0 menjadi mu = 0,25 (yang sesuai dengan eta = -1,098612), dan dengan demikian var_mu= 0,75 * 0,25 = 0,1875.

Namun, menarik untuk dicatat, bahwa saya mengubah benih dan memutar ulang semuanya dan mu = 0,75 untuk y = 1 dan mu = 0,25 untuk y = 0 (dan dengan demikian jumlah lainnya tetap sama). Dengan kata lain, mulai = NULL memunculkan yang sama wdan zterlepas dari apa ydan apa x, karena mereka menginisialisasi eta = 1.098612 (mu = 0.75) jika y = 1 dan eta = -1.098612 (mu = 0.25) jika y = 0.

Jadi nampak bahwa nilai awal untuk koefisien Intercept dan untuk koefisien-X tidak diatur untuk mulai = NULL, melainkan nilai awal diberikan kepada eta tergantung pada nilai-y dan tidak tergantung pada nilai-x. Dari sana wdan zdihitung, lalu dikirim bersama xke qr.solver.

Kode untuk dijalankan sebelum potongan di atas:

set.seed(123)

x <- rnorm(100)
p <- 1/(1 + exp(-x))
y <- rbinom(100, size = 1, prob = p)


glm.fit.truncated <- function(x, y, weights = rep.int(1, nobs), 
start = 0,etastart = NULL, mustart = NULL, 
offset = rep.int(0, nobs),
family = binomial(), 
control = list(), 
intercept = TRUE,
singular.ok = TRUE
){
control <- do.call("glm.control", control)
x <- as.matrix(x)
xnames <- dimnames(x)[[2L]]
ynames <- if(is.matrix(y)) rownames(y) else names(y)
conv <- FALSE
nobs <- NROW(y)
nvars <- ncol(x)
EMPTY <- nvars == 0
## define weights and offset if needed
if (is.null(weights))
  weights <- rep.int(1, nobs)
if (is.null(offset))
  offset <- rep.int(0, nobs)

## get family functions:
variance <- family$variance
linkinv  <- family$linkinv
if (!is.function(variance) || !is.function(linkinv) )
  stop("'family' argument seems not to be a valid family object", call. = FALSE)
dev.resids <- family$dev.resids
aic <- family$aic
mu.eta <- family$mu.eta
unless.null <- function(x, if.null) if(is.null(x)) if.null else x
valideta <- unless.null(family$valideta, function(eta) TRUE)
validmu  <- unless.null(family$validmu,  function(mu) TRUE)
if(is.null(mustart)) {
  ## calculates mustart and may change y and weights and set n (!)
  eval(family$initialize)
} else {
  mukeep <- mustart
  eval(family$initialize)
  mustart <- mukeep
}
if(EMPTY) {
  eta <- rep.int(0, nobs) + offset
  if (!valideta(eta))
    stop("invalid linear predictor values in empty model", call. = FALSE)
  mu <- linkinv(eta)
  ## calculate initial deviance and coefficient
  if (!validmu(mu))
    stop("invalid fitted means in empty model", call. = FALSE)
  dev <- sum(dev.resids(y, mu, weights))
  w <- sqrt((weights * mu.eta(eta)^2)/variance(mu))
  residuals <- (y - mu)/mu.eta(eta)
  good <- rep_len(TRUE, length(residuals))
  boundary <- conv <- TRUE
  coef <- numeric()
  iter <- 0L
} else {
  coefold <- NULL
  eta <-
    if(!is.null(etastart)) etastart
  else if(!is.null(start))
    if (length(start) != nvars)
      stop(gettextf("length of 'start' should equal %d and correspond to initial coefs for %s", nvars, paste(deparse(xnames), collapse=", ")),
           domain = NA)
  else {
    coefold <- start
    offset + as.vector(if (NCOL(x) == 1L) x * start else x %*% start)
  }
  else family$linkfun(mustart)
  mu <- linkinv(eta)
  if (!(validmu(mu) && valideta(eta)))
    stop("cannot find valid starting values: please specify some", call. = FALSE)
  ## calculate initial deviance and coefficient
  devold <- sum(dev.resids(y, mu, weights))
  boundary <- conv <- FALSE

  ##------------- THE Iteratively Reweighting L.S. iteration -----------
  for (iter in 1L:control$maxit) {
    good <- weights > 0
    varmu <- variance(mu)[good]
    if (anyNA(varmu))
      stop("NAs in V(mu)")
    if (any(varmu == 0))
      stop("0s in V(mu)")
    mu.eta.val <- mu.eta(eta)
    if (any(is.na(mu.eta.val[good])))
      stop("NAs in d(mu)/d(eta)")
    ## drop observations for which w will be zero
    good <- (weights > 0) & (mu.eta.val != 0)

    if (all(!good)) {
      conv <- FALSE
      warning(gettextf("no observations informative at iteration %d",
                       iter), domain = NA)
      break
    }
    z <- (eta - offset)[good] + (y - mu)[good]/mu.eta.val[good]
    w <- sqrt((weights[good] * mu.eta.val[good]^2)/variance(mu)[good])
    # ## call Fortran code via C wrapper
    # fit <- .Call(C_Cdqrls, x[good, , drop = FALSE] * w, z * w,
    #              min(1e-7, control$epsilon/1000), check=FALSE)
    # 

    #print(iter)
    #print(z)
    #print(w)
  }


  }
  return(list(z=z, w=w, mustart=mustart, etastart=etastart, eta=eta, offset=offset, mu=mu, mu.eta.val=mu.eta.val,
              weight=weights, var_mu=variance(mu)))

}