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 bCodificación parcial y python
2019
ciencia de datos
estadística
R
python
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
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 1Y 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 1Si 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-16Y 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: int64convertimos 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 0x_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.