Introducción

En esta unidad se iniciará con el concepto de correlación lineal para variables cuantitativas, luego se revisarán los conceptos y métodos de estimación de la regresión lineal simple y el modelo múltiple, la inferencia sobre el modelo, interpretación de los coeficientes y la validación de supuestos del modelo.

Para el ejercicio práctico, se trabajá con la librería stats y una primera base de datos que tiene datos simulados sobre: porcentaje de personas con enfermedad cardiaca (heart.disease), porcentaje de personas que usa la bicicleta como medio de transporte (biking) y porcentaje de personas que fuman (smoking) en una muestra de ciudades independientes.

1. Correlación lineal

Cuando se quiere cuantificar la relación que tienen dos variables cuantitativas existen varios indicadores que permiten establecer si existe o no correlación y la magnitud de esta. El indice de correlación de Pearson, determina que tan relacionadas están dos variables en términos lineales, es decir, que tan alineada es su relación respecto a una recta. Este índice se calcula a través de la covarianza, la cual mide que tan alejadas están el par de observaciones en \(x\) y \(y\) respecto a sus respectivos promedios, sin embargo, esta es una medida que arroja valores en las unidades de medida de las variables, es por esto que se realiza una estandarización a través de la desviación estándar de las variables \(x\) y \(y\). El indice se calcula de la siguiente forma:

\(r_{xy}=\frac{S_{xy}}{S_{x}S_{y}},\)

donde: \(S_{xy}\) es la covarianza entre \(x\) y \(y\), \(S_{x}\) y \(S_{y}\) es la desviación estándar de \(x\) y \(y\), respectivamente. Este indicador tiene una escala entre -1 y 1, entre más cercano a -1 o 1, indica que la relación es fuerte y está alineada exactamente sobre una recta, si es negativa la relación es inversa (recta con pendiente negativa), si es positiva (recta con pendiente positiva).

2. Regresión Lineal simple

Tenemos un par de observaciones \((x_i, y_i)\), para \(i = 1, 2, \ldots n\), donde \(n\) es el tamaño de la muestra seleccionada. El indice \(i\) se utiliza por notación, para determinar que esa variable denotada por la letra \(x\) o \(y\) toma diferentes valores de acuerdo a la variación de \(i\). A \(x_i\) la llamaremos la variable predictora o explicativa, la cual ayuda a predecir o explicar de la respuesta o variable predicha \(y_i\). En el caso del modelo simple, se tiene una sola variable explicativa, pero luego se ampliará a un modelo múltiple donde se tienen varias \(x\).

Al modelar la relación entre \(X\) e \(Y\), se utilizará una función de la siguiente forma:

\[ Y = f(X) + \epsilon. \]

donde \(f\) describe la relación funcional entre dos variables y \(\epsilon\) determina el error, dado que este es un modelo estocástico, no exacto, es decir, al ingresar un valor de \(x\) al modelo, este estima el valor \(y\) con cierta cantidad de error. La función \(f\) puede ser de cuaquier forma, en este módulo se va a restringuir a una función lineal. El modelo queda de la siguiente forma:

\[ Y = \beta_0 + \beta_1 X + \epsilon. \]

donde \(\beta_0\) es el intercepto y \(\beta_1\) es la pendiente. Teniendo en cuenta además que,

\[ \epsilon_i \sim N(0, \sigma^2). \] Es decir, que \(\epsilon_i\) son variables aleatorias normales independientes e idénticamente distribuidas con media \(0\) y varianza \(\sigma^2\). Este modelo tiene 3 parámetros a ser estimados: \(\beta_0\), \(\beta_1\), y \(\sigma^2\).

Nota: En términos de notación, tener en cuenta que las letras en mayúscula indican variables aleatorias, \(Y\), y minúscula un valor potencial para una variable aleatoria. Cómo tenemos \(n\) observaciones, vamos a tener \(n\) variables aleatorias \(Y_i\) y sus posibles valores \(y_i\).

La variable respuesta \(Y_i\) es una variable aleatoria debido al comportamiento aleatorio del error \(\epsilon_i\), cada respuesta \(Y_i\) está enlazada a un valor observable \(x_i\) y a un valor aleatorio no observable \(\epsilon_i\). Es decir, \(Y_i\) tiene una distribución condicional dependiente del valor de of \(X_i\).

