Cómo programar una regresión logística en R desde 0

La regresión logística es uno de los algoritmos de machine learning más utilizados en los problemas de clasificación. En este post programaremos una regresión logística en R desde 0 para así aprender todo lo que hay que saber sobre este algoritmo. ¡Vamos a ello!

Entendiendo la regresión logística

La regresión logística tiene ciertas similitudes a la regresión lineal, la cual programamos desdeo 0 en R este post. Básicamente, la regresión lineal es una línea recta que para cada valor de x nos devuelve una predicción de nuestra variable y.

En el caso de la regresión logística, ya que también se trata de una regresión. Sin embargo, la regresión logística, al tratarse de un algoritmo de clasificación devuelve la probabilidad de pertenencia a una clase y. Además, en vez de ser una línea recta, la forma de la regresión logística es una S, como se ve en la imagen a continuación (más adelante explico por qué):

regresión logística vs regresión lineal

Así pues, mientras que la regresión lineal toma valores de -∞ a ∞, la regresión logística toma valores de 0 a 1 (el rango de la probabilidad).

Además, al igual que en el caso de la regresión lineal, la regresión logística puede haber dos tipos de modelos:

  • Un modelo simple, es decir, que solo se usa una variable predictora.
  • Un modelo complejo que incluye más de una variable predictora.

En cualquier caso, como la forma de la curva es diferente, es de esperar que el método para obtener esta curva también cambie. De hecho, mientras que la regresión lineal minimiza la suma de errores al cuadrado, la regresión logística busca la curva de máxima verosimilitud. ¡Veamos más sobre esto!

Buscando la curva de máxima verosimilitud

Al igual que en el caso de la regresión lineal, existen diferentes regresiones logísticas que podemos ajustar a nuestro modelo. Sin embargo, esto no significa que todas las regresiones sean igual de válidas. Veámos un ejemplo:

library(ggplot2)

x = c(1:10)
y = as.factor(c(rep(0,5), rep(1,5)))
data = data.frame(x,y)

sigmoid = function(x){
  1/(1+exp(-x))
}

data$function1 = sigmoid(data$x - 5) + 1
data$function2 = sigmoid(data$x -7) +1 

ggplot(data,aes(x)) + 
  geom_point(aes(y = y,col = y),size = 5) +
  geom_line(aes(y = function1), col = "green", linetype = "dashed") +
  geom_line(aes(y = function2), col = "blue", linetype = "dashed") +
  theme_minimal() + theme(legend.position = "bottom")
Ajuste de diferentes funciones logisticas

Como ves, podemos aplicar diferentes regresiones logísticas (las líneas azul y verdes). Si te fijas en la probabilidad que devuelve la regresión para el valor x > 5 verás que la regresión azul da mucha probabilidad de pertenencia al grupo 0. Por su parte, la regresión verde da probabilidades más altas. Viendo la gráfica podemos decir que la regresión verde se ajusta mejor que la azul, aunque esto no significa que sea la mejor regresión.

Básicamente, en el proceso lógico que he hecho en el párrafo anterior, inconscientemente, hemos calculado la máxima verosimilitud. Al fin y al cabo, la verosimilitud no es más que la probabilida de que, con los datos que tenemos, la verdadera regresión sea esa. La curva verde es más probable y, por tanto, más verosimil, y por eso la hemos elegido.

Así pues, la verosimilitud se obtiene multiplicando la probabilidad de observar cada una de las observaciones dada una regresión. Así pues, diferentes regresiones tendrán diferente verosimilitud, y nosotros nos quedaremos con aquella regresión que maximice la verosimilitud.

Nota: para facilitar la explicación matemática, supondremos dos cuestiones:

  • La función que sige la regresión logística es la función sigmoide, que representaremos así: \(σ(z)=\frac{1}{1+e^{-z}}\)
  • Tenedremos un vector de parámetros x y longitud m, que usaremos para calcular la verosimilitud: \(θ^T x = \sum^m_{i=1}θ_ix_i = θ_1x_1 + θ_2x_2+ ··· + θ_m x_m\)

Al tratarse de una predicción entre dos clases, podemos considerar que cada categoría sigue una distribución de Bernoulli, de tal forma que la probabilidad de pertenencia a la una clase es \(p = σ(θ^Tx)\).

Así pues, como la suma de proabilidades debe dar 1 (distribución de Bernoulli), con calcular la probabilidad de pertenencia a una clase nos serviría. Al fin y al cabo, la probabilidad de pertenencia a la otra clase sera lo que quede hasta llegar a 1 de la probabilidad de la primera clase. Así pues, la probabilidad de una observación se puede representar de la siguiente manera:

\[ P(Y = y|X = x) = σ(θ^Tx)^y · [1- σ(θ^Tx)]^{(1-y)} \]

Esta sería la función de verosimilitud y esta es la función que deberíamos maximizar. Sin embargo, el hecho de que sea una multiplicación de probabilidades hace que sea difícil encontrar el máximo.

Para solucionarlo, se toman logaritmos, de tal forma que en vez de maximizar la función de verosimilitud, se maximiza el logaritmo de la función de verosimilitud, también conocido como log likelihood.

