Regresión Lineal Simple con R: Fundamentos

Caso Práctico

Introducción

Las notas contenidas en este documento, al igual que algunas imágenes han sido tomadas en su mayoría de la serie de videos del canal de YouTube Brandon Foltz.

El objetivo es plasmar algunos fundamentos teóricos de la Regresión Lineal Simple haciendo uso del lenguaje R.

Los supuestos matemáticos de la regresión lineal simple no se presentan en el documento (linealidad, independencia, homocedasticidad, etc.).

Caso Práctico

Propinas por el servicio: supongamos que somos dueños de un pequeño restaurante y deseamos predecir el monto de las propinas que se generan.

El monto de la propina está relacionado con el valor total de la factura, es decir, una propina dejada por valor de una factura de 5 dólares será menor, que la esperada por una factura de 50 dólares.

Se decide tomar entonces durante una noche, los datos de las propinas recibidas de seis comidas aleatorias. Inicialmente sólo se toma la información del monto de las propinas en dólares, pero no el valor de la factura pagada.

Los datos se muestran en la siguiente tabla:

propinas <- data.frame(plato = c(1, 2, 3, 4, 5, 6),
                       propina = c(5, 17, 11, 8, 14, 5))

propinas

Pregunta: ¿cómo se podría predecir el valor de propinas futuras, sólo con esta información?.

Para dar respuesta a esa pregunta, la mejor herramienta que se tiene es calcular la media, ya que sólo disponemos de una sola variable (la propina).

El promedio de las propinas es:

\[\bar{y} = \frac{5+17+11+8+14+5}{6}=10\]

mean(propinas$propina)

[1] 10

Con una sola variable y sin ninguna información adicional, el mejor predictor para la siguiente medición es la media de la muestra.

Con una sola variable la mejor predicción posible de obtener es de 10 US$, para el valor de las propinas.

Por lo tanto la variabilidad está explicada por el valor de las propinas en si mismas.

Gráfico propinas

El gráfico muestra los valores de las propinas dadas por cada comida y la línea roja punteada señala la media de las propinas.

library(tidyverse)
propinas %>% 
  ggplot(mapping = aes(x = plato, y = propina))+
  geom_point(color = "cadetblue")+
  labs(x = "Plato",
       y = "Propina US$",
       title = "Valor Propina",
       caption = "La línea a trazos represeta el promedio.")+
  scale_x_continuous(n.breaks = 6)+
  scale_y_continuous(n.breaks = 5)+
  geom_hline(yintercept=10, linetype="dashed", color = "chocolate")+
  geom_text(aes(label=propina),hjust= -0.5, vjust=0.5)+
  theme_bw()

Residuales

¿Qué tan buena es esta predicción hecha por la línea que representa la media?.

Algunos de los valores están por encima y otros por debajo de esta línea.

Tomando la distancia medida desde la media (10 US$) a los valores observados, obtenemos los valores residuales y así verificar la bondad de ajuste del “modelo”.

\(Residuales = x_i - \bar{x}\)

luego,

\(Residual\ x_1 = 5-10 = -5\)

\(Residual\ x_2 = 17-10 = 7\)

\(Residual\ x_3 = 11-10 = 1\)

\(Residual\ x_4 = 8-10 = -2\)

\(Residual\ x_5 = 14-10 = 4\)

\(Residual\ x_6 = 5-10 = -5\)

Calculando los residuales con R:

residuales <- (propinas$propina - mean(propinas$propina))
residuales

[1] -5  7  1 -2  4 -5

Estos residuales son llamados también, el error, el cual indica cuánto se aleja el valor observado, de la linea de ajuste (media).

Si se suman los residuales que se encuentran por debajo de la media se obtiene un valor de -12.

Y si se suman los residuales que se encuentran por encima de la línea de ajuste el valor obtenido es 12.

Al sumar los residuales:

\[\sum(Residuales) = -12+12=0\]

sum(residuales)

[1] 0

La suma algebraica de los residuales siempre es cero.

Suma Cuadrada del Error

