Codificación parcial y python

2019
ciencia de datos
estadística
R
python
Author

José Luis Cañadas Reche

Published

July 15, 2019

O como se conoce en estos tiempos modernos one hot encoding. En realidad se trata simplemente de cómo codificar una variable categórica en un conjunto de números que un algoritmo pueda utilizar.

Ya hablé de esto mismo en el post codificación de variables categóricas I

Básicamente, la codificación parcial lo que hace es crearse tantas variables indicadoras como niveles tengo en mi variable menos 1.

Ejemplo. Construimos un conjunto de datos simple, con 3 variables

set.seed(155)

x1 <- rnorm(n = 100, mean = 4, sd = 1.5 )
x2_cat <- factor(rep(c("a","b","c","d"),  25 ))
y <- numeric(length = 100)

# Construimos artificialmente relación entre x e y
y <- 2 +  4 * x1 + rnorm(25, 0, 1)

# cambiamos "el intercept" según la variable categórica

y[x2_cat == "a"] <- y[x2_cat == "a"] + 8
y[x2_cat == "b"] <- y[x2_cat == "b"] - 5
y[x2_cat == "c"] <- y[x2_cat == "c"] + 3
y[x2_cat == "d"] <- y[x2_cat == "d"] - 3

dat <- data.frame(y, x1, x2_cat)
head(dat)
#>          y       x1 x2_cat
#> 1 31.00555 5.200100      a
#> 2 19.13571 5.061407      b
#> 3 20.49049 3.888438      c
#> 4 17.93157 4.978832      d
#> 5 25.23501 3.989047      a
#> 6 21.86234 6.272138      b

En R al definir x2_cat como un factor él ya sabe que para ciertos métodos (por ejemplo una regresión) hay que codificar esa variable y por defecto utiliza la codificación parcial. Con la función contrasts vemos como lo hace.

contrasts(dat$x2_cat)
#>   b c d
#> a 0 0 0
#> b 1 0 0
#> c 0 1 0
#> d 0 0 1

Y en las columnas tenemos las 3 variables indicadoras que ha construido. ¿Por qué 3? pues muy fácil, puesto que la “a” se puede codificar con el valor 0 en las tres variables indicadoras. Esto que parece una obviedad evita problemas de colinealidad en los algoritmos de regresión por ejemplo.

fit1 <-  lm(y ~ x1 + x2_cat, data = dat)
summary(fit1)
#> 
#> Call:
#> lm(formula = y ~ x1 + x2_cat, data = dat)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -1.55549 -0.67355 -0.00943  0.59814  1.99480 
#> 
#> Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  10.12924    0.28393   35.67   <2e-16 ***
#> x1            3.94536    0.06034   65.39   <2e-16 ***
#> x2_catb     -12.95740    0.26148  -49.55   <2e-16 ***
#> x2_catc      -4.97712    0.25845  -19.26   <2e-16 ***
#> x2_catd     -10.98359    0.25785  -42.60   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.9094 on 95 degrees of freedom
#> Multiple R-squared:  0.9855, Adjusted R-squared:  0.9849 
#> F-statistic:  1616 on 4 and 95 DF,  p-value: < 2.2e-16

¿Qué hubiera pasado si hubiéramos tratado con 4 variables indicadoras?

dat2 <- dat
dat2$ind1 <- ifelse(dat$x2_cat == "a", 1, 0)
dat2$ind2 <- ifelse(dat$x2_cat == "b", 1, 0)
dat2$ind3 <- ifelse(dat$x2_cat == "c", 1, 0)
dat2$ind4 <- ifelse(dat$x2_cat == "d", 1, 0)
head(dat2[dat2$x2_cat=="d", ],3)
#>           y       x1 x2_cat ind1 ind2 ind3 ind4
#> 4  17.93157 4.978832      d    0    0    0    1
#> 8   9.30838 2.736958      d    0    0    0    1
#> 12 12.31765 2.943479      d    0    0    0    1

Si metemos ahora esas variables en el modelo

fit2 <-  lm(y ~ x1 +  ind2 + ind3 + ind4 + ind1, data = dat2)
summary(fit2)
#> 
#> Call:
#> lm(formula = y ~ x1 + ind2 + ind3 + ind4 + ind1, data = dat2)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -1.55549 -0.67355 -0.00943  0.59814  1.99480 
#> 
#> Coefficients: (1 not defined because of singularities)
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  10.12924    0.28393   35.67   <2e-16 ***
#> x1            3.94536    0.06034   65.39   <2e-16 ***
#> ind2        -12.95740    0.26148  -49.55   <2e-16 ***
#> ind3         -4.97712    0.25845  -19.26   <2e-16 ***
#> ind4        -10.98359    0.25785  -42.60   <2e-16 ***
#> ind1               NA         NA      NA       NA    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.9094 on 95 degrees of freedom
#> Multiple R-squared:  0.9855, Adjusted R-squared:  0.9849 
#> F-statistic:  1616 on 4 and 95 DF,  p-value: < 2.2e-16

