Apuntes de R

Parte 3: Introducción a regresión lineal simple

Author

David Torres Irribarra, Alén Amigo Quintana e Isidora Naranjo López

Published

April 4, 2024

Regresión lineal simple

Para el introducir el uso de regresión lineal simple, se presentrá el uso del modelo compacto y el modelo ampliado de regresión. Trabajaremos con la media de neuroticismo y la media de extraversión de la base de datos de los “Big five”.

# Leemos y guardamos la base de datos en el objeto big5
big5 <- read.csv("https://david-ti.github.io/introstats/data/big_five.csv")

# Itéms invertidos
big5$n2i <- (5 + 1) - big5$n2
big5$n4i <- (5 + 1) - big5$n4

# Creamos la variable a partir de los promedios de neuroticismo
big5$media_neuro <- rowMeans(big5[,c("n1","n2i","n3","n4i","n5","n6","n7","n8","n9","n10")]
                             , na.rm = TRUE)

# Creamos la variable a partir de los promedios de extroversión
big5$media_extro <- rowMeans(big5[,c("e1","e2","e3","e4","e5","e6","e7","e8","e9","e10")]
                             , na.rm = TRUE)

Modelo compacto

A continuación, generamos un resumen de un modelo lineal en el cual predecimos los datos de la variable neuroticismo en función del intercepto. Para calcular un modelo de regresión lineal haremos uso de la función lm.

# Visaulizamos el resumen
summary(lm(big5$media_neuro ~ 1))

Call:
lm(formula = big5$media_neuro ~ 1)

Residuals:
   Min     1Q Median     3Q    Max 
-2.039 -0.639  0.061  0.661  1.961 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.03900    0.03777   80.46   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8445 on 499 degrees of freedom

¿Cuales son los datos relevantes que debemos considerar?

Dato Valor Significado
Estimate 3.03900 Promedio
Residual standard error 0.8445 Desviación estándar (basada en la suma de cuadrados)

Modelo ampliado

Ahora, generamos un resumen de los resultados de un modelo lineal en el cual predecimos los datos de la variable neuroticismo en función de la variable extraversión.

# Visaulizamos el resumen
summary(lm(big5$media_neuro ~ 1+big5$media_extro))

Call:
lm(formula = big5$media_neuro ~ 1 + big5$media_extro)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.1004 -0.6268  0.0119  0.6242  1.9732 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)        3.4192     0.3367  10.156   <2e-16 ***
big5$media_extro  -0.1226     0.1079  -1.137    0.256    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8443 on 498 degrees of freedom
Multiple R-squared:  0.002587,  Adjusted R-squared:  0.0005841 
F-statistic: 1.292 on 1 and 498 DF,  p-value: 0.2563

También, podemos visualizar el resumen del siguiente modo:

# Visaulizamos el resumen
summary(lm(big5$media_neuro ~ big5$media_extro))

Call:
lm(formula = big5$media_neuro ~ big5$media_extro)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.1004 -0.6268  0.0119  0.6242  1.9732 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)        3.4192     0.3367  10.156   <2e-16 ***
big5$media_extro  -0.1226     0.1079  -1.137    0.256    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8443 on 498 degrees of freedom
Multiple R-squared:  0.002587,  Adjusted R-squared:  0.0005841 
F-statistic: 1.292 on 1 and 498 DF,  p-value: 0.2563

¿Cuales son los datos relevantes que debemos considerar?

Dato Valor Significado
Estimate / (Intercept) 3.4192 Coeficiente del intercepto
Estimate / big5$media_extro -0.1226 Coeficiente de la pendiente
Residual standard error 0.8443 Desviación estándar de los residuos
Multiple R-squared 0.002587 Coeficiente de determinación (Reducción de la variabilidad no explicada del modelo)
Std. Error / (Intercept) 0.3367 Error estándar de la estimación del intercepto
Std. Error / big5$media_extro 0.1079 Error estándar de la estimación de la pendiente

Interpretación del modelo de regresión

La regresión predice un nivel promedio de 3.4 puntos en la escala de neuroticismo para aquellas personas un nivel promedio de 0 puntos en la escala de extraversión. Recuerden que nos indica donde se ubica la recta

La pendiente tiene un valor de -0.12, lo cual indica una relación negativa entre los puntajes en la escala de extraversión y los puntajes en la escala de neuroticismo. Nos inidca que se espera un cambio en el promedio de la variable neuroticismo de -0.12 puntos antes un aumento de 1 unidad en promedio en la variable extraversión.

