Standardstartwerte, die die logistische Regression mit glm anpassen

Ich frage mich, wie die Standardstartwerte in angegeben werden glm.

Dieser Beitrag schlägt vor, dass Standardwerte als Nullen festgelegt werden. Das man sagt , dass es ein Algorithmus dahinter jedoch relevante Verbindung unterbrochen wird.

Ich habe versucht, ein einfaches logistisches Regressionsmodell mit einem Algorithmus-Trace zu versehen:

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

Erstens ohne Angabe von Anfangswerten:

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

Im ersten Schritt sind Anfangswerte NULL.

Zweitens setze ich Startwerte auf Nullen:

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

Und wir können sehen, dass sich die Iterationen zwischen dem ersten und dem zweiten Ansatz unterscheiden.

Um die von angegebenen Werte zu sehen, habe glmich versucht, das Modell mit nur einer Iteration anzupassen:

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

Schätzungen von Parametern entsprechen (nicht überraschend) Schätzungen des ersten Ansatzes in der zweiten Iteration, dh das [1] 0.386379 1.106234 Festlegen dieser Werte als Anfangswerte führt zu derselben Iterationssequenz wie im ersten Ansatz:

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

Die Frage ist also, wie diese Werte berechnet werden.

TL; DR

• start=c(b0,b1)initialisiert eta auf b0+x*b1(mu auf 1 / (1 + exp (-eta)))
• start=c(0,0) initialisiert eta auf 0 (mu auf 0,5) unabhängig vom y- oder x-Wert.
• start=NULL initialisiert eta = 1,098612 (mu = 0,75), wenn y = 1 ist, unabhängig vom x-Wert.

• start=NULL initialisiert eta = -1,098612 (mu = 0,25), wenn y = 0 ist, unabhängig vom x-Wert.

• Sobald eta (und folglich mu und var (mu)) berechnet wurde wund zim Geiste von berechnet und an einen QR-Löser gesendet wird qr.solve(cbind(1,x) * w, z*w).

Lange Form

Aufbauend auf Rolands Kommentar: Ich habe einen gemacht glm.fit.truncated(), wo ich glm.fitden C_CdqrlsAnruf angenommen und ihn dann auskommentiert habe. glm.fit.truncatedgibt die zund w-Werte (sowie die Werte der zur Berechnung von zund verwendeten Mengen w) aus, die dann an den C_CdqrlsAufruf übergeben werden:

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

Mehr dazu C_Cdqrls hier . Glücklicherweise greift die Funktion qr.solvein Basis R direkt auf die LINPACK-Versionen zu, die in aufgerufen werden glm.fit().

Wir laufen also glm.fit.truncatedfür die verschiedenen Startwertspezifikationen und rufen dann qr.solvemit den w- und z-Werten auf, und wir sehen, wie die "Startwerte" (oder die ersten angezeigten Iterationswerte) berechnet werden. Wie Roland angedeutet hat, wirkt sich die Angabe von start=NULLoder start=c(0,0)in glm () auf die Berechnungen für w und z aus, nicht für start.

Für den Start = NULL: zist ein Vektor, bei dem die Elemente den Wert 2.431946 oder -2.431946 haben, und wist ein Vektor, bei dem alle Elemente 0.4330127 sind:

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 

Für den Start = c (0,0): zist ein Vektor, bei dem die Elemente den Wert 2 oder -2 haben, und wist ein Vektor, bei dem alle Elemente 0,5 sind:

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

Das ist alles schön und gut, aber wie berechnen wir das wund z? Nahe dem Boden sehen glm.fit.truncated()wir

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

Sehen Sie sich die folgenden Vergleiche zwischen den ausgegebenen Werten der zur Berechnung verwendeten Mengen an zund 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)

Es ist zu beachten, dass start.is.00der Vektor munur die Werte 0,5 hat, da eta auf 0 gesetzt ist und mu (eta) = 1 / (1 + exp (-0)) = 0,5. start.is.nullsetzt diejenigen mit y = 1 auf mu = 0,75 (was eta = 1,098612 entspricht) und diejenigen mit y = 0 auf mu = 0,25 (was eta = -1,098612 entspricht) und damit var_mu= 0,75 * 0,25 = 0,1875.

Es ist jedoch interessant festzustellen, dass ich den Samen geändert und alles neu interpretiert habe und mu = 0,75 für y = 1 und mu = 0,25 für y = 0 (und somit die anderen Mengen gleich geblieben sind). Das heißt, start = NULL führt zu demselben wund zunabhängig davon, was yund xsind, weil sie eta = 1,098612 (mu = 0,75) initialisieren, wenn y = 1 und eta = -1,098612 (mu = 0,25), wenn y = 0.

Es scheint also, dass ein Startwert für den Intercept-Koeffizienten und für den X-Koeffizienten nicht für start = NULL festgelegt wird, sondern dass eta abhängig vom y-Wert und unabhängig vom x-Wert Anfangswerte erhalten. Von dort aus wund zberechnet wird , dann zusammen mit geschickt xan der qr.solver.

Code, der vor den obigen Chunks ausgeführt werden soll:

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

}