Análisis de regresión lineal

setwd("~/Documents/pye1pm")

library(pacman)
p_load("DT","prettydoc","xfun","readr")

Tenemos 3 variables en nuestro conjunto de datos:

Peso en Kg. de la persona
Edades de las personas
Contenido total de Grasas en sangre

grasas <- read.table("http://verso.mat.uam.es/~joser.berrendero/datos/EdadPesoGrasas.txt", header = TRUE)
names(grasas)

## [1] "peso"   "edad"   "grasas"

Con el fin de conocer las relaciones existentes entre cada par de variables podemos representar una matriz de diagramas de dispersión. Al parecer existe una relación lineal bastante clara entre la edad y las grasas, pero no entre los otros dos pares de variables. Por otra parte el fichero contiene un dato atípico.

pairs(grasas)

## Cálculo y representación de la recta de mínimos cuadrados

El comando básico es lm (linear models).

El primer argumento de este comando es una fórmula y ~ x en la que se especifica cuál es la variable respuesta o dependiente (y) y cuál es la variable regresora o independiente (x).

El segundo argumento, llamado data especifica cuál es el fichero en el que se encuentran las variables. El resultado lo guardamos en un objeto llamado regresion. Este objeto es una lista que contiene toda la información relevante sobre el análisis. Mediante el comando summary obtenemos un resumen de los principales resultados:

regresion <- lm(grasas ~ peso, data = grasas )
summary(regresion)

## 
## Call:
## lm(formula = grasas ~ peso, data = grasas)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -127.729  -53.686   -9.239   46.537  128.404 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  199.298     85.818   2.322   0.0294 *
## peso           1.622      1.229   1.320   0.2000  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 76.65 on 23 degrees of freedom
## Multiple R-squared:  0.07038,    Adjusted R-squared:  0.02996 
## F-statistic: 1.741 on 1 and 23 DF,  p-value: 0.2

Inferencia en el modelo de regresión simple

Suponemos ahora que los datos proceden de un modelo de regresión simple de la forma

\[ y_i = \beta_0 + \beta_1 x_i + \epsilon_i, \ \ \ \ i=1,\ldots,n, \] donde los errores aleatorios ϵi son independientes con distribución normal de media 0 y varianza σ2

Bajo este modelo:

Los errores típicos de los estimadores de los parámetros β0 y β1 se encuentran en la columna Std Error de la salida anterior. En el ejemplo, sus valores son 29.638 y 0.724 respectivamente.
La columna t value contiene el estadístico t, es decir, cociente entre cada estimador y su error típico. Estos cocientes son la base para llevar a cabo los contrastes H0:β0=0 y H0:β1=0 . Los correspondientes p-valores aparecen en la columna Pr(>|t|). En este caso son muy pequeños por lo que se rechazan ambas hipótesis para los niveles de significación habituales.
El estimador de la desviación típica de los errores σ aparece como Residual standard error y su valor en el ejemplo es 43.5
Los intervalos de confianza para los parámetros se obtienen con el comando confint. El parámetro level permite elegir el nivel de confianza (por defecto es 0.95):

confint(regresion)

##                  2.5 %     97.5 %
## (Intercept) 21.7696109 376.825392
## peso        -0.9209324   4.165618

confint(regresion, level=0.90)

##                    5 %       95 %
## (Intercept) 52.2166142 346.378389
## peso        -0.4847468   3.729432

Indice de correlación

de tipo pearson

cor(grasas)

##             peso      edad    grasas
## peso   1.0000000 0.2400133 0.2652935
## edad   0.2400133 1.0000000 0.8373534
## grasas 0.2652935 0.8373534 1.0000000

Los parámetros de la ecuación de la recta de mínimos cuadrados que relaciona la cantidad de grasas en la sangre en función del peso vienen dados por la columna ´Estimate´ de la tabla ´Coefficients´ de la salida anterior.

Por lo tanto, en este ejemplo la ecuación de la recta de mínimos cuadrados es:

\[ y = 199.298 + 1.622 x \]

y = 199.298 + 1.622*(20)
y

## [1] 231.738

Gráfica

Los siguientes comandos representan la nube de puntos (comando plot) y añaden la representación gráfica de la recta de mínimos cuadrados (comando abline aplicado al objeto generado por lm):

plot(grasas$peso, grasas$grasas, xlab="Peso en kilos de las personas", ylab="Contenido de grasa en sangre")
abline(regresion)

El coeficiente de determinación (es decir, el coeficiente de correlación al cuadrado) mide la bondad del ajuste de la recta a los datos. A partir de la salida anterior, vemos que su valor en este caso es Multiple R-squared: 0.07038

Estimación de predicciones

Supongamos que queremos utilizar la recta de mínimos cuadrados para predecir la cantidad de grasas para individuos de edades 31,31,32,…,50 . Basta crear un fichero de datos que contenga las nuevas variables regresoras y usar el comando predict:

nuevos.pesos <- data.frame(peso = seq(20,90))
predict(regresion, nuevos.pesos)

##        1        2        3        4        5        6        7        8 
## 231.7444 233.3667 234.9890 236.6114 238.2337 239.8561 241.4784 243.1008 
##        9       10       11       12       13       14       15       16 
## 244.7231 246.3454 247.9678 249.5901 251.2125 252.8348 254.4572 256.0795 
##       17       18       19       20       21       22       23       24 
## 257.7018 259.3242 260.9465 262.5689 264.1912 265.8136 267.4359 269.0582 
##       25       26       27       28       29       30       31       32 
## 270.6806 272.3029 273.9253 275.5476 277.1700 278.7923 280.4146 282.0370 
##       33       34       35       36       37       38       39       40 
## 283.6593 285.2817 286.9040 288.5264 290.1487 291.7710 293.3934 295.0157 
##       41       42       43       44       45       46       47       48 
## 296.6381 298.2604 299.8828 301.5051 303.1274 304.7498 306.3721 307.9945 
##       49       50       51       52       53       54       55       56 
## 309.6168 311.2391 312.8615 314.4838 316.1062 317.7285 319.3509 320.9732 
##       57       58       59       60       61       62       63       64 
## 322.5955 324.2179 325.8402 327.4626 329.0849 330.7073 332.3296 333.9519 
##       65       66       67       68       69       70       71 
## 335.5743 337.1966 338.8190 340.4413 342.0637 343.6860 345.3083

Por ejemplo, para un individuo de 30 años, predecimos una cantidad de grasas de 262.2

Los intervalos de confianza para la respuesta media y los intervalos de predicción para la respuesta se pueden obtener usando el comando predict. Por ejemplo, el siguiente código calcula y representa los dos tipos de intervalos para el rango de edades que va de 20 a 60 años (los de predicción en rojo):

#Gráfico de dispersión y recta 
plot(grasas$peso, grasas$grasas, xlab="Peso en kilos de las personas", ylab="Contenido de grasa en sangre")
abline(regresion)

#Intervalos de confianza de la respuesta media:
#ic es una matriz con tres columnas: la primera es la #prediccion, las otras dos son los extremos del intervalo

ic <- predict(regresion, nuevos.pesos, interval = "confidence")
lines(nuevos.pesos$peso, ic[,2], lty=2)
lines(nuevos.pesos$peso, ic[,3], lty=2)

Descarga este código

xfun::embed_file("A5U1.Rmd")

Download A5U1.Rmd

Descarga los datos

xfun::embed_file("grasas.csv")

Download grasas.csv