logo

Transformaciones

En muchas ocasiones cuando se trabaja con datos reales, no se logra tener el supuesto de normalidad. Uns manera sencilla de solucionar este problema, es aplicar alguna transformación a las variables originales. Estas transformaciones pueden ayudar a que las variables se aproximen a la distribución normal.

Exite una lista muy larga transformaciones, a continuación, presentamos las más usadas en la literatura. Es importante resaltar que se debe siempre realizar una validación adecuada para asegurar que las transformaciones usadas fueron apropiadas.

Lista de transformaciones

  1. Logarítmica: \[y=\log_b(x+a)\]
  2. Radical: \[y=\sqrt[b]{x+a}\]
  3. Inverso: \[y=\dfrac{1}{x}\]
  4. Potencia: \[y=x^a\]
  5. Logit: \[y=\dfrac{x}{1-x}\]
  6. Fisher: \[y=\dfrac{1+x}{1-x}\]
  7. Exponencial: \[y=e^{x}\]
  8. Trigonométricas inversas: \[\arcsin(x), \,\, \arccos(x) \]
  9. Box-Cox: \[ y = \begin{cases} \dfrac{x^{\lambda}-1}{\lambda} & \lambda \not =0 \\ & \\ \log(x) & \lambda = 0 \end{cases} \]
  10. Yeo-Johson: \[ y(x) = \begin{cases} \frac{(x+1)^{\lambda} - 1}{\lambda} & \text{si } x \geq 0 \text{ y } \lambda \neq 0\\ \log(x+1) & \text{si } x \geq 0 \text{ y } \lambda = 0\\ -\frac{(-x+1)^{2-\lambda} - 1}{2-\lambda} & \text{si } x < 0 \text{ y } \lambda \neq 2\\ -\log(-x+1) & \text{si } x < 0 \text{ y } \lambda = 2 \end{cases} \]

Ejemplo 1:

Para este primer ejemplo, simularemos los datos y veremos que tan conveniente nos resulta aplicar algunas transformaciones.

n  <- 100
set.seed(2023)
x1 <- rgamma(n,5,1)
x2 <- rbeta(n,1,1)

# Histogramas
par(mfrow=c(1,2))
hist(x1,probability=T,main=" ",xlab="Datos gamma",ylab=" ")
hist(x2,probability=T,main=" ",xlab="Datos beta",ylab=" ")

# Pruebas gráficas
par(mfrow=c(1,2))
qqnorm(x1,main=" ",xlab="Cuantiles teoricos",ylab="Cuantiles muestrales")
qqline(x1,col="blue",lwd=2)
qqnorm(x2,main=" ",xlab="Cuantiles teoricos",ylab="")
qqline(x2,col="blue",lwd=2)

# Pruebas analíticas
SW1 <- shapiro.test(x1)
SW1
## 
##  Shapiro-Wilk normality test
## 
## data:  x1
## W = 0.9226, p-value = 1.94e-05
SW2 <- shapiro.test(x2)
SW2
## 
##  Shapiro-Wilk normality test
## 
## data:  x2
## W = 0.95497, p-value = 0.00179

Observamos claramente por medio del histograma la falta de normalidad y luego lo podemos comprobar de manera gráfica y analítica.

# Transformaciones sugeridas
x1_trans  <- log(x1)        
x2_trans  <- acos(sqrt(x2))  

# Histogramas
par(mfrow=c(1,2))
hist(x1_trans,probability=T,main=" ",xlab="Datos gamma_trans",ylab= " ")
x <- seq(min(x1_trans)-1, max(x1_trans)+1, length = 100)
y <- dnorm(x, mean(x1_trans),sd(x1_trans))
lines(x, y, col = "red",lwd=3)
hist(x2_trans,probability=T,main=" ",xlab="Datos beta_trans",ylab= " ")
x <- seq(min(x2_trans)-1, max(x2_trans)+1, length = 100)
y <- dnorm(x, mean(x2_trans),sd(x2_trans))
lines(x, y, col = "red",lwd=3)

# Pruebas gráficas
par(mfrow=c(1,2))
qqnorm(x1_trans,main=" ",xlab="Cuantiles teoricos",ylab="Cuantiles muestrales")
qqline(x1_trans,col="blue",lwd=2)
qqnorm(x2_trans,main=" ",xlab="Cuantiles teoricos",ylab=" ")
qqline(x2_trans,col="blue",lwd=2)

# Pruebas analíticas
library(nortest)
SW1 <- shapiro.test(x1_trans)
SW1
## 
##  Shapiro-Wilk normality test
## 
## data:  x1_trans
## W = 0.9886, p-value = 0.5536
SW2 <- shapiro.test(x2_trans)
SW2
## 
##  Shapiro-Wilk normality test
## 
## data:  x2_trans
## W = 0.98059, p-value = 0.1478

Logramos conseguir la normalidad en las dos variables, pero el problema estaría en como logramos encontrar cuales sugerir las transformaciones adecuadas.

Ejemplo 2:

Una familia de funciones importantes para transformar es Box-Cox. Esta transformación depende del parámetro \(\lambda\). El valor de \(\lambda\) se estimar a partir de los datos y su estimación nos proporciona la transformación más adecuada que podamos aplicar. Para valores exactos de \(\lambda\) tenemos:

  1. Si \(\lambda = 0\), utilice: \[y=\log(x)\]

  2. Si \(\lambda=\dfrac{1}{2}\), utilice: \[y=\sqrt{x}\]

  3. Si \(\lambda=-\dfrac{1}{2}\), utilice: \[y=\dfrac{1}{\sqrt{x}}\]

  4. Si \(\lambda= 2\), utilice: \[y=x^2\]

  5. Si \(\lambda= -1\), utilice: \[y=\dfrac{1}{{x}}\]

  6. Para otros casos, ingrese el valor estimado de \(\lambda\). Este valor tiene un rango de −2 y 2.

Para probar estas transformaciones, simularemos variables, observamos los valores de \(\lambda\) estimados y realizaremos las transformaciones sugeridas. Es importante resaltar que estas transformaciones solo se pueden aplicar a datos positivos.

n   <- 100
set.seed(2023)
x1  <- rexp(n,1)
x2  <- rf(n,10,1)

# Histogramas
par(mfrow=c(1,2))
hist(x1,probability=T,main=" ",xlab="Datos exp",ylab="Densidad")
hist(x2,probability=T,main=" ",xlab="Datos F",ylab="Densidad")

# Pruebas gráficas
par(mfrow=c(1,2))
qqnorm(x1,main=" ",xlab="Cuantiles teóricos",ylab="Cuantiles muestrales")
qqline(x1,col="blue",lwd=2)
qqnorm(x2,main=" ",xlab="Cuantiles teóricos",ylab=" ")
qqline(x2,col="blue",lwd=2)

# Pruebas analíticas
SW1 <- shapiro.test(x1)
SW1
## 
##  Shapiro-Wilk normality test
## 
## data:  x1
## W = 0.85056, p-value = 1.244e-08
SW2 <- shapiro.test(x2)
SW2
## 
##  Shapiro-Wilk normality test
## 
## data:  x2
## W = 0.16532, p-value < 2.2e-16

Observamos claramente por medio del histograma la falta de normalidad y luego lo podemos comprobar de manera gráfica y analítica. En R, podemos hacer uso de la función boxcox() de la librería MASS para estimar el parámetro de transformación \(\lambda\).

library(MASS)
# Debemos construir modelos lineales
modelo1<-lm(x1~1)
modelo2<-lm(x2~1)

# Estimación de lambda gráfico
par(mfrow=c(1,2))
b1<-boxcox(modelo1)
b2<-boxcox(modelo2)

# Estimación  de lambda 
lambda1 <- b1$x[which.max(b1$y)]
round(lambda1,1)
## [1] 0.3
lambda2 <- b2$x[which.max(b2$y)]
round(lambda2,1)
## [1] -0.2

Por lo tanto, debemos hacer las siguientes transformaciones para buscar conseguir una aproximación a la normalidad.

x1_trans <- (x1^lambda1-1)/lambda1  
x2_trans <- (x2^lambda2-1)/lambda2

Analicemos ahora si logran cumplir la normalidad estas variables transformadas.

# Histogramas
par(mfrow=c(1,2))
hist(x1_trans,probability=T,main=" ",xlab="Datos exp_trans",ylab="Densidad")
x <- seq(min(x1_trans)-1, max(x1_trans)+1, length = 100)
y <- dnorm(x, mean(x1_trans),sd(x1_trans))
lines(x, y, col = "red",lwd=3)
hist(x2_trans,probability=T,main=" ",xlab="Datos F_trans",ylab="Densidad")
x <- seq(min(x2_trans)-1, max(x2_trans)+1, length = 100)
y <- dnorm(x, mean(x2_trans),sd(x2_trans))
lines(x, y, col = "red",lwd=3)

# Q-Q plots
par(mfrow=c(1,2))
qqnorm(x1_trans, main = "", xlab = "Cuantiles teóricos", ylab = "Cuantiles muestrales")
qqline(x1_trans,col="blue",lwd=2)
qqnorm(x2_trans, main = "", xlab = "Cuantiles teóricos", ylab = "Cuantiles muestrales")
qqline(x2_trans,col="blue",lwd=2)

# Pruebas
CVM1_trans <- cvm.test(x1_trans)
CVM1_trans 
## 
##  Cramer-von Mises normality test
## 
## data:  x1_trans
## W = 0.01804, p-value = 0.982
CVM2_trans <- cvm.test(x2_trans)
CVM2_trans 
## 
##  Cramer-von Mises normality test
## 
## data:  x2_trans
## W = 0.031888, p-value = 0.8183

Ejemplo 3:

Si tienes datos que contienen valores negativos o ceros, entonces la transformación de Box-Cox no es apropiada para esas variables. Para manejar datos que incluyen ceros y valores negativos, podrías considerar la transformación Yeo-Johnson, que es una extensión de la transformación de Box-Cox y puede manejar una gama más amplia de valores.

n <- 100
set.seed(2023)
x1 <- rgamma(n,5,1)-3  # gamma con valores negativos
x2 <- rnorm(n,-10,1)


# Histogramas
par(mfrow=c(1,2))
hist(x1,probability=T,main=" ",xlab="Datos gamma",ylab="Densidad")
hist(x2,probability=T,main=" ",xlab="Datos normales",ylab="Densidad")

# Pruebas gráficas
par(mfrow=c(1,2))
qqnorm(x1,main=" ",xlab="Cuantiles teóricos",ylab="Cuantiles muestrales")
qqline(x1,col="blue",lwd=2)
qqnorm(x2,main=" ",xlab="Cuantiles teóricos",ylab=" ")
qqline(x2,col="blue",lwd=2)

# Pruebas analíticas
SW1 <- shapiro.test(x1)
SW1
## 
##  Shapiro-Wilk normality test
## 
## data:  x1
## W = 0.9226, p-value = 1.94e-05
SW2 <- shapiro.test(x2)
SW2
## 
##  Shapiro-Wilk normality test
## 
## data:  x2
## W = 0.9892, p-value = 0.6006

Observamos claramente por medio del histograma la falta de normalidad y luego lo podemos comprobar de manera gráfica y analítica. La librería car incluye la función powerTransform() que permite aplicar la transformación Yeo-Johnson a datos multivariados.

library(car)

# Aplicar la transformación Yeo-Johnson
trans1 <- powerTransform(x1, family = "yjPower")
trans2 <- powerTransform(x2, family = "yjPower")

# Valor estimado de lambda
lambda1 <- trans1$roundlam
round(lambda1,1)
##  x1 
## 0.5
lambda2 <- trans2$roundlam
round(lambda2,1)
## x2 
##  1

Por lo tanto, debemos hacer las siguientes transformaciones para buscar conseguir una aproximación a la normalidad. Es obvio que para la variable 2 no es necesaria ninguna transformación.

x1_trans <- ifelse(x1>=0,((x1+1)^(lambda1)-1)/lambda1, 
                  -(((-x1+1)^(2-lambda1)-1)/(2-lambda1)))

Analicemos ahora si logra cumplir la normalidad esta variable transformada.

# Histogramas
par(mfrow=c(1,2))
hist(x1_trans,probability=T,main=" ",xlab="Datos gamma_trans",ylab="Densidad")
x <- seq(min(x1_trans)-1, max(x1_trans)+1, length = 100)
y <- dnorm(x, mean(x1_trans),sd(x1_trans))
lines(x, y, col = "red",lwd=3)
hist(x2,probability=T,main=" ",xlab="Datos normales",ylab="Densidad")
x <- seq(min(x2)-1, max(x2)+1, length = 100)
y <- dnorm(x, mean(x2),sd(x2))
lines(x, y, col = "red",lwd=3)

# Pruebas gráficas
par(mfrow=c(1,2))
qqnorm(x1_trans,main=" ",xlab="Cuantiles teóricos",ylab="Cuantiles muestrales")
qqline(x1_trans,col="blue",lwd=2)
qqnorm(x2,main=" ",xlab="Cuantiles teóricos",ylab=" ")
qqline(x2,col="blue",lwd=2)

# Pruebas analíticas
SW1 <- shapiro.test(x1_trans)
SW1
## 
##  Shapiro-Wilk normality test
## 
## data:  x1_trans
## W = 0.98879, p-value = 0.5678
SW2 <- shapiro.test(x2)
SW2
## 
##  Shapiro-Wilk normality test
## 
## data:  x2
## W = 0.9892, p-value = 0.6006

Nota importante: La función powerTransform() también incluye la transformación Box-Cox. Para usarla debe establecer el argumento family="bcPower".

Ejemplo 4:

Finalmente, trabajaremos con la función bestNormalize() de la librería que lleva su mismo nombre, bestNormalize. Vale la pena resaltar que la última versión de esta librería fue publicada en 18 de agosto de 2023.

Esta diseñada para ayudar a encontrar la mejor transformación para normalizar un vector. Hay muchas técnicas que se han desarrollado con este objetivo, sin embargo, cada una ha estado sujeta a sus propias fortalezas y debilidades, y no está claro cómo decidir cuál funcionará mejor hasta que se observen los datos. Este paquete analiza una variedad de transformaciones posibles y devolve la mejor, es decir, la que hace consigue una mejor aproximación a la normal.

Para evaluar la eficacia de la técnica de normalización, la función bestNormalize() implementa validaciones cruzadas repetidas para estimar el estadístico de Pearson dividido por sus grados de libertad. Esto se llama “estadística de normalidad” y si es cercano a 1 (o menos), entonces se puede considerar que la transformación funciona bien. La función está diseñada para seleccionar la transformación que produce el valor \(P/df\) más bajo.

Además, esta última versión presenta una nueva adaptación de una técnica de normalización, que definen normalización de cuantiles ordenados orderNorm() o ORQ(). Transforma los datos basándose en un mapeo de rangos a la distribución normal, permitiendo garantizar datos transformados distribuidos normalmente. La adaptación utiliza una aproximación logit desplazada para realizar la transformación en datos recién observados fuera del dominio original. En datos nuevos dentro del dominio original, la transformación utiliza la interpolación lineal de la transformación ajustada.

Es importante tener en cuenta que estamos trabajando con la palabra normalizar, la cual estamos empleando para transformar un vector de datos de tal manera que los valores transformados sigan una distribución gaussiana (o equivalentemente, una curva de campana). Pero en literatura puede confundirse con “estandarizar, que se refiere a cuando tenemos un vector de datos normales, restamos por su media y dividimos por su desviación estándar, para conseguir datos normales estándar.

library(bestNormalize)
n <- 250
set.seed(100)
x <- rgamma(n,1,1)

par(mfrow=c(1,1))
hist(x,probability=T,main=" ",col="cyan",xlab="Datos gamma",ylab="Densidad")

Estos datos claramente no son normales. Usemos la función bestNormalize() para realizar un conjunto de transformaciones potenciales y veamos cómo funciona cada método.

# 1. Trnasformacion arcsinh 
arcsinh_trans <- arcsinh_x(x)
arcsinh_trans
## Standardized asinh(x) Transformation with 250 nonmissing obs.:
##  Relevant statistics:
##  - mean (before standardization) = 0.7383146 
##  - sd (before standardization) = 0.5458515
# 2. Transformacion Box-Cox 
boxcox_trans <- boxcox(x)
boxcox_trans
## Standardized Box Cox Transformation with 250 nonmissing obs.:
##  Estimated statistics:
##  - lambda = 0.3254863 
##  - mean (before standardization) = -0.3659267 
##  - sd (before standardization) = 0.9807881
# 3. Transformacion Yeo-Johnson
yeojohnson_trans <- yeojohnson(x)
yeojohnson_trans
## Standardized Yeo-Johnson Transformation with 250 nonmissing obs.:
##  Estimated statistics:
##  - lambda = -0.7080476 
##  - mean (before standardization) = 0.4405464 
##  - sd (before standardization) = 0.2592004
# 4. Trnasformacion OQR
orderNorm_trans <- orderNorm(x)
orderNorm_trans
## orderNorm Transformation with 250 nonmissing obs and no ties 
##  - Original quantiles:
##    0%   25%   50%   75%  100% 
## 0.001 0.268 0.721 1.299 4.161
# 5. Mejor transformacion automatica
best_trans <- bestNormalize(x)
best_trans
## Best Normalizing transformation with 250 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 1.7917
##  - Box-Cox: 1.0442
##  - Center+scale: 3.0102
##  - Double Reversed Log_b(x+a): 4.7638
##  - Exp(x): 9.5306
##  - Log_b(x+a): 1.7072
##  - orderNorm (ORQ): 1.1773
##  - sqrt(x + a): 1.144
##  - Yeo-Johnson: 1.1875
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##  
## Based off these, bestNormalize chose:
## Standardized Box Cox Transformation with 250 nonmissing obs.:
##  Estimated statistics:
##  - lambda = 0.3254863 
##  - mean (before standardization) = -0.3659267 
##  - sd (before standardization) = 0.9807881

Para examinar cómo se desempeña cada uno de ellos, podemos visualizar los valores transformados en un histograma.

par(mfrow = c(2,2))
hist(arcsinh_trans$x.t,main="Transformacion arcsinh",xlab="",ylab=" ")
hist(boxcox_trans$x.t,main="Transformacion Box-Cox",xlab="",ylab=" ")
hist(yeojohnson_trans$x.t,main="Transformacion Yeo-Johnson",xlab="",ylab=" ")
hist(orderNorm_trans$x.t,main="Transformacion OQR",xlab="",ylab=" ")

La mejor transformación se presenta a continuación, que en este caso fue Box-Cox.

# Histograma
hist(best_trans$x.t,probability=T,main="Mejor transformacion",col="cyan",xlab="",ylab=" ")
x <- seq(min(best_trans$x.t)-1, max(best_trans$x.t)+1, length = 100)
y <- dnorm(x, mean(best_trans$x.t),sd(best_trans$x.t))
lines(x, y, col = "blue",lwd=3)

# Prueba gráfica
qqnorm(best_trans$x.t,main="Mejor transformacion",xlab="Cuantiles teóricos",ylab="Cuantiles muestrales")
qqline(best_trans$x.t,col="blue",lwd=2)

# Prueba analítica
SW <- shapiro.test(best_trans$x.t)
SW
## 
##  Shapiro-Wilk normality test
## 
## data:  best_trans$x.t
## W = 0.99289, p-value = 0.2767

Podemos comprobar el buen desempeño de la transformación.

Algo muy interesante dentro de la función bestNormalize() es que guarda las estadísticas de normalidad estimadas que obtuvo para cada pliegue y repetición de la validación cruzada. Estos se pueden visualizarlos con relativa facilidad a través de boxplots, que pueden dar una idea de hasta qué punto la estadística de normalidad prefiere realmente la transformación, o si no importa tanto qué transformación se elija.

boxplot(log10(best_trans$oos_preds),las=2,cex.axis=0.7)

En este ejemplo, Box-Cox, ORQ, raíz cuadrada, parecen tener un desempeño similar entre sí, mientras que Yeo-Johnson tiene un rendimiento inferior.

Ejercicio de clase:

Para aplicar lo aprendido, trabajaremos con dos conjuntos de datos reales. Redacta un pequeño informe de la solución de los ejercicios. No se puede entregar solo el código.

  1. Utiliza el conjunto de datos obtenido de Gapminder que incluye información de 142 países sobre expectativa de vida, ingreso per cápita, y población de 1952 a 2007. Para ello debemos bajar el conjunto de datos que se encuentra disponible en el repositorio de Github de la Dra. Jenny Bryan con el siguiente código:
devtools::install_github("jennybc/gapminder")
# Cargamos los datos
gapminder <- gapminder::gapminder

Trabaja en la transformación de las variables población (pop), expectativa de vida (lifeExp) y el ingreso per cápita (gdpPercap) para un año cualquier en específico. Por ejemplo, con el siguiente código obtienes la población del año que quieras.

# Datos de población 
data_pob <- gapminder$pop[gapminder$year == XXX]

Apique todo lo aprendido en clase y transforme estos datos para conseguir normalidad.

  1. El conjunto autotrader de datos se extrajo del sitio web de autotrader como parte de la librería bestNormalize. Esta base de datos contiene 10 variables medidas en 6283 observaciones.
data("autotrader")

Trabaja en la transformación de las variables kilometraje (mileage), antigüedad (year) y precio (price).

Nota importante: Por cuestión de ejemplificar la clase, solo trabaje con la prueba de Shapiro-Wilk, pero es muy importante tener en cuenta que existen diferentes pruebas. A continuación sus funciones:

# Pruebas
library(nortest)
SW <- shapiro.test() # Shapiro-Wilk
KS <- ks.test()      # Kolmogorov-Smirnov
L  <- lillie.test()  # Lilli
SF <- sf.test()      # Shapiro-Francia
AD <- ad.test()      # Anderson-Darling
CVM <- cvm.test()    # Cramer-Von Mises

\[ \]