Domyślne wartości początkowe pasujące do regresji logistycznej z glm

Zastanawiam się, jak określono domyślne wartości początkowe w glm.

Ten post sugeruje, że wartości domyślne są ustawione na zera. Ten jeden mówi, że istnieje algorytm za nim, jednak istotne link jest uszkodzony.

Próbowałem dopasować prosty model regresji logistycznej ze śledzeniem algorytmu:

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)))

Po pierwsze, bez specyfikacji wartości początkowych:

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

W pierwszym kroku początkowe wartości to NULL.

Po drugie, ustawiam wartości początkowe na zero:

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

Widzimy, że iteracje między pierwszym i drugim podejściem są różne.

Aby zobaczyć wartości początkowe określone przez glm, próbowałem dopasować model tylko z jedną iteracją:

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

Oszacowania parametrów (co nie jest zaskakujące) odpowiadają oszacowaniom pierwszego podejścia w drugiej iteracji, tj. [1] 0.386379 1.106234 Ustawienie tych wartości jako wartości początkowych prowadzi do tej samej sekwencji iteracji, co w pierwszym podejściu:

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

Pytanie brzmi: jak te wartości są obliczane?

TL; DR

• start=c(b0,b1)inicjuje eta do b0+x*b1(mu do 1 / (1 + exp (-eta)))
• start=c(0,0) inicjuje eta do 0 (mu do 0,5) niezależnie od wartości y lub x.
• start=NULL inicjuje eta = 1,098612 (mu = 0,75), jeśli y = 1, niezależnie od wartości x.

• start=NULL inicjuje eta = -1,098612 (mu = 0,25), jeśli y = 0, niezależnie od wartości x.

• Po ETA (i co za tym idzie mu i var (il)), została obliczona wi zsą obliczane i wysyłane do QR Solver w duchu qr.solve(cbind(1,x) * w, z*w).

Długa forma

Budowanie off komentarz Rolanda: Zrobiłem glm.fit.truncated(), gdzie wziąłem glm.fitdół do C_Cdqrlsrozmowy, a następnie zauważył go. glm.fit.truncatedwyprowadza wartości zi w(a także wartości ilości użytych do obliczeniaz i w), które następnie zostaną przekazane do C_Cdqrlswywołania:

## 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) 

Więcej można przeczytać o C_Cdqrls tutaj . Na szczęście funkcja qr.solvew podstawowym R jest podłączana bezpośrednio do wywoływanych wersji LINPACKglm.fit() .

Sprawdzamy więc glm.fit.truncatedróżne specyfikacje wartości początkowej, a następnie wywołujemy qr.solvez wartościami w i z, i widzimy, jak obliczane są „wartości początkowe” (lub pierwsze wyświetlane wartości iteracji). Jak wskazał Roland, określenie start=NULLlub start=c(0,0)w glm () wpływa na obliczenia w i z, a nie dla start.

Na początek = NULL: zjest wektorem, w którym elementy mają wartość 2,431946 lub -2,431946 i wjest wektorem, w którym wszystkie elementy mają 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 

Dla początku = c (0,0): zjest wektorem, w którym elementy mają wartość 2 lub -2 i wjest wektorem, w którym wszystkie elementy mają 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 

To wszystko dobrze i dobrze, ale jak obliczyć wiz ? W dolnej części glm.fit.truncated()widzimy

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

Spójrz na następujące porównania między wartościami wyjściowymi ilości użytych do obliczenia z i 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)

Zauważ, że start.is.00 wektor będzie miał mutylko wartości 0,5, ponieważ eta jest ustawiona na 0, a mu (eta) = 1 / (1 + exp (-0)) = 0,5. start.is.nullustawia te zy = 1 na mu = 0,75 (co odpowiada eta = 1,098612), a te zy = 0 na mu = 0,25 (co odpowiada eta = -1,098612), a zatem var_mu= 0,75 * 0,25 = 0,1875.

Warto jednak zauważyć, że zmieniłem ziarno i przestawiłem wszystko, a mu = 0,75 dla y = 1 i mu = 0,25 dla y = 0 (a zatem pozostałe wielkości pozostały takie same). To znaczy, start = NULL daje to samow i zbez względu na to, co yi xsą, ponieważ inicjują one eta = 1,098612 (mu = 0,75), jeśli y = 1 i eta = -1,098612 (mu = 0,25), jeśli y = 0.

Wygląda więc na to, że wartość początkowa dla współczynnika przechwytywania i dla współczynnika X nie jest ustawiona dla start = NULL, ale raczej wartości początkowe są podawane eta w zależności od wartości y i niezależne od wartości x. Stamtąd wi zsą obliczane, a następnie przesyłane wraz z xdo qr.solver.

Kod uruchamiany przed powyższymi fragmentami:

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)))

}