Verificación de supuestos del modelo de regresión lineal múltiple

El diagnóstico desempeña un papel importante en el desarrollo y la evaluación de los modelos de regresión múltiple. La mayoría de los procedimientos de diagnóstico para la regresión lineal simple que se vieron con anterioridad se aplican directamente a la regresión múltiple. En este apartado se repasarán el procedimiento de diagnóstico. También se han desarrollado muchos diagnósticos especializados para la regresión múltiple

Ajustados \(\hat{y_i}\) y residuales \(e_i\)

Se utilizarán los valores ajustados \(\hat{y_i}\) (Ecuación 2) por el modelo de datos de la Ecuación 1 y los residuales \(e_i\) (Ecuación 3).

\[\begin{align} y_i = \beta_0 + \beta_1x_{i1} + \beta_2x_{i2} + ...+ \beta_{p-1}x_{i,~p-1} + \epsilon_i \end{align} \tag{1}\]

\[\begin{align} \hat{Y} = X\hat{\beta} \end{align} \tag{2}\]

\[\begin{align} e &= Y - \hat{Y}\\ e &= Y - X\hat{\beta} \end{align} \tag{3}\]

Verificación de la normalidad.

Normalidad gráfica

Para realizar la verificación gráfica de la normalidad es necesario construir un gráfico de dispersión cuantiles teóricos vs cuantíles de los residuales, o gráfico cuantil cuantil (QQplot). Es necesario poseer los cálculos de los residuales \(e_i\).

Se siguen los siguientes pasos:

Ordenar los \(n\) residuales \(e_i\) en orden ascendente y asignarles en ese orden los valores de \(1\) a \(n\). Denotamos \(r_i\) con \(i=1,2,3,..,n\) los datos de los residuales \(e_i\) en orden creciente.
Calcular una posición de graficación \(p_i\) para cada dato \(r_i\) en función de su rango y del total de observaciones como:

\[\begin{align} p_i=\frac{i-0.5}{n}~~;~~i=1,2,...,n \end{align}\]

Calcular la normal inversa \(Z_i\) de \(p_i\):

\[\begin{align} Z_i = \phi^{-1} \left( \frac{1-0.5}{n} \right) \end{align}\]

Hacer gráfico de dispersión: en el eje de las abscisas por \(Z_i\) y en el eje de las ordenadaslos \(r_i\)

Si los residuos siguen una distribución normal, al graficarlos tienden a quedar alineados en una línea recta; por lo tanto, si claramente no se alinean se concluye que el supuesto de normalidad está siendo violado.

Cabe enfatizar el hecho de que el ajuste de los puntos a una recta no tiene que ser perfecto. No deberían existir desviaciones importantes a la línea recta para sospechar sobre el incumplimiento del supuesto de normalidad. En todo caso, siempre se recomienda realizar el proceso de verificaicón formal o analítico

Ejemplo verificación normalidad gráfica.

Dwaine Studios, Inc. gestiona estudios de fotografía en \(21\) ciudades intermedias. Estos estudios están especializados en retratos de niños. La empresa está considerando expandirse a otras ciudades intermedias y desea investigar si las ventas en una ciudad pueden predecirse a partir del número de personas de \(16\) años o menos en la comunidad y el ingreso per cápita en la comunidad. En la Tabla 1 se muestran los datos de estas variables correspondientes al año más reciente de las \(21\) ciudades en las que opera actualmente Dwaine Studios.

Tabla 1: Datos ejemplo ventas de Dwaine Studios.

Observación \(i\)	\(x_1\): Población menor 16 años (miles de personas)	\(x_2\): Ingreso per cápita (miles de dólares)	\(y\): ventas (miles de dólares)
1	68.5	16.7	174.4
2	45.2	16.8	164.4
3	91.3	18.2	244.2
4	47.8	16.3	154.6
5	46.9	17.3	181.6
6	66.1	18.2	207.5
7	49.5	15.9	152.8
8	52.0	17.2	163.2
9	48.9	16.6	145.4
10	38.4	16.0	137.2
11	87.9	18.3	241.9
12	72.8	17.1	191.1
13	88.4	17.4	232.0
14	42.9	15.8	145.3
15	52.5	17.8	161.1
16	85.7	18.4	209.7
17	41.3	16.5	146.4
18	51.7	16.3	144.0
19	89.6	18.1	232.6
20	82.7	19.1	224.1
21	52.3	16.0	166.5

Ajuste del modelo

#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Modelo
modelo <- lm(ventas~poblacion+ingreso)

Cálculo de residuales

#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Modelo
modelo <- lm(ventas~poblacion+ingreso)

#Residuales
residuales <- modelo$residuals

Gráfico cuantil- cuantil en ggplot

#Librerías
library(ggplot2)
library(plotly)
#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Modelo
modelo <- lm(ventas~poblacion+ingreso)

#Residuales
residuales <- modelo$residuals

# Residuales organizados
ri <- sort(residuales)
# Posición de graficación
n <- length(ventas)
i <- c(1:length(ventas))
pi <- (i-0.5)/n

# Normal inversa de pi
zi <- qnorm(pi, mean=0, sd=1)

#Grafico
datos <- data.frame(ri,zi)
grafico <- ggplot(data=datos, aes(x=zi, y=ri)) + geom_point()
grafico <- grafico + xlab("Cuantiles teóricos") + ylab("Cuantiles datos")
grafico <- grafico + geom_smooth(method = "lm", se = FALSE)
grafico <- grafico + theme_classic()
ggplotly(grafico)

Figura 1: Gráfico Cuantil - Cuantil (QQPLOT)

En la Figura 1 se observa un patrón moderadamente lineal. Si se calcula el coeficiente de correlación, con la sentencia cor de R, entre los cuantiles teóricos y los cuantiles de los residuales se obtiene un coeficiente de correlación de \(0.9790896\), este valor alto de coeficiente de correlación puede ayudar a verificar la normalidad de forma gráfica, como se observa a continuación:

#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Modelo
modelo <- lm(ventas~poblacion+ingreso)

#Residuales
residuales <- modelo$residuals

# Residuales organizados
ri <- sort(residuales)
# Posición de graficación
n <- length(ventas)
i <- c(1:length(ventas))
pi <- (i-0.5)/n

# Normal inversa de pi
zi <- qnorm(pi, mean=0, sd=1)

#Correlación
correlacion<-cor(zi,ri)
correlacion

[1] 0.9790896

Se puede realizar la verificación son sentencias de R:

#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Modelo
modelo <- lm(ventas~poblacion+ingreso)

#Residuales
residuales <- modelo$residuals

# Gráfico cuantil-cuantil
qqnorm(residuales, xlab = "Cuantiles teóricos", ylab="Cuántiles residuales")
qqline(residuales, col = "blue", lwd = 2)

Verificación formal - analítica supuesto de normalidad mediante Test Kolmogorov-Smirnov

Consideremos una muestra aleatoria de datos \(x_1,x_2,...x_n\) que procede de cierta función desconocida denotada \(F(x)\). Se quiere verificar si dichos datos fueron generados por un proceso normal, mediante las hipótesis estadísticas.

Planteamiento de hipótesis

\[\begin{align} H_o&: e_{i}\sim N(\mu,\sigma^2)\\H_1&: e_{i} \nsim N(\mu,\sigma^2) \\ \forall ~i&=1,2,...,n \end{align}\]

Organización de los residuales \(e_{i}\) de manera ascendente, denotando los datos ordenados como:

\[\begin{align} r_i=x_{(1)},x_{(2)},...,x_{(n)},~con~i=1,2,3...n \end{align}\]

Calculos probabilidad teórica (percentil) \(P_i\) de la forma:

\[\begin{align} P_i=\frac{i}{n},~ \forall ~i&=1,2,...,n \end{align}\]

Estandarización los residuales organizados \(r_i\) de la siguiente manera:

\[\begin{align} Z_i=\frac{x_i-\bar{x}}{S(x_i)},~ \forall ~i&=1,2,...,n \end{align}\]

\(\bar{x}\) corresponde a la media muestra para los \(r_i\).
\(S(x_i)\) es la desviación estándar de la muestra de \(r_i\).

Calculo de la probabilidad \(P(Z_i)\) para una distribución normal estándar.

\[\begin{align} P(Z_i)=\phi \left[\frac{x_i-\bar{x}}{S(x_i)} \right], ~ \forall ~i&=1,2,...,n \end{align}\] 5. Calculo de distancias \(D_1\) y \(D_2\):

Distancia 1 \(D_1\).

\[\begin{align} D_1=|P(Z_i)-P_i|, ~ \forall ~i&=1,2,...,n \end{align}\]

Distancia 2 \(D_2\).

\[\begin{align} D_2=|P(Z_i)-P_{i-1}|,~ \forall ~i&=1,2,...,n \end{align}\]

Escojo el valor para el estadístico de prueba \(D\)

\[\begin{align} D=max\left\lbrace D_1,D_2 \right\rbrace \end{align}\]

Comparamos el estadístico de prueba \(D\) con el estadístico teórico \(KS\)

\[\begin{align} KS=\frac{C_\alpha}{K(n)} \end{align}\]

\(C_\alpha\) viene dado por tablas de la siguiente manera (Tabla 2).

Tabla 2: \(C_{\alpha}\)

\(C_{\alpha}\) por distribución	\(\alpha=0.1\)	\(\alpha=0.05\)	\(\alpha=0.01\)
Normal	0.819	0.895	1.035
Exponencial	0.990	1.094	1.308
Weibull	0.760	0.819	0.944

El valor de \(K(n)\) también se encuentra tabulado, y lo encontramos de la siguiente forma (Tabla 3):

Tabla 3: \(K(n)\).

\(K(n)\) por distribución	\(K(n)\)
Normal	\(\sqrt{n}-0.01+\frac{0.85}{\sqrt{n}}\)
Exponencial	\(\sqrt{n}+0.12+\frac{0.11}{\sqrt{n}}\)
Weibull	\(\sqrt{n}\)

\[Si~D>KS \rightarrow Rechazo~H_o\]

Ejemplo verificación normalidad Test Kolmogorov-Smirnov.

Ajuste del modelo y cálculo de residuales

#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Modelo
modelo <- lm(ventas~poblacion+ingreso)

#Residuales
residuales <- modelo$residuals

Cálculo de estadístico de prueba \(D\), estadístico teórico \(KS\) y comparación

#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Modelo
modelo <- lm(ventas~poblacion+ingreso)

#Residuales
residuales <- modelo$residuals

# Residuales organizados
ri <- sort(residuales)

# Probabilidad teórica
n <- length(ventas)
i <- c(1:length(ventas))
pi <- i/n

#Estandarización de residuales 
zi <- (ri-mean(ri))/sd(ri)

#Probabilidad normal estándar
p_zi <- pnorm(zi, mean=0, sd=1)

#Cálculo de distancias D1 Y D2
D1 <- abs(p_zi - pi)
D2_1 <- p_zi[1]-0 
D2_2<- abs(p_zi[2:length(p_zi)] - pi[1:(length(pi)-1)])
D2 <- c(D2_1, D2_2)

#Estadístico de prueba KS
D<- max(D1,D2)

# Estadistico teórico para alpha=0.05
KS <- 0.895/(sqrt(length(ventas))-0.01+(0.85/sqrt(length(ventas))))

# Comparación
if(D>KS){
print("Existe evidencia suficiente para rechazar H0, por lo tanto, los residuales no provienen de una distribución normal")
} else{
  
print("No existe evidencia suficiente para rechazar H0, por lo tanto, los residuales no provienen de una distribución normal")
  }

[1] "No existe evidencia suficiente para rechazar H0, por lo tanto, los residuales no provienen de una distribución normal"

Los cálculos anteriores se muestran de manera resumida en la Tabla 4.

Tabla 4: Verificación supuesto de normalidad Kolmogorov-Smirnov

Observación \(i\)	\(y_i\): ventas	\(\hat{y_i}\): ajustados	\(e_i\): residuales	\(r_i\): residuales organizados	\(P_i\): probabilidad teórica	\(Z_i\): residuales estandarizados	\(P(Z_i)\): probabilidad normal estándar	D1	D2
1	174.4	187.1841	-12.7841146	-18.4238900	0.0476190	-1.7643134	0.0388396	0.0087794	0.0388396
2	164.4	154.2294	10.1705737	-15.0013134	0.0952381	-1.4365597	0.0754216	0.0198165	0.0278025
3	244.2	234.3963	9.8036764	-13.1132116	0.1428571	-1.2557508	0.1046032	0.0382540	0.0093651
4	154.6	153.3285	1.2714690	-12.7841146	0.1904762	-1.2242357	0.1104317	0.0800445	0.0324255
5	181.6	161.3849	20.2150722	-12.3381967	0.2380952	-1.1815336	0.1186954	0.1193998	0.0717808
6	207.5	197.7414	9.7585779	-6.2160615	0.2857143	-0.5952641	0.2758335	0.0098808	0.0377382
7	152.8	152.0551	0.7449178	-6.0849209	0.3333333	-0.5827058	0.2800457	0.0532876	0.0056686
8	163.2	167.8666	-4.6666316	-4.6666316	0.3809524	-0.4468872	0.3274783	0.0534741	0.0058551
9	145.4	157.7382	-12.3381967	0.3539791	0.4285714	0.0338978	0.5135207	0.0849493	0.1325683
10	137.2	136.8460	0.3539791	0.6530062	0.4761905	0.0625333	0.5249309	0.0487405	0.0963595
11	241.9	230.3874	11.5126289	0.7449178	0.5238095	0.0713350	0.5284344	0.0046249	0.0522440
12	191.1	197.1849	-6.0849209	1.2714690	0.5714286	0.1217587	0.5484550	0.0229736	0.0246454
13	232.0	222.6857	9.3142995	1.6129777	0.6190476	0.1544624	0.5613774	0.0576702	0.0100512
14	145.3	141.5184	3.7815611	3.7815611	0.6666667	0.3621308	0.6413729	0.0252938	0.0223253
15	161.1	174.2132	-13.1132116	9.3142995	0.7142857	0.8919584	0.8137924	0.0995067	0.1471257
16	209.7	228.1239	-18.4238900	9.4356009	0.7619048	0.9035745	0.8168895	0.0549847	0.1026037
17	146.4	145.7470	0.6530062	9.7585779	0.8095238	0.9345035	0.8249779	0.0154541	0.0630731
18	144.0	159.0013	-15.0013134	9.8036764	0.8571429	0.9388222	0.8260890	0.0310539	0.0165652
19	232.6	230.9870	1.6129777	10.1705737	0.9047619	0.9739571	0.8349611	0.0698008	0.0221818
20	224.1	230.3161	-6.2160615	11.5126289	0.9523810	1.1024754	0.8648725	0.0875085	0.0398894
21	166.5	157.0644	9.4356009	20.2150722	1.0000000	1.9358410	0.9735564	0.0264436	0.0211755

Se puede realizar el test de Kolmogorov - Smirnov con sentencias de R:

#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Modelo
modelo <- lm(ventas~poblacion+ingreso)

#Residuales
residuales <- modelo$residuals

#Test de Kolmogorov-Smirnov
library(nortest)
prueba<-lillie.test(residuales)

#Estadístico de prueba y p-value
D <- prueba$statistic
p_value <- prueba$p.value

# Estadistico teórico para alpha=0.05
KS <- 0.895/(sqrt(length(ventas))-0.01+(0.85/sqrt(length(ventas))))

# Comparación estadístico de prueba - teórico
if(D>KS){
print("Existe evidencia suficiente para rechazar H0, por lo tanto, los residuales no provienen de una distribución normal")
} else{
  
print("No existe evidencia suficiente para rechazar H0, por lo tanto, los residuales no provienen de una distribución normal")
}

[1] "No existe evidencia suficiente para rechazar H0, por lo tanto, los residuales no provienen de una distribución normal"

# Comparación p-value - alpha
alpha <- 0.05
if(p_value<alpha){
print("Existe evidencia suficiente para rechazar H0, por lo tanto, los residuales no provienen de una distribución normal")
} else{
  
print("No existe evidencia suficiente para rechazar H0, por lo tanto, los residuales no provienen de una distribución normal")
}

[1] "No existe evidencia suficiente para rechazar H0, por lo tanto, los residuales no provienen de una distribución normal"

Verificación del supuesto de varianza constante

Verificación gráfica del supuesto de varianza constante

Una forma de verificar el supuesto de varianza constante (los términos del error tienen la misma varianza) es realizando un gráfico de dispersión con las siguientes variables:

Eje X \(\rightarrow\) ajustados o predichos \(\hat{y_{i}}\)
Eje Y \(\rightarrow\) residuales \(e_{i}\)

Si los puntos de este gráfico de dispersión se distribuyen de manera aleatoria en una banda horizontal (sin ningún patrón claro y contundente), entonces es señal de que se cumple el supuesto de que los tratamientos tienen igual varianza.

Ejemplo verificación varianza constante gráfica.

Ajuste del modelo

#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Modelo
modelo <- lm(ventas~poblacion+ingreso)

Cálculo de ajustados \(\hat{y_i}\) y residuales \(e_i\)

#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Modelo
modelo <- lm(ventas~poblacion+ingreso)

#Ajustados
ajustados <- modelo$fitted.values

#Residuales
residuales <- modelo$residuals

Gráfico de ajustados \(\hat{y_i}\) y residuales \(e_i\)

#Librerías
library(ggplot2)
library(plotly)
#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Modelo
modelo <- lm(ventas~poblacion+ingreso)

# Ajustados
ajustados <- modelo$fitted.values

#Residuales
residuales <- modelo$residuals

#Grafico
datos <- data.frame(ajustados,residuales)
minimo <- min(residuales)
maximo <- max(residuales)
grafico <- ggplot(data=datos, aes(x=ajustados, y=residuales)) + geom_point()
grafico <- grafico + geom_hline(yintercept = c(maximo+3, minimo-3))
grafico <- grafico +geom_hline(yintercept = 0, linetype="dashed", color="blue")
grafico <- grafico + xlab("Ajustados") + ylab("Residuales")
grafico <- grafico + theme_classic()
ggplotly(grafico)

Figura 2: Gráfico Ajustados vs Residuales

La Figura 2 no sugiere ningún patrón sistemático que permita sospechar sobre heterocedasticidad.

Verificación formal - analítica del supuesto de varianza constante.

El test de Breuch - Pagan para verificar la varianza constante asume que los términos del error son independientes y se distribuyen normalmente y que la varianza del término del error \(\epsilon_i\) denotada como \(\sigma^2\), está relacionada con el nivel de las \(k\) variables regresoras con \(k=1,2,...,p-1\), de la siguiente manera,

\[\begin{align} log_e \sigma_i^2 = \gamma_0 + \gamma_1x_{i1} + \gamma_2x_{i2}+...+\gamma_{p-1}x_{i~~p-1} \end{align} \tag{4}\]

1. Planteamiento de la hipótesis.

Note que Ecuación 4 implica que \(\sigma_i^2\) incrementa o disminuye con el cambio en alguna de las \(k\) variables regresoras, dependiendo del signo de \(\gamma_{p-1}\). La varianza constante se da en el caso que \(\gamma_{p-1} = 0\). Por lo anterior se pueden plantear las hipótesis de la siguiente forma:

\[\begin{align} H_0: \gamma_k &= 0 \\ H_1: \gamma_k &\ne 0~~;~~para~al~menos~un~k=1,2,...,p-1 \end{align}\]

2. Estadístico de prueba

Existen dos formas de cálculo del estadístico de prueba, el estadístico de prueba \(BP\) propuesto por los autores del test. y el estadístico de prueba \(X_{BP}^2\) propuesto por Cook y Weisberg como mejora del test.

Forma 1. Estadístico de prueba Breuch Pagan \(BP\).

El estadístico de prueba \(BP\) se calcula de la siguiente manera:

\[\begin{align} BP = nR^{2*} \end{align}\]

Donde:

\(R^{2*}\) corresponde al coeficiente de determinación de la regresión siendo \(e_i^2\) la variable respuesta y \(x_k\) la variables regresoras.

El Coeficiente de Determinación \(R^2\) se calcula de manera general de la siguiente forma:

\[\begin{align} R^2= \frac{SSR}{SST} = 1-\frac{SSE}{SST} \end{align}\]

Forma 2. Estadístico de prueba modificado Cook and Weisberg \(X_{BP}^2\).

El estadístico de prueba \(X_{BP}^2\) se calcula de la siguiente manera:

\[\begin{align} X_{BP}^2 =\frac{\frac{SSR^*}{2}}{\left[ \frac{SSE}{n} \right]^2} \end{align}\]

Donde:

\(SSR^*\) corresponde a la suma de cuadrados de la regresión siendo \(e_i^2\) la variable respuesta y \(x_i\) la variable regresora.
\(SSE\) corresponde a la suma de cuadrados del error siendo \(y_i\) la variable respuesta y \(x_i\) la variable regresora.

3. Estadístico teórico

Los estadísticos de prueba siguen una distribución Chi-Cuadrado con \(k\) grado de libertad ( \(k\) grados de libertad iguales a la cantidad de variables regresoras), esto es:

\[\begin{align} BP\sim X^2_{\alpha,~k} \end{align}\]

\[\begin{align} X_{BP}^2\sim X^2_{\alpha,~k} \end{align}\]

Rechazo \(H_0\) si:

\[\begin{align} X_{BP}^2 > X^2_{\alpha,~k} \end{align}\]

Ejemplo verificación formal-analítica Test Breuch-Pagan

Solución para estadístico \(BP\)

#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Modelo
modelo <- lm(ventas~poblacion+ingreso)

#Breuch Pagan estadistico BP
library(lmtest)
prueba_homocedasticidad<-bptest(modelo)
BP <- prueba_homocedasticidad$statistic
p_value_BP <- prueba_homocedasticidad$p.value

#Estadístico teórico para alpha=0.05 y 2 regresoras
teorico <- qchisq(0.05, 2, lower.tail = FALSE)

# Comparación estadístico de prueba-teórico
if(BP>teorico){
  print("Existe suficiente evidencia estadística para rechazar H0")} else{
  print("No existe suficiente evidencia estadística para rechazar H0")
}

[1] “No existe suficiente evidencia estadística para rechazar H0”

# Comparación p-value - alpha=0.05
alpha <- 0.05
if(p_value_BP<alpha){
  print("Existe suficiente evidencia estadística para rechazar H0")} else{
  print("No existe suficiente evidencia estadística para rechazar H0")
}

[1] “No existe suficiente evidencia estadística para rechazar H0”

Solución para estadístico \(X^2_{BP}\)

#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Modelo
modelo <- lm(ventas~poblacion+ingreso)

#Breuch Pagan estadistico BP
library(lmtest)
prueba_homocedasticidad<-bptest(modelo, studentize = FALSE)
XBP <- prueba_homocedasticidad$statistic
p_value_XBP <- prueba_homocedasticidad$p.value

#Estadístico teórico para alpha=0.05 y 2 regresoras
teorico <- qchisq(0.05, 2, lower.tail = FALSE)

# Comparación estadístico de prueba-teórico
if(XBP>teorico){
  print("Existe suficiente evidencia estadística para rechazar H0")} else{
  print("No existe suficiente evidencia estadística para rechazar H0")
}

[1] “No existe suficiente evidencia estadística para rechazar H0”

# Comparación p-value - alpha=0.05
alpha <- 0.05
if(p_value_XBP<alpha){
  print("Existe suficiente evidencia estadística para rechazar H0")} else{
  print("No existe suficiente evidencia estadística para rechazar H0")
}

[1] “No existe suficiente evidencia estadística para rechazar H0”

Verificación del supuesto de independencia

Verificación gráfica del supuesto de independencia.

La graficación de los residuales en orden temporal, de obtención o recolección de datos, en los casos en que estos sean tomados de forma experimental, es útil para detectar correlaciones entre los residuales \(e_i\). Una tendencia identificable en el gráfico de residuales en orden temporal indica una correlación entre ellos. Esto implicaría que el supuesto de independencia de los errores \(\epsilon_i\) ha sido violado.

Una forma gráfica de verificar el supuesto de independencia es realizando un gráfico del siguiente tipo:

Eje X \(\rightarrow\) secuencia de tiempo desde \(t=1\) hasta \(t=n\).
Eje Y \(\rightarrow\) residuales \(e_{i}\)

Cuando los términos de error son independientes, esperamos que los residuos en un gráfico de secuencia fluctúen en un patrón más o menos aleatorio alrededor de la línea base.

Ejemplo verificación independencia gráfica.

Ajuste del modelo

#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Modelo
modelo <- lm(ventas~poblacion+ingreso)

Cálculo de residuales \(e_i\)

#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Modelo
modelo <- lm(ventas~poblacion+ingreso)

#Residuales
residuales <- modelo$residuals

Gráfico de residuales \(e_i\) en secuencia temporal

#Librerías
library(ggplot2)
library(plotly)
#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Modelo
modelo <- lm(ventas~poblacion+ingreso)

# Ajustados
ajustados <- modelo$fitted.values

#Residuales
residuales <- modelo$residuals

# Secuencia temporal
t <- c(1:length(ventas))
#Grafico
datos <- data.frame(t,residuales)
minimo <- min(residuales)
maximo <- max(residuales)
grafico <- ggplot(data=datos, aes(x=t, y=residuales)) + geom_point()
grafico <- grafico + geom_hline(yintercept = c(maximo+3, minimo-3))
grafico <- grafico +geom_hline(yintercept = 0, linetype="dashed", color="blue")
grafico <- grafico + xlab("Tiempo") + ylab("Residuales")
grafico <- grafico + theme_classic()
ggplotly(grafico)

Figura 3: Gráfico residuales en secuencia temporal

La Figura 3 no sugiere ningún patrón sistemático que permita sospechar sobre violación del supuesto de independencia.

Verificación formal - analítica supuesto de independencia mediante Test de Durbin - Watson.

La verificación formal-analítica del supuesto de independencia se realiza de manera indirecta a través del Test de Durbin-Watson. Esta prueba permite diagnosticar la presencia de correlación (autocorrelación) entre los residuales consecutivos (ordenados en el tiempo), que es una posible manifestación de la falta de independencia. Es importante mencionar que si se ha comprobado la normalidad y se verifica la autocorrelación implica entonces independencia.

1. Planteamiento de la hipótesis.

Sea \(\rho\) el parámetro que representa la correlación entre residuales consecutivos (ordenados en secuancia temporal), es decir, \(\rho = corr(e_t,~e_{t+1})\). Las hipótesis en la prueba de Durbin-Watson se plantearían de la siguiente manera:

\[\begin{align} H_0&: \rho=0\\ H_1&: \rho \ne 0 \end{align}\]

2. Estadístico de prueba \(d_0\) y comparación.

El estadístico de prueba \(d_0\) para la prueba de Durbin - Watson se calcula como sigue:

\[\begin{align} d_0 =\frac{\sum_{t=2}^T (e_t - e_{t-1})^2}{\sum_{t=1}^T e_t^2}~;~~ \forall~~ t=1,2,..., T \end{align}\]

\(e_t\) corresponde a los residuales \(e_i\) en secuencia temporal.

Se compara el estadístico de prueba con límites teóricos \(d_L\) y \(d_U\) usando las reglas expuestas a continuación:

\[\begin{align} &Si~ d_0 < d_L~\rightarrow~Se~rechaza~H_0\\ &Si~ d_L<d_0<d_U ~\rightarrow~Prueba~no~concluyente\\ &Si~ d_U<d_0<4-d_U ~\rightarrow~No~se~rechaza~H_0\\ &Si~ 4-d_U<d_0<4-d_L ~\rightarrow~Prueba~no~concluyente\\ &Si~ d_0 > 4-d_L~\rightarrow~Se~rechaza~H_0\\ \end{align}\]

Ejemplo verificación independencia Test de Durbin - Watson.

Ajuste del modelo

#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Modelo
modelo <- lm(ventas~poblacion+ingreso)

Cálculo de residuales \(e_i\)

#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Modelo
modelo <- lm(ventas~poblacion+ingreso)

#Residuales
residuales <- modelo$residuals

Cálculo de estadístico \(d_0\)

#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Modelo
modelo <- lm(ventas~poblacion+ingreso)

#Residuales
residuales <- modelo$residuals

# Estadístico de prueba
d_0 <- sum((residuales[2:length(residuales)]-residuales[1:(length(residuales)-1)])^2)/sum(residuales^2)

Estadísticos teóricos \(d_L\) y \(d_U\) y comparación

#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Modelo
modelo <- lm(ventas~poblacion+ingreso)

#Residuales
residuales <- modelo$residuals

# Estadístico de prueba
d_0 <- sum((residuales[2:length(residuales)]-residuales[1:(length(residuales)-1)])^2)/sum(residuales^2)

# Estadísticos teóricos para dos regresores y 21 observaciones
dL <- 1.125
dU <- 1.538
cuatro_dL <- 4-dL
cuatro_dU <- 4-dU
# Comparación
if(d_0<dL){
  print("Se rechaza H0")
} else if (dL <d_0 & d_0< dU){
  print("Prueba no concluyente")
} else if (dU <d_0 & d_0< cuatro_dU){
  print("No se rechaza H0")
} else if (cuatro_dU <d_0 & d_0< cuatro_dL){
  print("Prueba no concluyente")
} else if (d_0> cuatro_dL){
  print("Se rechaza H0")
}

[1] “No se rechaza H0”

Se puede realizar la prueba usando sentencias de R

#Datos
ventas <- c(174.4,164.4,244.2,154.6,181.6,207.5,152.8,163.2,145.4,137.2,241.9,191.1,232.0,145.3,161.1,209.7,146.4,144.0,232.6,224.1,166.5)
poblacion <- c(68.5,45.2,91.3,47.8,46.9,66.1,49.5,52.0,48.9,38.4,87.9,72.8,88.4,42.9,52.5,85.7,41.3,51.7,89.6,82.7,52.3)
ingreso <- c(16.7,16.8,18.2,16.3,17.3,18.2,15.9,17.2,16.6,16,18.3,17.1,17.4,15.8,17.8,18.4,16.5,16.3,18.1,19.1,16)

#Modelo
modelo <- lm(ventas~poblacion+ingreso)

#Residuales
residuales <- modelo$residuals

#Prueba
library(car)
prueba<-durbinWatsonTest(modelo, alternative = "two.sided")
d_0 <- prueba$dw
d_0

[1] 1.653144

De la prueba se obtiene que \(d_0 = 1.653144\). De tabla, para \(2\) variables regresoras y \(n=21\) se obtiene \(d_L = 1.125\) y \(d_U = 1.538\). Por lo tanto, \(d_U<d_0<4-d_U\), se concluye que, no existe evidencia suficiente para rechazar \(H_0\) por lo que \(\rho = 0\), los residuales \(e_i\) no están autocorrelacionados. Sabiendo que se verificó el supuesto de normalidad, se concluye que se cumple el supuesto de independencia.