Valeurs de départ par défaut correspondant à la régression logistique avec glm

Je me demande comment les valeurs de départ par défaut sont spécifiées dans glm.

Ce message suggère que les valeurs par défaut sont définies comme des zéros. Celui- ci dit qu'il y a un algorithme derrière, cependant le lien pertinent est rompu.

J'ai essayé d'adapter un modèle de régression logistique simple avec une trace d'algorithme:

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

Tout d'abord, sans spécification des valeurs initiales:

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

Dans la première étape, les valeurs initiales sont NULL.

Deuxièmement, j'ai défini des valeurs de départ comme des zéros:

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

Et nous pouvons voir que les itérations entre la première et la seconde approche diffèrent.

Pour voir les valeurs initiales spécifiées par glmJ'ai essayé d'ajuster le modèle avec une seule itération:

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

Les estimations des paramètres (sans surprise) correspondent aux estimations de la première approche dans la deuxième itération, c'est-à-dire que la [1] 0.386379 1.106234 définition de ces valeurs comme valeurs initiales conduit à la même séquence d'itérations que dans la première approche:

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

La question est donc de savoir comment ces valeurs sont calculées?

TL; DR

• start=c(b0,b1)initialise eta à b0+x*b1(mu à 1 / (1 + exp (-eta)))
• start=c(0,0) initialise eta à 0 (mu à 0,5) quelle que soit la valeur y ou x.
• start=NULL initialise eta = 1,098612 (mu = 0,75) si y = 1, quelle que soit la valeur x.

• start=NULL initialise eta = -1,098612 (mu = 0,25) si y = 0, quelle que soit la valeur x.

• Une fois que eta (et par conséquent mu et var (mu)) a été calculé, wet zest calculé et envoyé à un solveur QR, dans l'esprit de qr.solve(cbind(1,x) * w, z*w).

Forme longue

S'appuyant sur le commentaire de Roland: J'ai fait un glm.fit.truncated(), où j'ai pris glm.fitl' C_Cdqrlsappel, puis je l'ai commenté. glm.fit.truncatedrenvoie les valeurs zet w(ainsi que les valeurs des quantités utilisées pour calculer zet w) qui seraient ensuite transmises à l' C_Cdqrlsappel:

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

Pour en savoir plus, C_Cdqrls cliquez ici . Heureusement, la fonction qr.solvede base R puise directement dans les versions de LINPACK sollicitées dans glm.fit().

Nous courons donc glm.fit.truncatedpour les différentes spécifications de valeur de départ, puis faisons un appel à qr.solveavec les valeurs w et z, et nous voyons comment les "valeurs de départ" (ou les premières valeurs d'itération affichées) sont calculées. Comme Roland l'a indiqué, spécifier start=NULLou start=c(0,0)dans glm () affecte les calculs pour w et z, pas pour start.

Pour le début = NULL: zest un vecteur où les éléments ont la valeur 2,431946 ou -2,431946 et west un vecteur où tous les éléments sont 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 

Pour le début = c (0,0): zest un vecteur où les éléments ont la valeur 2 ou -2 et west un vecteur où tous les éléments sont 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 

C'est bien beau, mais comment calculer le wet z? Près du bas de glm.fit.truncated()nous voyons

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

Regardez les comparaisons suivantes entre les valeurs produites des quantités utilisées pour calculer zet 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)

Notez que start.is.00le vecteur mun'aura que les valeurs 0,5, car eta est défini sur 0 et mu (eta) = 1 / (1 + exp (-0)) = 0,5. start.is.nulldéfinit ceux avec y = 1 à mu = 0,75 (ce qui correspond à eta = 1,098612) et ceux avec y = 0 à mu = 0,25 (ce qui correspond à eta = -1,098612), et donc le var_mu= 0,75 * 0,25 = 0,1875.

Cependant, il est intéressant de noter que j'ai tout changé et réanimé et le mu = 0,75 pour y = 1 et mu = 0,25 pour y = 0 (et donc les autres quantités sont restées les mêmes). C'est-à-dire que start = NULL donne la même chose wet zindépendamment de ce yque xsont et parce qu'ils initialisent eta = 1.098612 (mu = 0.75) si y = 1 et eta = -1.098612 (mu = 0.25) si y = 0.

Il apparaît donc qu'une valeur de départ pour le coefficient d'interception et pour le coefficient X n'est pas définie pour start = NULL, mais plutôt des valeurs initiales sont données à eta en fonction de la valeur y et indépendamment de la valeur x. De là wet zsont calculés, puis envoyés avec xle qr.solver.

Code à exécuter avant les morceaux ci-dessus:

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

}