Se elevará el valor de los residuales al cuadrado, como lo muestra la siguiente tabla, con el objetivo de obtener valores positivos de los mismos, de la misma manera que en la desviación estándar.

\[Des.Est = \sqrt{\frac{\sum(x-\bar{x})^2}{n}}\]

data.frame(plato = propinas$plato,
           residuales = residuales,
           "resid^2" = (residuales^2))

La suma de estos residuales que estan elevados al cuadrado, se conoce como la Suma de Cuadrados del Error o la Suma de Cuadrados de los Residuales (SCE). En inglés inglés Sum of Squares Error (SSE).

Obteniendo la suma de cuadrados del error:

\[SCE = 25+49+1+4+16+25=120\]

sum(residuales^2)

[1] 120

Objetivo

El objetivo de la regresión lineal es crear un modelo lineal que minimice la suma de cuadrados del error.

Su propósito principal es estimar la función de regresión poblacional con base en la función de regresión muestral.

Dicho de otra manera, el objetivo es obtener un modelo (línea de ajuste) agregando una variable independiente que minimice el tamaño de los cuadrados tanto como sea posible.

Se obtiene una línea que se ajusta mejor a los datos tomados del problema y que miminiza los residuales.

Estadística Bivariada

El estudio del análisis de varianza ANOVA y de la correlación, se relacionan con la regresión lineal simple.

“Repaso de Álgebra”: Línea Recta

Modelo de Regresión Lineal

Para el modelo de regresión lineal se hace necesario la existencia de dos variables, a diferencia del caso anterior en el cual sólo se contaba con información de una sola variable (propina).

• Modelo matemático

\[y = mx+b\] El modelo de regresión lineal para la población de datos que se considere está dada por:

• Modelo estadístico

\[y = \beta_{0}+\beta_{1}x\]

\[\hat{y} = \hat{\beta_0}+\hat\beta_{1}X_i+\epsilon\]

En esencia es igual a la recta pendiente intercepto de arriba.

• \(\beta_{0} =\) intercepto con el eje \(y\) del parámetro de la población.

• \(\beta_{1}=\) pendiente del parámetro de la población.

• \(\epsilon=\) término para el error, variación que no se puede explicar de la variable \(y\).

La función o ecuación que describe la regresión lineal simple es:

\[E(y)= \beta_{0}+\beta_{1}x\] El valor esperado de \(y\) es la media, para un valor dado de \(x\).

El valor esperado de \(y\) es realmente la media de una distribución alrededor de los valores de \(y\).

A manera de ejemplo y tomando la base de datos cars de R, se presenta el modelo de regresión lineal para el conjunto de datos:

library(ggpmisc)
cars %>% 
  ggplot(mapping = aes(x = dist, y = speed))+
  geom_point()+
  geom_smooth(method = "lm", se = F)+
  labs(title = "Línea de Regresión", y = "Distancia", x = "Velocidad", subtitle = "Base de Datos cars")+
  stat_poly_eq(aes(label = paste(..eq.label.. ,
                                 ..rr.label.. ,
                                 sep = "~~~~")),
               formula = y ~ x, parse = TRUE,
               label.x.npc = 0.5,
               color = "blue")

Ecuación de regresión con estimaciones

Si se conocen los valores para los parámetros \(\beta_{0}\) y \(\beta_{1}\), es posible usar la ecuación de regresión lineal simple: \[E(y)=\beta_{0}+\beta_{1}x\]

Pero en realidad casi nunca se conocen los valores de los parámetros de la población, por lo tanto se deben estimar usando una muestra de los datos.

Cuando usamos una muestra de los datos la ecuación se escribe como:

\[\hat{y}= b_{0}+b_{1}x\]

donde \(\hat{y}\) es el estimador para el valor esperado de \(y\).

\(E(y)\) es la media del valor \(y\) para un valor dado de \(x\).

Percepción Geométrica

NOTA: en la literatura, en varios casos la diferencia entre las ecuaciones, se refleja en el uso de las letras b minúsculas, para remarcar que se trabaja con una muestra de la población.

