Codificación de la función de verosimilitud para regresión logística

Question

Agradecería ayuda para entender si hice una interpretación y codificación correcta de la función de verosimilitud para la regresión logística.

Antecedentes: Para una tarea voy a escribir una función en R de la función de verosimilitud para un modelo de regresión logística. Se ha dado que $p(x_i)=\frac{1}{1+\exp(-x_i\beta)}$ donde $\beta=(\beta_1,...,\beta_k)$ y $x_i=(x_{i,1},...,x_{i,k})$ y que la matriz de diseño $\mathbf{x}$ es una matriz $N \times k$.

Mi entendimiento: Investigué y encontré que la función de verosimilitud para una regresión logística está dada por $L(\beta)=\prod p_i(\beta)^{y_i}(1-p_i(\beta))^{1-y_i}$ por lo que escribí el siguiente código para esta fórmula:

L <- function(beta,y,X){
  p <- numeric()
  for(i in 1:nrow(X)) {             
    p <- c(p, 1/(1+exp(as.vector(-X[i,])%*%beta)))
  }
  return(prod((p^y)%*%(1-p)^(1-y)))
  }

Primero obtuve todos los valores de p de las filas en la matriz X con la ayuda de un bucle for, y los guardé en p. Esto me permitió luego insertar todos los vectores en la función de verosimilitud dada anteriormente. Observar que $\mathbf{y}=(y_1,...,y_N)^T$ y que $\mathbf{p}=(p_1,...,p_N)^T.$

Pregunta: No tengo una forma de verificar si mi función es correcta, y he basado mi código para la verosimilitud en la regresión logística en la literatura del curso, por lo que agradecería mucho ayuda con respecto a si mi interpretación es correcta. ¿Es verdad que una regresión logística tiene la función de verosimilitud dada? ¿Y es mi código una forma correcta de implementar esto?

Answer 1

Necesitas hacer cálculos de probabilidad de esta manera en el espacio logarítmico

Tu cálculo está utilizando un método ingenuo que no es cómo codificamos funciones de probabilidad (u otras funciones de probabilidad) en la práctica. Al realizar cálculos con funciones de probabilidad, es mucho mejor trabajar en el espacio logarítmico para evitar problemas de subdesbordamiento aritmético y problemas con la inestabilidad numérica. Esto requiere que trabajes exclusivamente con log-probabilidades en todos tus pasos intermedios y luego transformes de nuevo para obtener la función de probabilidad solo al final de tu cálculo.

Para hacer esto, necesitas calcular la función de log-probabilidad utilizando log-probabilidades en todos los cálculos intermedios. La función de log-probabilidad para la regresión logística es:

$$\ell_\mathbf{\mathbf{x},\mathbf{y}}(\boldsymbol{\beta}) = \sum_{i=1}^n [y_i \log(p_i(\boldsymbol{\beta})) + (1-y_i) \log(1-p_i(\boldsymbol{\beta}))],$$

donde podemos escribir las log-probabilidades en esta expresión (usando la función hiperbólica log1pexp) como:

$$\begin{align} \log(p_i(\boldsymbol{\beta})) &= \log \Bigg( \frac{1}{1 + \exp(-\mathbf{x}_i \boldsymbol{\beta})} \Bigg) \\[6pt] &= -\log \Bigg( 1 + \exp(-\mathbf{x}_i \boldsymbol{\beta}) \Bigg) \\[12pt] &= -\text{log1pexp} (- \mathbf{x}_i \boldsymbol{\beta}), \\[12pt] \log(1-p_i(\boldsymbol{\beta})) &= \log \Bigg( 1 - \frac{1}{1 + \exp(-\mathbf{x}_i \boldsymbol{\beta})} \Bigg) \\[6pt] &= \log \Bigg( \frac{\exp(-\mathbf{x}_i \boldsymbol{\beta})}{1 + \exp(-\mathbf{x}_i \boldsymbol{\beta})} \Bigg) \\[6pt] &= -\mathbf{x}_i \boldsymbol{\beta} - \log \Bigg( 1 + \exp(-\mathbf{x}_i \boldsymbol{\beta}) \Bigg) \\[12pt] &= -\mathbf{x}_i \boldsymbol{\beta} -\text{log1pexp} (- \mathbf{x}_i \boldsymbol{\beta}). \\[6pt] \end{align}$$

Esto nos permite escribir todo el modelo en el espacio logarítmico como:

$$\ell_\mathbf{\mathbf{x},\mathbf{y}}(\boldsymbol{\beta}) = - \sum_{i=1}^n [\text{log1pexp} (- \mathbf{x}_i \boldsymbol{\beta}) + (1-y_i) \mathbf{x}_i \boldsymbol{\beta}].$$

A continuación he creado una función de R que calcula la función de probabilidad para la regresión logística utilizando el cálculo en el espacio logarítmico. Este método de programación de la función será mucho mejor que lo que tienes en tu pregunta, ya que evitará el subdesbordamiento aritmético que ocurre al multiplicar muchas pequeñas probabilidades. En esta función he incluido una opción lógica loglike para permitir al usuario especificar si quieren devolver la función de probabilidad o la función de log-probabilidad; observa que incluso si el usuario solicita la función de probabilidad, toda la computación intermedia se realiza en el espacio logarítmico y solo se transforma de nuevo a la función de probabilidad en el último paso.

likelihood.logistic.regression <- function(beta, y, x, loglike = FALSE) {

  #Obtener parámetros
  n <- length(y)

  #Calcular la log-probabilidad
  TERMS <- rep(-Inf, n)
  for (i in 1:n) {
    XB <- sum(x[i, ]*beta)
    TERMS[i] <- - VGAM::log1pexp(-XB) - (1-y[i])*XB }
  LOGLIKE <- sum(TERMS)

  #Retornar la salida
  if (loglike) { LOGLIKE } else { exp(LOGLIKE) } }

Como se señala en los comentarios, una forma de verificar el código es simular datos de una regresión logística y luego usar la función de log-probabilidad en un optimizador para calcular el MLE del vector de parámetros. Si tu función está escrita correctamente (y si todas las demás funciones que estás usando en tu análisis están escritas correctamente) entonces el MLE que calculas a partir de una gran simulación debería ser razonablemente cercano a los valores de parámetro reales utilizados para generar la simulación. Dejaré esto como un ejercicio para ti.