\[ \ell(θ) = \sum_{i =1}^n y^ilog\;σ(θ^Tx^i) + (1-y^i)log\;σ(1-θ^Tx^i)\]

Ahora que tenemos esta función, debemos encontrar el valor de θ que maximiza la verosimilitud. Para ello no podríamos usar el método tradicional de derivar la función e igualarla a cero, ya que de este método no nos serviría (explicación).

En su lugar, deberemos utilizar un algoritmo de optimización que nos permita encontrar el óptimo de la función. Uno de los algoritmos de optimización típicos, utilizado también en otros algoritmos como las redes neuronales, es el ascenso del gradiente o gradient ascent.

Para poder utilizar este optimizador, deberíamos calcular la derivadas parciales de la función sobre θ, lo cual tras operar queda así:

\[ \frac{∂\ell(θ)}{∂θ_j}= [y-σ(θ^Tx)]x_j \]

Ahora que ya sabemos qué es lo que debemos optimizar, vamos a aplicar toda esa teoría para crear nuestra regresión logística en R desde 0. ¡Vamos allá!

Programar regresión logística en R desdeo 0

De cara a programar la regresión logística desde 0 en R, primero debemos entender cómo funcionan los algoritmos de optimización en R. Sobre ellos, hay que destacar dos cuestiones:

  • Los algoritmos de optimización en R no buscan máximos, sino mínimos. En nuestro caso buscamos máximos, por lo que de buenas a primeras no servirían. Para solucionar esto, simplemente tendremos que optimizar por el negativo del log likelihood.
  • Aunque podríamos programar la función del gradiente y pasársela al algoritmo, elpropio algoritmo es capaz de calcularnos la aproximación del gradiente. Por tanto, podemos ahorrarnos programar esta función.

Por otro lado, de cara a hacer la función lo más flexible posible, voy a considerar las siguientes cuestiones:

  1. Que el algoritmo funcione para cualquier dataset, sin tener que modificar nada internamente en la función.
  2. Que el algoritmo sea flexible y nos permita incluir, o no, el término alpha o intercept.

Teniendo esto en cuenta, para programar el negative log likelihood simplemente hay que:

  1. Calcular θT x, que como hemos visto anteriormente, simplemente se trata de la multiplicación de cada una de las variables por su parámetro y la suma del alfa, en caso de que lo haya. Esto habrá que hacerlo para todas las observaciones.
  2. Calculamos σθT , es decir, la probabilidad de cada una de las observaciones de pertenencia a una clase.
  3. Para cada una de las observaciones, calcular el logaritmo de theta, utilizando la fórmula que hemos mencionado anteriormente. A este resultado hay que incluirle el símbolo negativo (ya que en vez de maximizar, la función minimiza).

Así pues, nuestra función que devuelve el negativo del log likelihood es el siguiente:


sigmoid = function(x){
  1/(1+exp(-x))
}

neg_log_likelihood = function(par, data, y, include_alpha = T){
  x = data[,names(data) != y]
  y_data = data[,y]

  # 1. Calculamos theta
  if(include_alpha){

    # Multiplicamos cada valor por su parámetro
    multiplied =  mapply("*",x,par[2:length(par)])

    # Lo sumamos para cada observacion y sumamos el valor alpha
    theta =  rowSums(multiplied) + par[1]
  }else{
    theta =  rowSums(mapply("*",x,par))
  }

  # 2. Calculamos p
  p = sigmoid(theta)
  #p = exp(theta) / (1 + exp(theta))

  # 3. Calculamos el -log likelihood
  val = -sum(y_data * log(p) + (1-y_data)*log(1-p)) 

  return(val)
}

Ahora, simplemente debemos optimizar esta función. Para ello, en R viene por defecto la función optim. Aunque con esa función también sirva, usaré la función optimx del paquete optimx, ya que permite ver las estimaciones para diferentes modelos de optimización

Nota: como no existen mínimos locales en la función, podemos inicializar los parámetros con los valores que nosotros queramos (link a explicación).

library(optimx)

opt = optimx(
  par = c(0,0,0),
  fn = neg_log_likelihood,
  data = data,
  y = "obese",
  include_alpha = T,
  control = list(trace = 0, all.methods = TRUE)
)

opt[,1:5]
                       p1          p2        p3         value fevals
BFGS                   NA          NA        NA 8.988466e+307     NA
CG           2.792292e-03 -0.09950912 0.1877614  7.926620e+01    496
Nelder-Mead  3.271456e+01 -0.38368078 0.3679554  3.657498e+01    256
L-BFGS-B               NA          NA        NA 8.988466e+307     NA
nlm          3.265841e+01 -0.38340782 0.3680640  3.657492e+01     NA
nlminb       3.365225e+01 -0.38829738 0.3661810  3.664568e+01     70
spg                    NA          NA        NA 8.988466e+307     NA
ucminf                 NA          NA        NA 8.988466e+307     NA
Rcgmin      -2.075240e-24 -0.10142197 0.1913561  7.926331e+01   1013
Rvmmin       0.000000e+00  0.00000000 0.0000000  3.465736e+02      2
newuoa                 NA          NA        NA 8.988466e+307     NA
bobyqa                 NA          NA        NA 8.988466e+307     NA
nmkb                   NA          NA        NA 8.988466e+307     NA
hjkb                   NA          NA        NA 8.988466e+307     NA