Aplicado al ejemplo en el cual sólo se tiene información de una variable, la propina, el valor de la pendiente sería cero, ya que de hecho la variable \(x\) no existe, por lo tanto:

El valor de \(\hat{y}\) es \(10\) para todo valor de \(x\).

Siguiendo con el mismo ejemplo del restaurante si se agrega la información del total de la factura y el valor de la propina; en este caso se dice que el valor de la propina depende del monto total de la factura (US$).

Agregando una nueva columna a los datos:

propinas$factura <- c(34, 108, 64, 88, 99, 51)

propinas

Hipotéticamente se dice que el valor de la propina dependerá del valor total de la factura, una factura por un valor bajo dará como resultado una propina más baja comparada con el monto de una propina dejada por una factura de más alto valor.

Si la línea de tendencia resultante observada, reduce la suma de cuadrados del error, obtenida en el caso en el que se tenía sólo la variable propinas que usaba el modelo de la media de \(10 US\$\), se dice que el modelo de regresión es mejor cualitativamente comparado con el modelo de la media.

Se podrá calcular de manera cuantitativa cuánto es mejor un modelo en comparación con otro.

Método de Mínimos Cuadrados

Ahora que se tiene la información en pares, sobre el valor total de la factura y la propina dejada de esa factura; se desea conocer en qué grado el monto de la propina puede ser predicha por el valor total de la factura.

Se puede afirmar entonces que la variable propina es la variable dependiente y la variable independiente es la factura.

Este planteamiento es importante hacerlo para determinar el rol de las variables en nuestro modelo, no tiene sentido alguno afirmar que el valor total de la factura está dado por el monto de la propina.

Criterio de Mínimos Cuadrados

\[min\sum(y_{i}-\hat{y}_i)^2\]

• \(y_{i}=\) es el valor observado de la variable dependiente (propina).

• \(\hat{y}_i=\) valor estimado (predicho) de la variable dependiente (valor predicho de la propina). Basado en un modelo de regresión.

Se tendrán dos valores para cada valor de \(x\), el valor observado (real) y el valor predicho.

El objetivo es minimizar la suma de las diferencias cuadradas entre el valor observado para la variable dependiente (\(y_i\)) y el valor predicho o estimado de la variable dependiente (\(\hat{y}_i\)) que esta dado por la línea de regresión. Suma de los cuadrados del error.

Entonces la suma de los cuadrados residuales debe ser mucho menor que la obtenida sólo con la variable dependiente; \(\beta_1=0\), \(\hat{y}=10\) para cada valor de \(x\). La suma de los cuadrados debe ser mucho menor que \(120\) que fué la obtenida.

Gráfico de dispersión

Un gráfico de dispersión arroja información acerca de el comportamiento que siguen los datos, hacer uso de una escala acorde a los datos, facilita la interpretación del gráfico.

Infinitas rectas pueden ajustarse en menor o mayor medida a los datos, la idea es obtener aquella que minimice entonces, la suma de los cuadrados del error.

En comparación con los residuales obtenidos con la media de las propinas (línea horizonal), la suma de los cuadrados residuales de las rectas con pendiente diferente a cero, es menor.

propinas %>% 
  ggplot(mapping = aes(x = factura, y = propina))+
  geom_point(color = "blue")+
  geom_smooth(method = "lm", se = FALSE, color = "orange", linetype = "dashed")+
  labs(x = "Factura",
       y = "Propina",
       title = "Posibles Líneas de Regresión")+
  geom_abline(slope = 0.146,
              intercept = 1,
              linetype = "dashed",
              color = "forestgreen",
              size = 1)+
  geom_abline(slope = 0.146,
              intercept = -2,
              linetype = "dashed",
              color = "purple",
              size = 1)+
  geom_text(aes(label=propina),hjust=-0.5, vjust=0)+
  theme_bw()

Correlación

El coeficiente de correlación de Pearson es una medida lineal entre dos variables aleatorias cuantitativas. A diferencia de la covarianza, la correlación es independiente de la escala de medida.

Este coeficiente puede ser de dos tipos:

• Paramétrico: sujeto a distribución normal o gausiana.

• No paramétrico: no está sujeto a distribución normal o gausiana.

Útil para obtener el coeficiente de correlación de el modelo de regresión lineal.

\[r(x,y) = \frac{Cov_{x,y}}{\sigma_x \sigma_y} = \frac{\sum(x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum(x_i - \bar{x})^2}*\sqrt{\sum(y_i - \bar{y})^2}}\]

El coeficiente de correlación \(r\) permite verificar si la relación obtenida en los datos es lineal, positiva o negativa; también indica que no hay correlación lineal, en caso de ser cero. Entre mas cercano es a 1 o -1 es mas fuerte, entre mas cercano a 0 es débil hasta llegar hacerse nula.

• Correlación perfecta negativa = \(-1\) .
• No existe correlación = \(0\).
• Correlación perfecta positiva = \(1\).

Calculando la correlación de las variables con R:

cor(propinas$propina, propinas$factura)

[1] 0.865665

También se puede obtener la matriz de correlación de las variables del conjunto de datos:

cor(propinas)

             plato    propina   factura
plato    1.0000000 -0.1309307 0.1511229
propina -0.1309307  1.0000000 0.8656650
factura  0.1511229  0.8656650 1.0000000

Estadística Descriptiva de las Variables

Conocer las variables de manera gŕafica y numérica, es una buena práctica que arroja información importante para validar las condiciones que permitan la aplicación del modelo.

Obtener las medias para cada una de las variables:

• \(\bar{x}\) = \(74\) (valor promedio factura)

• \(\bar{y}\) = \(10\) (valor promedio propina)

mean(propinas$factura)

[1] 74

mean(propinas$propina)

[1] 10

El punto que se ubica en la coordenada \(\bar{x} = 74\) , \(\bar{y} = 10\), se conoce como centroide, y es el punto formado por las medias de las variables.

La línea de regresión que mejor se ajuste pasará por el centroide.

Resumen tabular de los valores numéricos necesarios para el modelo de regresión lineal:

propinas %>% 
  mutate(desvprop = propina-mean(propina),
         desvfact = factura-mean(factura),
         proddesv = desvfact*desvprop,
         cuadesvp = desvprop^2,
         cuadesvf = desvfact^2)

NA

Otros valores numéricos bastante útiles son:

• Mediana.

median(propinas$factura)

[1] 76

• Desviación Estándar.

sd(propinas$factura)

[1] 29.00345

• Varianza

var(propinas$factura)

[1] 841.2

• Covarianza

cov(propinas$propina, propinas$factura)

[1] 123

Obtener la Línea de Regresión

\[\hat{y}_i=b_{0}+b_{1}x_{i}\]

Pendiente: \[b_{1}=\frac{\sum(x_{i}-\bar{x})(y_{i}-\bar{y})}{\sum(x_{i}-\bar{x})^2} = 0.1462\]

donde,

• \(\bar{x}=\) media de la variable independiente.

• \(\bar{y}=\) media de la variable dependiente.

• \(x=\) valor de la variable independiente.

• \(y=\) valor de la variable dependiente.

Intercepto:

\[b_{0}= \bar{y}-b_{1}\bar{x} = -0.8188\]

Entonces, la fórmula que describe la línea de regresión lineal en términos de las variables es:

\[propina = -0.8188+(0.1462)factura\]

Para obtener la representación geométrica de la línea de regresión con R, hacemos uso de el método lm.

propinas %>% 
  ggplot(mapping = aes(x = factura, y = propina))+
  geom_point(color = "cadetblue")+
  geom_text(aes(label=propina),hjust=0, vjust=0)+
  geom_smooth(method = "lm", se = FALSE, color = "firebrick2")+
  geom_point(mapping = aes(x = 74, y = 10), color = "blue", size = 2)+
  annotate("text",
           label = "Centroide",
           x = 83,
           y = 10,
           color = "blue")+
  labs(x = "Factura",
       y = "Propina",
       title = "Línea de Regresión")+
   stat_poly_eq(aes(label = paste(..eq.label.. ,
                                 sep = "~~~~")),
               formula = y ~ x, parse = TRUE,
               label.x.npc = 0.5,
               color = "firebrick")+
  theme_bw()

Interpretación

\[\hat{y_i} = 0.1462x-0.8188\]

El valor de la pendiente para este caso, significa, que por cada dólar que aumenta el valor de la factura, esperaríamos que el monto para la propina aumente en \(0.1462\) \(US\$\), aproximadamente \(15\) centavos de dólar.

La interpretación del intercepto, puede o nó, tener un significado en el “mundo real”, dependerá entonces de el fenómeno que se esté estudiando. En este caso en particular, si el valor de la factura es cero, el modelo predice un valor de propina de \(-0.8188\) \(US\$\), lo cual no tiene sentido en la vida real ya se obtiene un valor negativo de propina.

Es importante anotar, que el modelo de regresión es único para el conjunto de datos que representa; adicionar o cambiar datos generará un cambio en el modelo.

Ajuste y Coeficiente de Determinación

¿Qué tan bien se ajusta la línea de regresión a los datos, comparado cuando se usa sólo el promedio de la variable dependiente?.

Recordemos que cuando calculamos la suma de los cuadrados del error para la variable dependiente se obtuvo: \(SSE=120\). La suma total de los cuadrados del error (SST = SSE) nunca será mayor a 120 para este conjunto de datos.

Tengamos presente también, que el objetivo es reducir el SSE a través de la línea que mejor se ajuste al cojunto de datos, ya que cuanto mejor se ajuste el modelo, menor será el valor del SSE.

Vamos a obtener la suma de cuadrados del error debido a la regresión (SSR), la cual es la diferencia entre SST y SSE.

Obtengamos con el modelo, los valores que predice para el monto de las propinas:

\[\hat{y_i} = 0.1462x_i-0.2188\]

donde,

\(x_i = factura_i\),

\(\hat{y} = predicción\ propina\)

Por ejemplo para el plato 3 se tiene una factura por valor de \(64US\$\). Por lo tanto el valor que predice el modelo es:

\[\hat{y_i} = 0.1462(64)-0.2188 = 8.5365\]

Ahora veamos en la tabla los valores predichos por el modelo:

mod1 <- lm(propina ~ factura, data = propinas)

propinas %>% 
  mutate(prop_pred = predict(mod1))

Error Predicciones

La diferencia que existe entre el valor observado y el valor predicho, es el error, es decir la distancia que hay entre el la predicción y el valor observado. También se les conoce como residuales.

propinas %>% 
  mutate(prop_pred = predict(mod1)) %>%
  ggplot(mapping = aes(y = prop_pred, x = factura))+
  geom_smooth(method = "lm", se = F)+
  geom_point(color = "darkred")+
  geom_text(aes(label=round(prop_pred, digits = 2)),hjust=0, vjust=1.5, color = "darkred")+
  geom_point(mapping = aes(x = factura, y = propina), color = "black")+
  labs(title = "Error Modelo de Regresión Lineal",
       x = "Factura",
       y = "Propina")+
  stat_poly_eq(aes(label = paste(..eq.label.. ,
                                 sep = "~~~~")),
               formula = y ~ x, parse = TRUE,
               label.x.npc = 0.5,
               color = "blue")+
  theme_bw()

• Tabla de error y cuadrado del error

propinas %>% 
  mutate(prop_pred = predict(mod1),
         error = propina-prop_pred,
         cuad_error = (error)^2)

La suma de los cuadrados del error SSE es 30.074893:

sum((propinas$propina - predict(mod1))^2)

[1] 30.07489

Podemos apreciar el comparativo de los valores de SSE cuando sólo se toma como predictor el promedio de la variable y cuando aplicamos regresión lineal simple; esto resume el objetivo de la regresión lineal, el cual es reducir el valor de la suma de los cuadrados del error.

El modelo redujo el valor del error, de 120 a 30.075, es decir que se redujo en 89.925 unidades.

\[SST = SSR + SSE\]