La desviación estándar de los residuos es de 0.84 puntos, indicándonos que en general, al hacer una predicción en base a la recta de regresión esperamos que los residuos pueden estar típicamente entre 1.68 puntos más arriba o más abajo que nuestra predicción de la media condicional.

Por último, la introducción del predictor (en este caso, la variable extraversión) ha reducido los residuos del modelo base (medidos usando la suma de cuadrados de los residuos) en un 0.2%.

Dispersograma con recta de regresión

Para analizar gráficamente la relación entre estas variables, podemos generar un diagrama de dispersión que permita visualizar las medias de neuroticismo y las medias de extraversión. Adicionalmente, es posible graficar la recta de regresión correspondiente con la función abline, de manera que nos muestre la tendencia general de los datos.

# Graficamos
with(big5,plot(media_extro,jitter(media_neuro), 
               main="Dispersiograma", 
               xlab="Extraversión", 
               ylab="Neuroticismo"))

# Agregamos la recta de regresión
abline(lm(big5$media_neuro ~ big5$media_extro))

Comparación de modelos

Al momento de querer comparar disitntos modelos de regresión estaremos trabajando con comparaciones de modelos anidados, es decir, que uno de los modelos sea un caso especial del otro. En este caso trabajaremos con el modelo compacto y el modelo ampliado presentados anteriormente.

Nos interesa es evaluar que tanto cambia el error (la suma de cuadrados de los residuos) entre estos dos modelos, sabiendo que ese cambio se debe a los parámetros adicionales en el modelo ampliados.

R2 como indicador de comparación de ajuste

Tal como vimos en clases, podemos comparar modelos a partir del cálculo de la reducción proporcional del error, es decir, a partir del cálculo de R2. Este valor se encuentra denominado como Multiple R-squared en el resumen que nos brinda R del cálculo de la regresión lineal simple.

¿Cuáles son los pasos para llegar este valor?

  1. Calculamos la suma la cuadrados totales

\[ SST \]

SST = sum(resid(lm(big5$media_neuro ~ 1))^2)

SST
[1] 355.9095
  1. Calculamos la suma la cuadrados del modelo ampliado

\[ SSE_A \]

SSEA = sum(resid(lm(big5$media_neuro ~ big5$media_extro))^2)

SSEA
[1] 354.9888
  1. Calculamos la reducción en la suma la cuadrados

\[ SSR = {SST - SSE_A} \]

SSR = SST - SSEA

SSR
[1] 0.9207307
  1. Calculamos la reducción proporcional del error

\[ R^2 = \dfrac{SST - SSE_A}{SST} \]

R_cuadrado = (SST - SSEA)/SST

R_cuadrado
[1] 0.00258698

Vemos que este equivale al Multiple R-squared en el resumen que nos brinda R.

En este caso vemos que el error cambia un 0.2% entre los dos modelos.

Relación del Coeficiente de Regresión estandarizado y el Coeficiente de Correlación

El coeficiente de regresión estandarizado y el coeficiente de correlación están relacionados, ya que ambos indican la fuerza y la dirección de la relación entre dos variables. Mientras que el coeficiente de correlación mide la asociación lineal entre las variables en su escala original, el coeficiente de regresión estandarizado proporciona esta medida después de estandarizar las variables, lo que facilita la comparación entre diferentes conjuntos de datos al tener una escala común (media de 0 y desviación estándar de 1).

### Forma 1

#Convertimos a puntaje Z los valores de ambas variables utilizadas en el análisis.
big5$Puntajes_Z_me <- scale(big5$media_extro)
big5$Puntajes_z_mn <- scale(big5$media_neuro)

#Hacemos la regresión para extroversión como variable dependiente y neuroticismo como variable independiente, con los puntajes estandarizados en vez de los originales.
Regresion1 = lm(big5$Puntajes_z_mn~big5$Puntajes_Z_me)
summary(Regresion1)

Call:
lm(formula = big5$Puntajes_z_mn ~ big5$Puntajes_Z_me)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.48700 -0.74216  0.01409  0.73905  2.33644 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)
(Intercept)        -1.964e-16  4.471e-02   0.000    1.000
big5$Puntajes_Z_me -5.086e-02  4.475e-02  -1.137    0.256

Residual standard error: 0.9997 on 498 degrees of freedom
Multiple R-squared:  0.002587,  Adjusted R-squared:  0.0005841 
F-statistic: 1.292 on 1 and 498 DF,  p-value: 0.2563
#Utilizamos el análisis de correlación para comprobar que el valor obtenido es igual al obtenido en el coeficiente de regresión del análisis anterior.
cor.test(big5$media_neuro,big5$media_extro, use=pairwise.complete.obs) 

    Pearson's product-moment correlation

data:  big5$media_neuro and big5$media_extro
t = -1.1365, df = 498, p-value = 0.2563
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.13793773  0.03699322
sample estimates:
        cor 
-0.05086236

Forma 2

#Hacemos la regresión para la media de neuroticismo como variable dependiente y la media de extroversión como variable independiente, con los puntajes estandarizados en vez de los originales.
Regresion2 = lm(scale(big5$media_neuro)~scale(big5$media_extro))
summary(Regresion2)

Call:
lm(formula = scale(big5$media_neuro) ~ scale(big5$media_extro))

Residuals:
     Min       1Q   Median       3Q      Max 
-2.48700 -0.74216  0.01409  0.73905  2.33644 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)
(Intercept)             -1.964e-16  4.471e-02   0.000    1.000
scale(big5$media_extro) -5.086e-02  4.475e-02  -1.137    0.256

Residual standard error: 0.9997 on 498 degrees of freedom
Multiple R-squared:  0.002587,  Adjusted R-squared:  0.0005841 
F-statistic: 1.292 on 1 and 498 DF,  p-value: 0.2563
#Utilizamos el análisis de correlación para comprobar que el valor obtenido es igual al obtenido en el coeficiente de regresión del análisis anterior.
cor.test(big5$media_neuro,big5$media_extro, use=pairwise.complete.obs) 

    Pearson's product-moment correlation

data:  big5$media_neuro and big5$media_extro
t = -1.1365, df = 498, p-value = 0.2563
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.13793773  0.03699322
sample estimates:
        cor 
-0.05086236

Forma 3

#Creamos el objeto que contiente los resultados del análisis de regresión lineal para el neuroticismo como variable dependiente y la extroversión como variable independiente.
Regresion3 = lm(big5$media_neuro~big5$media_extro)

#Utilizamos la función lm.beta sobre el objeto Regresion1 para obtener el coeficiente de regresión estandarizado. Para poder utilizar esta función debemos descargar y cargar la librería QuantPsyc.
#install.packages("QuantPsyc")
library(QuantPsyc)
Loading required package: boot
Loading required package: dplyr

Attaching package: 'dplyr'
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union
Loading required package: purrr
Loading required package: MASS

Attaching package: 'MASS'
The following object is masked from 'package:dplyr':

    select

Attaching package: 'QuantPsyc'
The following object is masked from 'package:base':

    norm
#Utilizamos la función lm.beta para obtener el coeficiente de regresión estandarizado.
lm.beta(Regresion3)
big5$media_extro 
     -0.05086236 
#Utilizamos el análisis de correlación para comprobar que el valor obtenido es igual al obtenido en el coeficiente de regresión estandarizado.
cor.test(big5$media_neuro,big5$media_extro, use=pairwise.complete.obs)

    Pearson's product-moment correlation

data:  big5$media_neuro and big5$media_extro
t = -1.1365, df = 498, p-value = 0.2563
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.13793773  0.03699322
sample estimates:
        cor 
-0.05086236

Como podemos ver en los outputs, el coeficiente de regresión estandarizado para la variable de extroversión es aproximadamente -0.05, mientras que el coeficiente de correlación de Pearson entre las variables de neuroticismo y extroversión es también aproximadamente -0.05. Ambos valores son coherentes y muestran una relación débil y negativa entre las dos variables.

Aclaración práctica

Es importante que sepamos que cuando realizamos un modelo de regresión en R y lo guardamos en un objeto no solo se alojan en este los datos que vemos cuando usamos la función “summary”, sino que además se guardan una serie de datos adicionales. Para entenderlo mejor vamos a realizar un modelo de regresión con la variable media neuro como variable de respuesta y la variable edad como variable independiente y lo guardaremos en el objeto “regresion1”. Para esto haremos uso de la función lm.

# Guardamos el modelo de regresión
regresion1 <- lm(big5$age ~ 1 + big5$media_neuro)

Ahora que hemos guardado el modelo de regresión en un objeto, revisemoslo.

# Guardamos el modelo de regresión
str(regresion1) 
List of 13
 $ coefficients : Named num [1:2] 32.12 -2.22
  ..- attr(*, "names")= chr [1:2] "(Intercept)" "big5$media_neuro"
 $ residuals    : Named num [1:497] -0.672 -5.561 -2.672 13.217 -11.672 ...
  ..- attr(*, "names")= chr [1:497] "1" "2" "3" "4" ...
 $ effects      : Named num [1:497] -565.46 -41.92 -2.62 13.42 -11.36 ...
  ..- attr(*, "names")= chr [1:497] "(Intercept)" "big5$media_neuro" "" "" ...
 $ rank         : int 2
 $ fitted.values: Named num [1:497] 27.7 26.6 25.7 26.8 27.7 ...
  ..- attr(*, "names")= chr [1:497] "1" "2" "3" "4" ...
 $ assign       : int [1:2] 0 1
 $ qr           :List of 5
  ..$ qr   : num [1:497, 1:2] -22.2935 0.0449 0.0449 0.0449 0.0449 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:497] "1" "2" "3" "4" ...
  .. .. ..$ : chr [1:2] "(Intercept)" "big5$media_neuro"
  .. ..- attr(*, "assign")= int [1:2] 0 1
  ..$ qraux: num [1:2] 1.04 1.03
  ..$ pivot: int [1:2] 1 2
  ..$ tol  : num 1e-07
  ..$ rank : int 2
  ..- attr(*, "class")= chr "qr"
 $ df.residual  : int 495
 $ na.action    : 'omit' Named int [1:3] 8 208 225
  ..- attr(*, "names")= chr [1:3] "8" "208" "225"
 $ xlevels      : Named list()
 $ call         : language lm(formula = big5$age ~ 1 + big5$media_neuro)
 $ terms        :Classes 'terms', 'formula'  language big5$age ~ 1 + big5$media_neuro
  .. ..- attr(*, "variables")= language list(big5$age, big5$media_neuro)
  .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:2] "big5$age" "big5$media_neuro"
  .. .. .. ..$ : chr "big5$media_neuro"
  .. ..- attr(*, "term.labels")= chr "big5$media_neuro"
  .. ..- attr(*, "order")= int 1
  .. ..- attr(*, "intercept")= int 1
  .. ..- attr(*, "response")= int 1
  .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..- attr(*, "predvars")= language list(big5$age, big5$media_neuro)
  .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
  .. .. ..- attr(*, "names")= chr [1:2] "big5$age" "big5$media_neuro"
 $ model        :'data.frame':  497 obs. of  2 variables:
  ..$ big5$age        : int [1:497] 27 21 23 40 16 22 29 23 16 19 ...
  ..$ big5$media_neuro: num [1:497] 2 2.5 2.9 2.4 2 2.2 2.9 4.4 2.2 2.9 ...
  ..- attr(*, "terms")=Classes 'terms', 'formula'  language big5$age ~ 1 + big5$media_neuro
  .. .. ..- attr(*, "variables")= language list(big5$age, big5$media_neuro)
  .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:2] "big5$age" "big5$media_neuro"
  .. .. .. .. ..$ : chr "big5$media_neuro"
  .. .. ..- attr(*, "term.labels")= chr "big5$media_neuro"
  .. .. ..- attr(*, "order")= int 1
  .. .. ..- attr(*, "intercept")= int 1
  .. .. ..- attr(*, "response")= int 1
  .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. .. ..- attr(*, "predvars")= language list(big5$age, big5$media_neuro)
  .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
  .. .. .. ..- attr(*, "names")= chr [1:2] "big5$age" "big5$media_neuro"
  ..- attr(*, "na.action")= 'omit' Named int [1:3] 8 208 225
  .. ..- attr(*, "names")= chr [1:3] "8" "208" "225"
 - attr(*, "class")= chr "lm"

Podemos observar que el objeto regresion1 contiene una serie de datos. Algunas de las partes más importantes son:

  • coefficients: Un vector con los coeficientes estimados (intercepto y pendientes).
  • residuals: Un vector con los residuos (diferencias entre valores observados y predichos).
  • fitted.values: Un vector con los valores predichos por el modelo.
  • rank: El rango del modelo, que indica cuántos coeficientes fueron estimados.
  • df.residual: Los grados de libertad residuales.
  • call: La llamada original a lm(), útil para reproducir el modelo.
  • formula: La fórmula usada en la regresión.
  • model: Los datos originales usados en la regresión.

Para efectos del curso es importante que podamos observar como se distribuyen los residuos. Para ello podemos realizar un histrograma de la siguiente manera.

# Histograma de residuos de regresion1
hist(regresion1$residuals)