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.
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)
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}\).