\[ Y_i \mid X_i \sim N(\beta_0 + \beta_1 x_i, \sigma^2) \] Cómo \(Y_i\) es una función de \(x_i\), entonces podemos escribir su media (valor esperado) como una función de \(x_i\),

\[ \text{E}[Y_i \mid X_i = x_i] = \beta_0 + \beta_1 x_i. \]

Sin embargo, su varianza se mantiene constante,

\[ \text{Var}[Y_i \mid X_i = x_i] = \sigma^2. \] En la siguiente imagen podemos ver que para un valor dado \(x\), el valor esperado de \(Y\) es \(\beta_0 + \beta_1 x\). Cada valor de \(x\), \(Y\) tienen la misma varianza \(\sigma^2\).

Supuestos del modelo de regresión Lineal

  • La relación entre \(Y\) y \(x\) es lineal con la forma \(\beta_0 + \beta_1 x\).
  • Los errores \(\epsilon\) son independientes.
  • Los errores \(\epsilon\) se distribuyen normal.
  • Para cada valor de \(x\), la varianza de \(Y\) es la misma, \(\sigma^2\).

Además, no se supone ninguna distribución de probabilidad sobre \(X\) dado que se asume que son fijos, no aletorios.

2.1 Métodos de estimación para los coeficientes

Los modelos de regresión lineal (RL), definen a \(Y\) como una función de \(X\), pero la pregunta que surge es cómo definimos una buena linea recta?, dado que por un conjunto de puntos \((x,y)\) pueden pasar varias rectas. La respuesta está en tratar de encontrar la recta que tenga la menor suma de errores, entendiendo el error cómo la diferencia entre los valores observados y los estimados o predichos por el modelo \((y_i - \hat{y})\), ver imagen.

Mínimos Cuadrados Ordinarios

Para determinar los valores para los \(betas\) en el modelo de regresión, existen varios métodos de estimación. En regresión lineal uno de los más usados es el de mínimos cuadrados ordinarios OLS, el cual busca estimar los parámetros tal que la suma de cuadrados de la diferencias entre las observaciones y la línea recta estimada sea mínima. Este método en RL permite obtener estimadores insesgados y de mínima varianza para los \(betas\).

dadas las observaciones \((x_i, y_i)\), para \(i = 1, 2, \ldots n\), queremos encontrar los valores de \(\beta_0\) y \(\beta_1\), es decir, \(\hat{\beta}_0\) y \(\hat{\beta}_1\) que minimicen

\[ f(\beta_0, \beta_1) = \sum_{i = 1}^{n}(y_i - (\beta_0 + \beta_1 x_i))^2 = \sum_{i = 1}^{n}(y_i - \beta_0 - \beta_1 x_i)^2. \]

Resolviendo el sistema de ecuaciones normales:

\[ \begin{aligned} n \beta_0 + \beta_1 \sum_{i = 1}^{n} x_i &= \sum_{i = 1}^{n} y_i\\ \beta_0 \sum_{i = 1}^{n} x_i + \beta_1 \sum_{i = 1}^{n} x_i^2 &= \sum_{i = 1}^{n} x_i y_i \end{aligned} \] Se obtienen las siguientes expresiones para \(\hat{\beta}_0\) y \(\hat{\beta}_1\):

\[ \begin{aligned} \hat{\beta}_1 &= \frac{\sum_{i = 1}^{n} x_i y_i - \frac{(\sum_{i = 1}^{n} x_i)(\sum_{i = 1}^{n} y_i)}{n}}{\sum_{i = 1}^{n} x_i^2 - \frac{(\sum_{i = 1}^{n} x_i)^2}{n}} = \frac{S_{xy}}{S_{xx}}\\ \hat{\beta}_0 &= \bar{y} - \hat{\beta}_1 \bar{x} \end{aligned} \]

los parámetros \(\beta_0\) y \(\beta_1\), indican el valor promedio de \(Y\) cuando \(x\) toma el valor cero y el cambio en \(Y\) por cada unidad de aumento en \(X\).

Descomposición de la varianza

En un modelo de regresión, la variación total se puede expresar en términos de la variación debida la regresión y al error, a través de la siguiente expresión:

\[ y_i - \bar{y} = (y_i - \hat{y}_i) + (\hat{y}_i - \bar{y}). \] Aqui, \(y_i - \hat{y}_i\) mide la desviación de una observación respecto a la linea de regresión y \(\hat{y}_i - \bar{y}\) mide la desviación de la línea de regresión ajustada del promedio muestral. Al elevar al cuadrado y agregar la sumatoria a ambos lados de la ecuación anterior:

\[ \sum_{i=1}^{n}(y_i - \bar{y})^2 = \sum_{i=1}^{n}(y_i - \hat{y}_i)^2 + \sum_{i=1}^{n}(\hat{y}_i - \bar{y})^2. \]

Lo anterior, muestra que la suma de cuadrados totales \(SST\) la cual representa la variación total de los valores observados \(y\), es igual a la suma de cuadrados de la regresión \(SSReg\) más la suma de cuadrados del error \(SSE\). Estas dos última muestran la avriación explicada por los valores observados \(y\) y la variación no explicada por los valores de \(y\).

Coeficiente de determinación

El coeficiente de determinación \(R^2\) está definido como:

\[ \begin{aligned} R^2 &= \frac{\text{SSReg}}{\text{SST}} = \frac{\sum_{i=1}^{n}(\hat{y}_i - \bar{y})^2}{\sum_{i=1}^{n}(y_i - \bar{y})^2} \\[2.5ex] &= \frac{\text{SST} - \text{SSE}}{\text{SST}} = 1 - \frac{\text{SSE}}{\text{SST}} \\[2.5ex] &= 1 - \frac{\sum_{i=1}^{n}(y_i - \hat{y}_i)^2}{\sum_{i=1}^{n}(y_i - \bar{y})^2} = 1 - \frac{\sum_{i = 1}^{n}e_i^2}{\sum_{i=1}^{n}(y_i - \bar{y})^2} \end{aligned} \] este se interpreta como la proporción de varianza de \(y\) que es explicada por el modelo.

2.2 Inferencia de la regresión lineal

Para realizar inferencia sobre los parámetros del modelo \(\beta_0\) y \(\beta_1\), se deben calcular los errores estándar y estadísticos de prueba como los vistos en el curso anterior cuando se realizaban sobre el promedio, la proporción o la varianza, pero en este caso lo que queremos probar es si los betas son significativos (diferentes de cero) en el modelo. De esta forma es importante determinar el error estándar de cada coeficiente. Teniendo en cuenta que los beta´s son variables,

\[ \begin{aligned} \hat{\beta}_0 &\sim N\left( \beta_0, \sigma^2 \left(\frac{1}{n} + \frac{\bar{x}^2}{S_{xx}}\right) \right) \\ \hat{\beta}_1 &\sim N\left( \beta_1, \frac{\sigma^2}{S_{xx}} \right). \end{aligned} \]

estandarizando, se obtiene que

\[ \frac{\hat{\beta}_0 - \beta_0}{\text{SD}[\hat{\beta}_0]} \sim N(0, 1) \] \[ \frac{\hat{\beta}_1 - \beta_1}{\text{SD}[\hat{\beta}_1]} \sim N(0, 1) \] De esta forma, los errores estándar por cada coeficienten sería:

\[ \text{SE}[\hat{\beta}_0] = s_e\sqrt{\frac{1}{n} + \frac{\bar{x}^2}{S_{xx}}} \]

\[ \text{SE}[\hat{\beta}_1] = \frac{s_e}{\sqrt{S_{xx}}} \]

Recordemos que el error estándar es la desviación estándar estimada en la distribución muestral de cada beta, teniendo en cuenta que, provienen de una muestra aleatoria. Diviendo por el error estándar, se obtienen los siguientes resultados que permiten estimar intervalos de confianza y contrastar hipótesis usando una distribución \(t\):

\[ \frac{\hat{\beta}_0 - \beta_0}{\text{SE}[\hat{\beta}_0]} \sim t_{n-2} \]