Y vemos que como hay colinealidad R no estima el coeficiente de una de las variables indicadoras y hasta nos avisa con el mensaje Coefficients: (1 not defined because of singularities)

Pues la verdad es que mola que R sepa como tratar las categóricas si las has definido como factor pero también hace que la gente se olvide de que lo que en realidad hace es la codificación parcial.

Hablando de esto con un colega salió a colación que en python hay que explicitar la codificación y que quizá eso sea bueno porque así se sabe lo que se está haciendo y no hay lugar a dudas. Hasta aquí todo correcto, salvo que leyendo la documentación de pandas get_dummies resulta que por defecto construye tantas variables indicadoras como categorías y sólo tiene como opcional lo de quitar la primera con el parámetro drop_first, total me dije, no pasa nada, veamos como lo hace scikit learn y nada, resulta que por defecto también deja todas OneHotEncoder.

Reflexionando me dije, bueno, pues entonces cuando haga una regresión lineal con sklearn si uso las opciones por defecto de codificar las categóricas pues me debe saltar lo mismo que en R, es decir que hay un coeficiente que no puede estimar, pero resulta que sklearn hace un pelín de trampa y no salta el error, y no salta porque en sklearn la regresión lineal no ajusta una regresión lineal clásica, sino que por defecto y sin que tú lo pidas te hace una regresión regularizada y entonces no salta ese problema.

Pues la verdad , ¿qué puedo decir? no me hace gracia que por defecto no me quite la variable indicadora que sobra ni que haga regresión con regularización sin yo decirle nada.

En fin, veamos el ejemplo con python, aprovecho que escribo en un rmarkdown y puedo pasar objetos de R a python entre chunks sin muchos problemas.

import pandas as pd
import numpy as np
from sklearn.linear_model import LinearRegression
dat_py = r.dat
dat_py.describe()
#>                 y          x1
#> count  100.000000  100.000000
#> mean    18.633855    3.988010
#> std      7.402175    1.540500
#> min     -1.163152    0.315987
#> 25%     14.283915    2.815231
#> 50%     19.325282    3.885684
#> 75%     22.438970    5.062777
#> max     43.075931    8.304381
dat_py.x2_cat.value_counts()
#> a    25
#> b    25
#> c    25
#> d    25
#> Name: x2_cat, dtype: int64

convertimos a dummies con pandas por ejemplo

dat_py = pd.get_dummies(data=dat_py)
print(dat_py.head())
#>            y        x1  x2_cat_a  x2_cat_b  x2_cat_c  x2_cat_d
#> 0  31.005546  5.200100         1         0         0         0
#> 1  19.135715  5.061407         0         1         0         0
#> 2  20.490494  3.888438         0         0         1         0
#> 3  17.931571  4.978832         0         0         0         1
#> 4  25.235006  3.989047         1         0         0         0
x_variables = ['x1', 'x2_cat_a', 'x2_cat_b','x2_cat_c','x2_cat_d']
# Selecciono y convierto a numpy array
X = dat_py[x_variables].values  
y = dat_py['y'].values
X[0:3]
#> array([[5.20010037, 1.        , 0.        , 0.        , 0.        ],
#>        [5.06140731, 0.        , 1.        , 0.        , 0.        ],
#>        [3.88843797, 0.        , 0.        , 1.        , 0.        ]])
y[0:3]
#> array([31.00554596, 19.1357146 , 20.49049385])

lm = LinearRegression()
fit_python = lm.fit(X,y)
print('Intercept: ',fit_python.intercept_)
#> Intercept:  2.899716060541241
print('Coef: ',fit_python.coef_)
#> Coef:  [ 3.94536095  7.22952708 -5.72787734  2.25241099 -3.75406074]

Y vemos que si estima todos los coeficientes cuando no debería haber podido, esto tiene que ver como he dicho antes con que LinearRegression de sklearn no es la regresión lineal al uso sino que mete regularización.

Otro día veremos la librería statmodels de python cuya salida nos da una información más rica de los modelos y bastante parecida a lo que estamos acostumbrados con R.

Nota: Leyendo la docu de LinearRegression en ningún sitio dice que use regularización así que no alcanzo a entender por qué ha podido estimar todos los coeficientes. A ver si alguno de mis amigos pythonisos me lo aclara.