\[120 = 89.925 + 30.075\] Como se había mencionado anteriormente, el valor total máximo que se puede obtener es de 120, pero para el caso en el que se aplica la regresión lineal se obtiene un valor para \(SSE=30.075\) y para la suma de los cuadrados del error debido a la regresión, es de \(SSR=89.925\).

Coeficiente de Determinación

¿Qué tan bien se ajusta la línea de regresión a los datos?.

En este punto es donde la regresión comienza a tener similitud con ANOVA; la suma total de los cuadrados SST es dividida entre SSR y SSE y luego se mide la relación entre SSR y SST.

Si SSR es grande, le corresponde una mayor parte de SST y, por lo tanto, SSE es más pequeño en relación con SST.

El coeficiente de determinación \(r^2\) está dado por:

\[r^2 = \frac{SSR}{SST} = \frac{89.925}{120} = 0.749375\]

Interpretación

El valor permite medir si el modelo es, o no , estadísticamente siginificativo.

De la misma manera el valor se puede dar como porcentaje, en este caso \(74.93\%\).

Es posible concluir, que el \(74.93\%\) de la suma total de los cuadrados, puede ser explicado aplicando la ecuación de regresión para predecir el monto de la propina. El resto (\(25.07\%\)) es debido al error.

De esta información podemos decir que el ajuste que se obtiene es bueno.

propinas %>% 
  ggplot(mapping = aes(x = factura, y = propina))+
  geom_point(color = "orange")+
  geom_text(aes(label=propina),hjust=0, vjust=0)+
  geom_smooth(method = "lm", se = FALSE, color = "forestgreen")+
  labs(x = "Factura",
       y = "Propina",
       title = "Coeficiente de Determinación")+
   stat_poly_eq(aes(label = paste(..eq.label.. ,
                                  ..rr.label.. ,
                                 sep = "~~~~")),
               formula = y ~ x, parse = TRUE,
               label.x.npc = 0.5,
               color = "steelblue")+
  theme_bw()

Diferencias entre la suma de cuadrados

\[SSE = \Sigma(y-\hat{y_i})^2\]

\[SST = \Sigma(y_i - \bar{y})^2\]

\[SSR = \Sigma(\hat{y_i}- \bar{y})^2\]

Resumen del Modelo

La función lm() de R arroja un resumen del modelo de la regresión lineal:

mod1 <- lm(propina ~ factura, data = propinas)
resumen_modelo <- summary(mod1)
resumen_modelo

Call:
lm(formula = propina ~ factura, data = propinas)

Residuals:
      1       2       3       4       5       6 
 0.8488  2.0285  2.4622 -4.0471  0.3445 -1.6369 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) -0.82026    3.32297  -0.247   0.8172  
factura      0.14622    0.04228   3.458   0.0259 *
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.742 on 4 degrees of freedom
Multiple R-squared:  0.7494,    Adjusted R-squared:  0.6867 
F-statistic: 11.96 on 1 and 4 DF,  p-value: 0.02586

Es conveniente también acceder a los estadísticos del modelo de manera individual:

# predicciones del modelo
predict(mod1)

        1         2         3         4         5 
 4.151213 14.971469  8.537803 12.047076 13.655492 
        6 
 6.636947 

# coeficientes
coefficients(mod1)

(Intercept)     factura 
 -0.8202568   0.1462197 

# residuales
residuals(mod1)

         1          2          3          4          5 
 0.8487874  2.0285307  2.4621969 -4.0470756  0.3445078 
         6 
-1.6369472 

# error estandar del ajuste
summary(mod1)$sigma

[1] 2.742029

# cuadrado del error
summary(mod1)$r.squared  

[1] 0.7493759

Apliquemos ANOVA a nuestro modelo:

anova(mod1)

Analysis of Variance Table

Response: propina
          Df Sum Sq Mean Sq F value  Pr(>F)  
factura    1 89.925  89.925   11.96 0.02586 *
Residuals  4 30.075   7.519                  
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Hasta acá se tienen las herramientas básicas para aplicar el modelo de regresión lineal simple y algunas funciones de R para obtener algunos resultados descriptivos.