Normalidad


Ejemplo de verificación de la normalidad

Comenzamos por cargar y ver los datos

library(latex2exp)

Datos = read.csv("ejemplo4.csv")
head(Datos)
##   x          y
## 1 1  0.9395244
## 2 2  2.7698225
## 3 3  7.0587083
## 4 4  9.0705084
## 5 5 13.6292877
## 6 6 20.7150650

donde el conjunto de datos ya había sido trabajado con anterioridad.

Linealidad y homocedasticidad

Vemos el scatterplot inicial

# Realizamos el ajuste del modelo
fit=lm(y~x, data=Datos)

# Scatterplot
par(mfrow=c(1,1)) 
par(mar=c(4, 5, 1, 1))
plot(Datos$x, Datos$y, xlab = TeX("$x$"), ylab=TeX("$y$"))
abline(fit, col="red")

Procedemos a transformar los datos para arreglar los problemas de linealidad y homocedasticidad

#Se transforma X a X^2 para lograr linealidad
Datos$Xprima=Datos$x^2

# Ajuste del modelo con los datos transformados
fit2=lm(y~Xprima, data=Datos)


#gráficas para linealidad y homocedasticidad
par(mar=c(4, 5, 1, 1))
par(mfrow=c(1,2))
plot(fit2, 1)
plot(fit2, 3)

Normalidad

Como en el ejemplo pasado de normalidad, analizaremos los casos cuando \(n\) es grande y cuando es pequeño.

  • Consideremos cuando \(n\) es grande. Comenzamos por ver la Q-Qplot que R trae por defecto

    #gráfica QQplot directa de R
    
    par(mar=c(4, 5, 1, 1))
    par(mfrow=c(1,1))
    plot(fit2, 2)

    notemos que en este caso no hay evidencia en contra de la normalidad, al menos por la información que nos arroja la Q-Qplot. Procedemos pues al análisis de las pruebas de hipótesis

    # Se usan los residuales estandarizados y los cuantiles de una normal
    
    # R tiene una función para crear errores de forma automatizada
    library(broom)
    
    # Agregamos los errores calculados al dataframe original
    Datosfit2=augment(fit2)

    y efectuamos

    shapiro.test(Datosfit2$.std.resid)
    ## 
    ##  Shapiro-Wilk normality test
    ## 
    ## data:  Datosfit2$.std.resid
    ## W = 0.99368, p-value = 0.9256

    donde no se está rechazando \(H_{0}\) (que es justamente lo que buscamos). Después, con las restantes dos pruebas

    library(nortest)
    library(normtest)
    
    nortest::lillie.test(Datosfit2$.std.resid)
    ## 
    ##  Lilliefors (Kolmogorov-Smirnov) normality test
    ## 
    ## data:  Datosfit2$.std.resid
    ## D = 0.05564, p-value = 0.6265
    normtest::jb.norm.test(Datosfit2$.std.resid)
    ## 
    ##  Jarque-Bera test for normality
    ## 
    ## data:  Datosfit2$.std.resid
    ## JB = 0.19216, p-value = 0.905

    en las que también se rechaza \(H_{0}\), lo cual nos indicaría con mayor fuerza que no hay evidencia en contra de la normalidad.

    Podemos ver lo que pasa con el modelo en el que no transformamos a la \(x\):

    # Q-Qplot
    par(mar=c(4, 5, 1, 1))
    par(mfrow=c(1,1))
    plot(fit, 2)

    y las pruebas de hipótesis

    Datosfit=augment(fit)
    shapiro.test(Datosfit$.std.resid)
    ## 
    ##  Shapiro-Wilk normality test
    ## 
    ## data:  Datosfit$.std.resid
    ## W = 0.92316, p-value = 2.08e-05
    nortest::lillie.test(Datosfit$.std.resid)
    ## 
    ##  Lilliefors (Kolmogorov-Smirnov) normality test
    ## 
    ## data:  Datosfit$.std.resid
    ## D = 0.1273, p-value = 0.0004008
    normtest::jb.norm.test(Datosfit$.std.resid)
    ## 
    ##  Jarque-Bera test for normality
    ## 
    ## data:  Datosfit$.std.resid
    ## JB = 8.3681, p-value = 0.0205

    en las que se rechaza \(H_{0}\), es decir, se rechaza que \(e_{st_{1}},\ldots, e_{st_{n}}\) provienen de una normal.

    Con base en lo anterior tenemos evidencia contundente en contra de la normalidad en el modelo no transformado.
  • Procedemos a ver el caso en se considera a la \(n\) pequeña. Comenzamos por obtener los errores studentizados del ajuste del modelo con los datos transformados

    library(MASS)
    ResStudent=studres(fit2)

    y después vemos la Q-Q plot

    # Para esta qqplot requerimos de la librería car
    library(car)
    
    qqPlot(ResStudent, dist="t", df=length(ResStudent)-3,  envelope=FALSE)

    ## [1] 72 44

    o laternativamente podemos utilizar

    qqnorm(ResStudent, pch = 1, frame = FALSE)
    qqline(ResStudent, col = "steelblue", lwd = 2)

    siendo las funciones qqnorm y qqline incluidas como base en R. Podemos coparar dichas gráficas

    par(mar=c(4, 5, 1, 1))
    par(mfrow=c(1,2))
    qqnorm(ResStudent, pch = 1, frame = FALSE)
    qqline(ResStudent, col = "steelblue", lwd = 2)
    qqPlot(ResStudent, dist="t", df=length(ResStudent)-3,  envelope=FALSE)

    ## [1] 72 44

    Finalmente, la prueba de hipótesis de Kolmogorov_Smirnov:

    ks.test(ResStudent, "pt", length(ResStudent)-3)
    ## 
    ##  One-sample Kolmogorov-Smirnov test
    ## 
    ## data:  ResStudent
    ## D = 0.055383, p-value = 0.9189
    ## alternative hypothesis: two-sided

    en el que no se rechaza \(H_{0}\). Así, nuevamente no estamos encontrando evidencia en contra de la normalidad en el modelo con los datos transformados. No obstante, ene l modelo con los datos no transformados tenemos que

    # Errores studentizados
    ResStudent2=studres(fit)
    
    # Gráfica que trae por defecto R
    qqPlot(ResStudent2, dist="t", df=length(ResStudent2)-3,  envelope=FALSE)

    ## [1] 70 11

    y la prueba de hipótesis

    ks.test(ResStudent2, "pt", length(ResStudent)-3)
    ## 
    ##  One-sample Kolmogorov-Smirnov test
    ## 
    ## data:  ResStudent2
    ## D = 0.12593, p-value = 0.08386
    ## alternative hypothesis: two-sided
    en la que se rechaza \(H_{0}\).