\[ \frac{\hat{\beta}_1 - \beta_1}{\text{SE}[\hat{\beta}_1]} \sim t_{n-2} \]

Finalmente, para los procesos de estimación mencionados anteriormente, se utilizan las siguientes expresiones:

\[ \hat{\beta}_0 \pm t_{\alpha/2, n - 2} \cdot \text{SE}[\hat{\beta}_0] \quad \quad \quad \hat{\beta}_0 \pm t_{\alpha/2, n - 2} \cdot s_e\sqrt{\frac{1}{n}+\frac{\bar{x}^2}{S_{xx}}} \]

\[ \hat{\beta}_1 \pm t_{\alpha/2, n - 2} \cdot \text{SE}[\hat{\beta}_1] \quad \quad \quad \hat{\beta}_1 \pm t_{\alpha/2, n - 2} \cdot \frac{s_e}{\sqrt{S_{xx}}} \] donde \(t_{\alpha/2, n - 2}\) es el valor crítico tal que \(P(t_{n-2} > t_{\alpha/2, n - 2}) = \alpha/2\).

3. Regresión lineal múltiple

Ampliando los conceptos revisados para la regresión simple, el modelo múltiple tiene en cuenta más de una variable explicativa o covariable, quedando de la siguiente forma:

\[ Y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} +....+\beta_j x_{ij} + \epsilon_i, \qquad i = 1, 2, \ldots, n \]

donde \(\epsilon_i \sim N(0, \sigma^2)\). El objetivo ahora es obtener la distribución del vector \(\hat{\beta}\)

\[ \hat{\beta} = \begin{bmatrix} \hat{\beta}_0 \\ \hat{\beta}_1 \\ \hat{\beta}_2 \\ \vdots \\ \hat{\beta}_{p-1} \end{bmatrix} \] Teniendo en cuenta lo mencionado en el modelo simple, \(\hat{\beta}\) es una variable aleatoria, ahora sería un vector aleatorio. Cada uno de los \(\hat{\beta}_j\) sigue una distribución normal,

\[ \hat{\beta}_j \sim N\left(\beta_j, \sigma^2 C_{jj} \right). \] donde \[C = \left(X^\top X\right)^{-1}\] Con base en esto la distribución de probabilidad para los beta´s sería entonces:

\[ \frac{\hat{\beta}_j - \beta_j}{s_e \sqrt{C_{jj}}} \sim t_{n-p}. \]

Teniendo en cuenta lo anterior, se pueden calcular intervalos de confianza para cada beta y probar hipótesis respecto a su efecto en el modelo de regresión de la siguiente forma:

\[ H_0: \beta_j = 0 \quad \text{vs} \quad H_1: \beta_j \neq 0 \]

donde el estadístico de prueba:

\[ \text{TS} = \frac{\text{EST} - \text{HYP}}{\text{SE}}. \]

Para este caso:

\[ t = \frac{\hat{\beta}_j - \beta_j}{\text{SE}[\hat{\beta}_j]} = \frac{\hat{\beta}_j-0}{s_e\sqrt{C_{jj}}}, \]

Bajo la hipótesis nula, sigue una distribución \(t\) con \(n-p\) grados de libertad.

4. Aplicación del modelo de Regresión Lineal Simple y Múltiple

Ahora revisaremos los conceptos mencionados anteriormente a través de una ejercicio aplicado. Para revisar cómo se estima el modelo en \(R\), se inciará cargando las siguientes librerías:

library(ggplot2)
library(dplyr)
library(stats)

Se exporta la base de datos que está en .csv, cómo ya se mencionó esta contiene datos sobre: porcentaje de personas con enfermedad cardiaca (heart.disease), porcentaje de personas que usa la bicicleta como medio de transporte (biking) y porcentaje de personas que fuman (smoking), en una muestra de 498 ciudades independientes.

Inicalmente se realizará un análisis exploratorio descriptivo de las variables en la base de datos.