Como vemos, diferentes modelos de optimización nos ofrecen diferentes valores para nuestros parámetros, aunque varios modelos (Nelder-Mead, nlm, nlm inb) parece que devuelven parámetros similares.

Comprobación de los valores de los parámetros del modelo

Aunque varios optimizadores hayan devuelto valores similares, eso no significa que esos parámetros estén bien, ya que puede que hayamos errado en algún momento a la hora de programar nuestra función.

Así pues, vamos a comprobar cuál es la estimación de los parámetros de la función base de R para calcular la regresión logística (por que sí, para usar esta función en R no hace falta programarla desde 0).

glm_result$coefficients
(Intercept)      Height      Weight 
 36.9508088  -0.4328929   0.4149990 

Como vemos las estimaciones obtenidas por nuestro algoritmo son muy similares a las obtenidas por el algoritmo base, lo cual es muy buen indicador.

Como todos los parámetros se parecen, nos quedaremos con los resultados del modelo nlminb .

unlist(opt["nlminb";,1:3])
        p1         p2         p3 
33.6522524 -0.3882974  0.3661810 

Ahora que ya sabemos que nuestro modelo es correcto.. ¡Veámos cómohacer las predicciones!

Cómo hacer predicciones regresión logística

Para hacer la predicción, simplemente tendremos que calcular la theta y multiplicarlo por la función sigmoide. De esta forma, obtendremos la probabilidad de cada una de las observaciones de obtener a una clase.

Así pues, vamos a programar esta función en R desde 0:


prediccion = function(x, par){

  # Hay alpha
  if(ncol(x) < length(par)){
    theta = rowSums(mapply("*",x,par[2:length(par)])) +  par[1]
  }else{
    theta = rowSums(mapply("*",x,par))
  }

  prob = sigmoid(theta)

  return(prob)

}

Ahora podemos calcular y ver cuál es la probabilidad de ser obeso de cada una de las observaciones y compararlo con el valor real mediante una matriz de confusión. Para ello, usaremos la función confusionMatrix del paquete caret.

confusionMatrix(as.factor(round(pediccions_obese)), as.factor(data$obese))
Confusion Matrix and Statistics

          Reference
Prediction   0   1
         0 168   6
         1   4 322

               Accuracy : 0.98            
                 95% CI : (0.9635, 0.9904)
    No Information Rate : 0.656           
    P-Value [Acc > NIR] : <2e-16          

                  Kappa : 0.9558          

 Mcnemar's Test P-Value : 0.7518          

            Sensitivity : 0.9767          
            Specificity : 0.9817          
         Pos Pred Value : 0.9655          
         Neg Pred Value : 0.9877          
             Prevalence : 0.3440          
         Detection Rate : 0.3360          
   Detection Prevalence : 0.3480          
      Balanced Accuracy : 0.9792          

       'Positive' Class : 0               
                                          

Como vemos, la regresión logística ha alcanzado un nivel de accuracy (sobre train) de un 98%, lo cual indica que efectivamente nuestro modelo ha funcionado correctamente.

Personalmente, una de las cuestiones que me gusta ver es cuál es el decission boundary de nuestro algoritmo y en qué casos está fallando:

library(ggplot2)
library(dplyr)

results =  cbind(data, prediction = round(pediccions_obese))

results %>%
  mutate(
    opacity = ifelse(prediction == obese,0.8,1)
  ) %>%
  ggplot(aes(Weight, Height, col = as.factor(obese),
             alpha = opacity)) +
  geom_point() +
  theme_minimal() +
  theme(legend.position = "bottom", panel.grid.minor = element_blank()) +
  scale_alpha_continuous(range =c(0.3,1), guide=FALSE)  +
  scale_x_continuous(breaks = seq(40,160,10)) +
  scale_y_continuous(breaks = seq(140,200,10)) 
Resultados Regresion Logistica R

Como vemos, las predicciones que han fallado, la mayoría lo han hecho por muy poco, excepto por dos casos: la persona que pesa más de 110kg, mide menos de 170cm y no es obesa la persona que pesa casi 80kg midiendo 150 y poco y tampoco es obesa.

Conclusión

Como vemos, la regresión logística es un algoritmo muy simple pero muy potente que tiene mucha más miga de lo que parece.

Sin duda alguna, programar un algoritmo desde 0 (en este caso la regresión logística en R), te ofrece tener un mayor conocimiento de qué es lo que hace el algoritmo, qué limitaciones tiene y cómo funciona.

Como siempre, espero que este post te haya gustado. Si es así, te animo a que te suscribas a la newsletter para que estés al tanto cuando suba nuevos posts. En cualquier caso, ¡nos vemos en el siguiente!

Blog patrocinado por: