Ejemplos de tranformaciones a la normalidad

Primero veamos la transformacion de Box-Cox vista en clase

modelo<-lm(Volume~., trees)
summary(modelo)

## 
## Call:
## lm(formula = Volume ~ ., data = trees)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.4065 -2.6493 -0.2876  2.2003  8.4847 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -57.9877     8.6382  -6.713 2.75e-07 ***
## Girth         4.7082     0.2643  17.816  < 2e-16 ***
## Height        0.3393     0.1302   2.607   0.0145 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.882 on 28 degrees of freedom
## Multiple R-squared:  0.948,  Adjusted R-squared:  0.9442 
## F-statistic:   255 on 2 and 28 DF,  p-value: < 2.2e-16

hist(modelo$res)

par(mfrow=c(2,2))
plot(modelo)

Graficamos la log-likelihood en función de lambda

library(MASS)
par(mfrow=c(1,1))
b<-boxcox(modelo)

Como el máximo no se localiza a simple vista, debemos localizarlo

b$x[which(b$y==max(b$y))]

## [1] 0.3030303

Una vez que conocemoes el coeficiente que maximizó la log-lik, aplicamos una transformación al modelo

modelo2<-lm((Volume^(0.3)-1)/0.3~.,trees)
par(mfrow=c(2,2))
plot(modelo2)

bestNormalize

Como se mencionó en clase, tambien es poisble usar el paquete bestNormaliza que contiene un conjunto de funciones para encontrar transformaciones que se usan para ‘normalizar’ los datos. Recuerde que en nuestro caso, queremos que se cumple el supuesto sobre los errores.

Hay varias opciones de transformación a la normalidad, cada una con sus propias ventajas y limitaciones. Además el paquete permite analizar transformaciones individualemte o en grupo. Algunas transformaciones incluidas son:

• Box-Cox (datos positivos)
• ArcSinh
• Yeo-Johnson
• Cuantiles -log10(a+x)

En el ejemplo anterior, se podría aplicar alguna de las anteriores:

library(bestNormalize)

## 
## Attaching package: 'bestNormalize'

## The following object is masked from 'package:MASS':
## 
##     boxcox

arcsinh_x(modelo$residuals)

## Standardized asinh(x) Transformation with 31 nonmissing obs.:
##  Relevant statistics:
##  - mean (before standardization) = -0.07244041 
##  - sd (before standardization) = 1.767576

orderNorm(modelo$res)

## orderNorm Transformation with 31 nonmissing obs and no ties 
##  - Original quantiles:
##     0%    25%    50%    75%   100% 
## -6.406 -2.649 -0.288  2.200  8.485

o podemos usar la funcion que nos dice cuales son las mejores transformaciones para nuestros datos:

bestNormalize(modelo$residuals)

## Best Normalizing transformation with 31 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 1.7133
##  - Center+scale: 1.42
##  - Exp(x): 2.9133
##  - Log_b(x+a): Inf
##  - orderNorm (ORQ): 1.62
##  - sqrt(x + a): Inf
##  - Yeo-Johnson: 1.2867
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##  
## Based off these, bestNormalize chose:
## Standardized Yeo-Johnson Transformation with 31 nonmissing obs.:
##  Estimated statistics:
##  - lambda = 0.896797 
##  - mean (before standardization) = -0.3195476 
##  - sd (before standardization) = 3.724466

por ejemplo, en este caso podemos observar que es la transformación de Yeo-Johnson, muy similar a BoxCox

modelo3<- lm((Volume^(0.89)-1)/0.89 ~., data=trees)
par(mfrow=c(2,2))
plot(modelo2)

Y podemos comparar los errores de modelos a través de graficas o pruebas de hipótesis

par(mfrow=c(1,2))
hist(modelo2$residuals)
hist(modelo3$residuals)