summary(base_heart)
##        id            biking          smoking        heart.disease    
##  Min.   :  1.0   Min.   : 1.119   Min.   : 0.5259   Min.   : 0.5519  
##  1st Qu.:125.2   1st Qu.:20.205   1st Qu.: 8.2798   1st Qu.: 6.5137  
##  Median :249.5   Median :35.824   Median :15.8146   Median :10.3853  
##  Mean   :249.5   Mean   :37.788   Mean   :15.4350   Mean   :10.1745  
##  3rd Qu.:373.8   3rd Qu.:57.853   3rd Qu.:22.5689   3rd Qu.:13.7240  
##  Max.   :498.0   Max.   :74.907   Max.   :29.9467   Max.   :20.4535

Para observar la distribución de la variable respuesta: porcentaje de personas con enfermedad cardiaca, se realizará un histograma:

hist(base_heart$heart.disease, main="Histograma % personas con enfermedad cardíaca")

Utilizaremos la función \(cor\), la cual permite calcular la correlación de pearson u otros métodos como kendall o spearman, estos últimos cuando no se cumple la normalidad.

Se prueba primero la normalidad de \(y\) a través de la prueba de Shapiro-Wilk.

shapiro.test(base_heart$heart.disease)
## 
##  Shapiro-Wilk normality test
## 
## data:  base_heart$heart.disease
## W = 0.98047, p-value = 3.158e-06

Debido a que no se cumple, se calcula la correlación usando un método no paramétrico como kendall.

cor(base_heart$heart.disease, base_heart$biking, method = c("kendall"))
## [1] -0.7745186
# en método se puede seleccionar: pearson, kendall o spearman.
cor(base_heart, method = c("kendall"))
##                        id      biking    smoking heart.disease
## id             1.00000000  0.03806776 0.03713041   -0.03435068
## biking         0.03806776  1.00000000 0.00698973   -0.77451860
## smoking        0.03713041  0.00698973 1.00000000    0.19845175
## heart.disease -0.03435068 -0.77451860 0.19845175    1.00000000

También se pueden construir métodos gráficos para observar la relación entre variables de la base de datos.

#install.packages("PerformanceAnalytics")
library("PerformanceAnalytics")
chart.Correlation(base_heart, histogram=TRUE, pch=19)

Los asteriscos es el gráfico anterior muestran la significancia de la correlación a través de valores p que están asociados con los simbolos: p-values(0, 0.001, 0.01, 0.05, 0.1, 1) <=> simbolos(“”, “”, “”, “.”, ” “).

También se pueden realizar correlogramas usando la librería \(corrplot\).

#install.packages("corrplot")
library(corrplot)
## corrplot 0.92 loaded
cor=cor(base_heart, method = c("kendall"))
corrplot(cor, type = "upper", order = "hclust", 
         tl.col = "black", tl.srt = 45)

Aunque a través de la función chart.correlation se puede observar gráficamente la relación entre variables, también se puede hacer uso de la función \(plot\)

par(mfrow=c(1,2))  #permite tener una ventana con una fila y dos columnas para los gráficos.
plot(heart.disease ~ biking, data=base_heart)
plot(heart.disease ~ smoking, data=base_heart)

Para la estimación del modelo lineal simple se usa la función \(lm\). Ajustaremos el siguiente modelo: Heart_disease = \(\beta_0\) + \(\beta_2\) smoking

model<-lm(heart.disease ~ smoking, data = base_heart)
summary(model)
## 
## Call:
## lm(formula = heart.disease ~ smoking, data = base_heart)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.7065 -3.7069  0.5007  3.6597  8.5434 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7.54311    0.41251  18.286  < 2e-16 ***
## smoking      0.17048    0.02355   7.239 1.73e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.352 on 496 degrees of freedom
## Multiple R-squared:  0.09556,    Adjusted R-squared:  0.09374 
## F-statistic: 52.41 on 1 and 496 DF,  p-value: 1.729e-12

Tenieno en cuenta la salida anterior, vamos a iniciar con el intercepto \(\beta_0\), este nos indica el valor promedio de \(y\) cuando las variables explicativas son iguales a 0, es decir, el valor esperado del porcentaje de personas con enfermedad cardíaca es 7.54 cuando el porcentaje de quienes fuman es cero. En el caso de \(\beta_1\), este nos indica que por cada unidad de aumento en el porcentaje de personas que fuman, el porcentaje de personas con enfermedad cardíaca aumenta en 0.17, con un valor p menor al nivel de significancia de 0.05, es decir, es una variable significativa en la explicación o predicción de Y (porcentaje de personas con enfermedad cardíaca).

Asi mismo, para la estimación del modelo lineal múltiple, se usa la función \(lm\). Ajustaremos el siguiente modelo: Heart_disease = \(\beta_0\) + \(\beta_1\) biking + \(\beta_2\) smoking

modelo<-lm(heart.disease ~ biking + smoking, data = base_heart)
summary(modelo)
## 
## Call:
## lm(formula = heart.disease ~ biking + smoking, data = base_heart)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1789 -0.4463  0.0362  0.4422  1.9331 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 14.984658   0.080137  186.99   <2e-16 ***
## biking      -0.200133   0.001366 -146.53   <2e-16 ***
## smoking      0.178334   0.003539   50.39   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.654 on 495 degrees of freedom
## Multiple R-squared:  0.9796, Adjusted R-squared:  0.9795 
## F-statistic: 1.19e+04 on 2 and 495 DF,  p-value: < 2.2e-16

La salida anterior muestra un análisis descriptivo de los residuales del modelo, los coeficientes beta´s, el error estándar, la prueba t para cada coediciente y los valores p. Además la estimación de la prueba de bondad y ajuste del modelo \(R^2\) y la prueba F, del análisis de varianza (anova), que plantea como hipótesis nula que todos los beta´s son iguales a cero.

Vamos a iniciar la interpretación con el intercepto \(\beta_0\), este nos indica el valor promedio de \(y\) cuando las variables explicativas toman el valor de 0, es decir, sería que el valor esperado del porcentaje de personas con enfermedad cardíaca es 14.98 cuando el porcentaje de quienes usan bicileta es cero y de los que fuman igual. En el caso de \(\beta_1\), este nos indica que por cada unidad de aumento en el porcentaje de personas que usan la bicicleta, el porcentaje de personas con enfermedad cardíaca disminuye en 0.20 independiente del efecto de las demás covariables presentes en el modelo, es decir, independiente del porcentaje de personas que fuman. Finalmente para \(\beta_2\), este indica que por por cada unidad de aumento en el porcentaje de personas que fuman en una ciudad, el porcentaje de personas con enfermedad cardíaca aumenta en 0.1,independiente del porcentaje de personas que usan bicicleta.

En cuanto a la significancia estadística, el valor p representado a través de la columna de probabilidad \(Pr(>|t|)\), debido a que son valores muy pequeños, indica que se rechaza la hipótesis nula \(H_0: \beta_j = 0\), por lo tanto, podemos decir que las covariables: porcentaje de personas que usan bicileta y porcentaje de personas que fuma es significativo en el modelo.

Además, podemos decir que, el porcentaje de variabilidad de \(y\) explicado por el modelo es 97.9% que corresponde a la medida de bondad de ajuste \(R^{2}\).

La estimación de intervalos de confianza para los coeficientes se obtienen de la siguiente forma:

confint(modelo, level = 0.95)
##                  2.5 %     97.5 %
## (Intercept) 14.8272075 15.1421084
## biking      -0.2028166 -0.1974495
## smoking      0.1713800  0.1852878

Si buscamos predecir valores de \(y\) para valores específicos de las \(x´s\), se crean primero los valores de x y se utiliza la función predict. Por ejemplo, nos interesa conocer el valor predicho del porcentaje de personas con enfermedad cardíaca para una ciudad donde el porcentaje de biking sea 0.50 y el de smoking 0.60, asi mismo, para una ciudad con 0.70 y 0.80, en cada covariable.

new_data = data.frame(biking = c(0.50, 0.70), smoking = c(0.60, 0.80)) # se deben tener en cuenta los mismos nombres de las variables en el dataframe.
new_data
##   biking smoking
## 1    0.5     0.6
## 2    0.7     0.8

luego se hace la predicción para estos nuevos datos:

predict(modelo, newdata = new_data, interval = "confidence", level = 0.95)
##        fit      lwr      upr
## 1 14.99159 14.83778 15.14541
## 2 14.98723 14.83468 15.13979

Para otros criterios de bondad de ajuste del modelo, se pueden usar el criterio akaike (AIC) o el criterio de información bayesiano (BIC), los cuales castigan por la cantidad de variables explicativas en el modelo buscando cumplir con el principio de parsimonia, donde lo que se busca es que el modelo explique \(y\) pero que sea sencillo. Dado que el \(R^{2}\) tiende a aumentar de acuerdo a las covariables ingresadas, en algunos casos es mejor estas dos métricas. Sin embargo, se utilizan más cuando se comparan modelos, y entre menor sea el valor, mejor el ajuste del modelo:

AIC(modelo)
## [1] 995.3534
BIC(modelo)
## [1] 1012.196

5. Validación de supuestos

Para este modelo de regresión lineal se deben validar los supuestos respecto a los residuales como la normalidad, homocedasticidad de la varianza e independencia.

El siguiente gráfico, nos muestra diferentes plots que ayudan con la validación de supuestos. El primero Residuals vs Fitted se utiliza para revisar los supuestos de linealidad y varianza constante, lo que se espera es que para cualquier valor ajustado (eje x) los residuales se encuentren alrededor de cero. El Q-Q plot es otro método visual que nos permite evaluar la normalidad de los errores, comparando los cuantiles de una distribución teórica normal con la distribución de los residuales, lo que se espera es que los puntos estén sobre la linea recta, que sería una distribución normal.

En el caso del gráfico Scale location, este también permite evaluar la homocedasticidad, lo que se espera es no encontrar ninguna tendencia en el comportamineto de los valores ajustado versus los residuales estandarizados. Por último, el gráfico de Residuales vs Leverage permite determinar si existen observaciones que pueden influir en las estimaciones del modelo, a través de la distancia de Cook, esta es una medida de la distancia cuadrática entre, el estimador por mínimos cuadrados basado en las n observaciones y el estimador obtenido eliminando la i-ésima observación. En este gráfico se identifican aquellos puntos (observaciones) que estén muy alejados de cero, se identifican y revisan para determinar que está pasando con estos datos.

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

El gráfico anterior permite evaluar la tendencia de residuales, la normalidad (QQ plot) y posibles valores atípicos.

Para determinar homocedasticidad de la varianza, se puede aplicar la prueba analítica Breush-Pagan (Hipótesis nula: residuales distribuidos con igual varianza).

library(lmtest)
lmtest::bptest(modelo)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo
## BP = 5.7775, df = 2, p-value = 0.05564

Esta prueba plantea que los residuales se distribuyen con la misma varianza, es decir, hay homocedasticidad, debido a que el valor p es mayor por ejemplo a 0.05 (significancia) vamos a decir que no se rechaza la hipótesis nula y se cumple, sin embargo, se debe tener cuidado dado que es un valor cercano al limite establecido por el error.

Para la normalidad de los residuales, primero calculamos los residuales del modelo usando la función resid, luego se aplica la prueba de shapiro que plantea como hipótesis nula la normalidad. En este caso, se observa que teniendo en cuenta el valor p no se rechaza la hipótesis nula que plantea noramilidad en los errores.

residuales=resid(modelo)
ks.test(residuales, 'pnorm')
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  residuales
## D = 0.11967, p-value = 1.278e-06
## alternative hypothesis: two-sided
shapiro.test(residuales)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales
## W = 0.997, p-value = 0.4935

5.1. Identificación de observaciones extremas

Además de la validación de supuestos en los errores de un modelo de RLM, se debe chequear la presencia de observaciones extremas, tales como: observaciones atípicas (outliers), puntos de balanceo y observaciones influyentes.

Una observación atípica (o outlier) es aquella que está separada (en su valor de la respuesta \(y\)) del resto de las observaciones y por tanto puede afectar los resultados del ajuste del modelo de regresión.

Interesa identificarlas para luego, si es posible analizar si se tratan de observaciones malas (por errores de registro o medición) que pueden ser descartadas, o si realmente son datos correctos pero extraños que no deben ser eliminados del conjunto de datos.

Para detectar observaciones atípicas se usan los residuales escalados, los cuales estandarizan los residuales para tener media cero y varianza igual a uno. Se considera que una observación es atípica cuando su residual \(r_{i}\) , es tal que:\(|r_i| > 3\)

# Residuales estandarizados
resid_std <- rstandard(modelo)

# Plot
plot(resid_std, type = "h", main = "Residuales estandarizados",
     ylab = "Residuales estandarizados", xlab = "Observación")
abline(h = c(-3, 3), col = "red", lty = 2) # umbral típico |r_i| > 3
# Etiquetas
outliers_res <- which(abs(resid_std) > 3)
text(outliers_res, resid_std[outliers_res], labels=outliers_res, pos=4, col="blue")

Un punto de balanceo o apalancamiento es una observación en el espacio de las predictoras, alejada del resto de la muestra y que puede controlar ciertas propiedades del modelo ajustado.

Este tipo de observaciones posiblemente no afecte los coeficientes de regresión estimados pero sí las estadísticas de resumen como el \(R^{2}\) y los errores estándar de los coeficientes estimados. Los puntos de balanceo son detectados mediante el análisis de los elementos de la diagonal principal de la matriz \(H\), los \(h_{i}\), que proporcionan una medida estandarizada de la distancia de la i-ésima observación al centro del espacio definido por las predictoras.

Se asume que la observación \(i\) es un punto de balanceo si \(h_{i} > 2p/n\), dodne \(p\) es el número de parámetros y \(n\) el número de observaciones.Pero si \(2p/n > 1\) este criterio no funciona pues los \(h_{i}\) siempre son menores que 1. Observaciones con \(h_{i}\) grandes y residuales \(r_{i}\) también grandes probablemente serán influyentes.

hat_values <- hatvalues(modelo)

plot(hat_values, type = "h", main = "Puntos de apalancamiento (Leverage)",
     ylab = "Leverage", xlab = "Observación")
abline(h = 2*mean(hat_values), col = "red", lty = 2)
# Etiquetas para puntos que exceden el umbral
outliers_lev <- which(hat_values > 2*mean(hat_values))
text(outliers_lev, hat_values[outliers_lev], labels=outliers_lev, pos=4, col="blue")

Una observación es influyente si tiene un impacto notable sobre los coeficientes de regresión ajustados, esto es, una observación influyente se dice hala al modelo en su dirección, es decir, si su exclusión del modelo causa cambios importantes en la ecuación de regresión ajustada.

Estas observaciones se caracterizan por tener un valor moderadamente inusual tanto en el espacio de las predictoras como en la respuesta.

Después de identificar las observaciones que están alejadas con respecto a los valores de \(Y\) (atípicas) y/o con respecto a sus valores en \(X\) (puntos de balanceo) evaluamos si éstas son influenciales.

Para la evaluación se cuenta con medidas como la distancia de Cook, donde, \(D_{i} > 4/n\) se considera una observación influyente.

# Distancia de Cook
cook <- cooks.distance(modelo)

# Plot
plot(cook, type = "h", main = "Distancia de Cook",
     ylab = "Cook's distance", xlab = "Observación")
abline(h = 4/length(cook), col = "red", lty = 2) # umbral sugerido
# Etiquetas
outliers_cook <- which(cook > 4/length(cook))
text(outliers_cook, cook[outliers_cook], labels=outliers_cook, pos=4, col="blue")

Finalmente, el diagnóstico de puntos influyentes o outliers, permiten identificar posibles observaciones que puedan influir en las estimaciones, una vez reconocidos estos datos, se procede a revisar si los datos son reales y plausibles, si hay errores se corrigen, de lo contrario se pueden mostrar las estimaciones cuando los datos están presenten y cuando no.

6. Referencias

Dalpiaz David, Applied Statistics with R.

Fox & Weisberg (2019) – An R Companion to Applied Regression

James, Witten, Hastie & Tibshirani (2013) – An Introduction to Statistical Learning.

Weisberg (2005) – Applied Linear Regression

Zhang & Wang (2017-2025) – Advanced Statistics using R