Los datos

Trabajaremos con datos que corresponden a tres variables medidas en 20 personas que se encuentran en un programa de salud: edad, peso y cantidad de grasa en la sangre. Los datos son:

edad  <- c(46,20,52,30,57,25,28,36,57,44,24,31,52,23,60,48,34,51,40,34)
peso  <- c(84,73,65,70,76,69,63,72,79,75,27,89,65,57,59,69,60,79,75,82)
grasa <- c(354,190,405,263,451,302,288,325,402,365,209,290,346,254,395,434,220,374,308,220)
datos <-data.frame(edad,peso,grasa)

El código necesario para ajustar un modelo de regresión en R es bastante sencillo, pero es muy importante entender e interpretar correctamente todas las salidas.

Antes de realizar cualquier proceso de regresión, debemos conocer la relación existente entre cada par de variables, para poder concuir si vale la pena o no calcular un modelo de regresión.

pairs(datos)

Para cuantificar el grado de asociación lineal, calculamos el coeficiente de correlación lienal, el cual, cuando se tiene más de dos variables, se puede calcular como la matriz de coeficientes de correlación.

cor(datos)
##            edad      peso     grasa
## edad  1.0000000 0.2705334 0.8741798
## peso  0.2705334 1.0000000 0.3009460
## grasa 0.8741798 0.3009460 1.0000000
Graficamente y con la cuantificación hecha con el coeficiente de correlación, notamos una fuerte relación entre las variables grasas-edad \((r\approx 0.88)\), pero una muy baja relación entre peso-edad \((r\approx0.27)\) y peso-grasa \((r\approx0.30)\). Entonces, la única relación a la que vale la pena analizar un modelo de regresión es entre las variables grasas y edad.

Modelo de regresión lineal

La función lm() (por sus siglas en inglés, 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. Dos salidas importantes se pueden observar en la impresión del modelo y en la función summary() .

reg <-lm(grasa~edad,data=datos)
reg
## 
## Call:
## lm(formula = grasa ~ edad, data = datos)
## 
## Coefficients:
## (Intercept)         edad  
##     107.633        5.356

En esta primer salida, podemos observar dos cosas importantes Call: el modelo aplicado y Coefficients: la estimación de los coeficientes \(\beta_0=107.633\) y \(\beta_1=5.356\).

summary(reg)
## 
## Call:
## lm(formula = grasa ~ edad, data = datos)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -69.75 -25.37  -2.68  23.51  69.25 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 107.6328    29.1057   3.698  0.00165 ** 
## edad          5.3565     0.7013   7.638 4.71e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 38.94 on 18 degrees of freedom
## Multiple R-squared:  0.7642, Adjusted R-squared:  0.7511 
## F-statistic: 58.33 on 1 and 18 DF,  p-value: 4.712e-07

En esta segunda salida, podemos observar una infomración muy completa de todo el modelo. Veamos la interpretación de cada una de ellas.

  1. Call: como ya mencionamos, es el modelo aplicado.

  2. Residuals: un resumen de la distribución de los residuales.

  3. Coefficients: en esta salida, además de obtener la estimación de los coeficientes \(\beta_0\) y \(\beta_1\), también aparece en cada fila, los valores necesarios para realizar la prueba de hipótesis de cada coeficiente, donde se contrasta si el coeficiente es igual a cero contra que es diferente a cero. Para estas pruebas se calcula la desviación estándar (Std. Error) , el valor \(t\) calculado (t value) y el \(p\)-valor (Pr(>|t|)).

Para nuestro ejemplo, podemos observar que ambos \(p\)-valor son muy pequeños o cercanos a cero, por lo tanto, rechazariamos en ambas la hipótesis nula, es decir, los dos coeficientes son significativamente distintos de cero.

  1. La última salida esta compuesta por varios valores muy importantes:

  1. Residual standard error: la desviación estándar de los errores, con sus respectivos grados de libertad, que para nuestro ejemplo son \(n-2=18\).

  2. Multiple R-squared: coeficiente de determinación. Para nuestro ejemplo, obtenemos un \(R^2=0.7642\), es decir, tenemos una buena relación entre las variables de interés. También aparece Asjusted R-squared, el cual es relevante para regresión lineal múltiple y por ello para este tema no lo analizamos.

Por todo lo anterior, nuestro modelo de regresión será: \[\hat{y}_i = 107.633+ 5.356x_i \]

Podemos visualizar ahora los datos con la recta de regresión ajustada.

plot(edad,grasa,xlab="Edad",ylab="Grasa en la sangre")
abline(reg, col="blue")

Cálculo de predicciones

Con esta recta de regresión estimada, podríamos por ejemplo predecir la cantidad de grasa en la sangre para cualquier edad. Por ejemplo, podemos predecir para las edades: \(21,22,\dots,60\). Esto lo podemos hacer de forma muy sencilla con la función predict() .

pred_edades <- data.frame(edad=seq(20,60))
predict(reg,pred_edades)
##        1        2        3        4        5        6        7        8 
## 214.7627 220.1192 225.4757 230.8322 236.1887 241.5452 246.9017 252.2582 
##        9       10       11       12       13       14       15       16 
## 257.6147 262.9712 268.3277 273.6842 279.0406 284.3971 289.7536 295.1101 
##       17       18       19       20       21       22       23       24 
## 300.4666 305.8231 311.1796 316.5361 321.8926 327.2491 332.6056 337.9621 
##       25       26       27       28       29       30       31       32 
## 343.3186 348.6751 354.0316 359.3881 364.7446 370.1010 375.4575 380.8140 
##       33       34       35       36       37       38       39       40 
## 386.1705 391.5270 396.8835 402.2400 407.5965 412.9530 418.3095 423.6660 
##       41 
## 429.0225

Entonces, vemos por ejemplo que un invididuo de \(30\) años tiene una predicción de cantidad de grasa de \(268.3277\).

También se podría predecir la cantidad de grasa en la sangre para una persona con estas características pero con edad diferente al rango que hemos trabajado. Por ejemplo si queremos predecir para las edades de \(10\) y \(80\) años. Cabe resaltar que si nos alejamos demasiado al rango de edades donde estamos trabajando, la predicción puede tener mayor variabilidad.

pred_edad <- data.frame(edad=c(10,80))
predict(reg,pred_edad)
##        1        2 
## 161.1978 536.1524

Entonces, encontramos que un individuo de \(10\) años tendrá \(161.1978\) y un individuo de \(80\) años tendrá \(536.1524\) de catidad de grasa en la sangre.

Inferencia en el modelo de regresión simple

Teorícamente los datos proceden de un modelo de regresión simple de la forma: \[y_i = \beta_0 + \beta_1x_i +\varepsilon_i, \quad i=1,\dots,n.\] donde los errores aleatorios \(\varepsilon_i\) son independientes con distribución normal de media cero y varianza \(\sigma^2\) constante.

Recordemos que podemos encontrar el valor de los residuales del modelo así:

error <-reg$residuals

El valor de la media cero en los residuales, se consigie por la misma construcción, entonces no es necesaria verificarla. Al calcularla, nos debe dar un valor muy cercano a cero. Esto es:

mean(error)
## [1] -1.508516e-15

La estimación de la desviación estándar de los residuales, como ya se habia mencionado aparece Residual standard error y su valor para este ejemplo es 38.94.

Podemos calcular adicionalmente, intervalos de confianza. Estos se obtienen con la función confint() . Por ejemplo:

  • intervalos de confianza para los coeficientes con un nivel de \(98\%\):
confint(reg,level=0.98)
##                   1 %       99 %
## (Intercept) 33.343914 181.921753
## edad         3.566426   7.146563

Esto es, con un \(98\%\) de confianza los coeficientes se encuentran entre estos valores:

\[(33.34 < \beta_0 < 181.92) \quad \text{y} \quad (3.57 < \beta_1 < 7.15)\] - Intervalos de confianza para la respuesta media y los intervalos de predicción. Estos los podemos calcular con la función predic() .

# Intervalos para la respuesta media 
ic1 <- predict(reg,pred_edades,interval='confidence',level=0.98)
# Intervalos de prediccion
ic2 <- predict(reg,pred_edades,interval='prediction',level=0.98)

# Gráfico del modelo de regresión
plot(edad, grasa, xlab='Edad', ylab='Grasa en la sangre')
abline(reg,col="blue")

# Gráfico de los intervalos para la respuesta media 
lines(pred_edades$edad,ic1[,2],lty=2,col="purple")
lines(pred_edades$edad,ic1[,3],lty=2,col="purple")
# Gráfico de los intervalos de prediccion
lines(pred_edades$edad,ic2[,2],lty = 2,col = 'red')
lines(pred_edades$edad,ic2[,3],lty = 2,col = 'red')

Supuestos del modelo

Recordemos que los valores ajustados \(\hat{y}_i\) y los residuos \(\varepsilon_i=\hat{y}_i−y_i\) se pueden obtener del modelo de regresión calculado con los comandos reg$fitted.values y reg$residuals respectivamente. Revisamos los supuestos del modelo.

1. Supuesto de normalidad

residuales <- reg$residuals
#Gráfico
qqnorm(residuales,xlab="Cuantiles teóricos",ylab="Cuantiles muestrales")
qqline(residuales,col="blue")

# Pruebas 
shapiro.test(residuales)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales
## W = 0.97369, p-value = 0.83
library(nortest)
ad.test(residuales)
## 
##  Anderson-Darling normality test
## 
## data:  residuales
## A = 0.20161, p-value = 0.8604

Podemos concluir que tanto gráficamente como por medio de las pruebas estadísticas, los residuales el supuesto de normalidad.

2. Supuesto de varianza constante

aj <- reg$fitted.values
#Gráfico
plot(aj,residuales,xlab="Valores ajustados",ylab="Residuales",col="blue")
abline(0,0,col="blue")
abline(-60,0,col="red",lty=2)
abline(60,0,col="red",lty=2)
# Prueba
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.0.4
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.0.4
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

bptest(grasa~edad)
## 
##  studentized Breusch-Pagan test
## 
## data:  grasa ~ edad
## BP = 0.16464, df = 1, p-value = 0.6849

Podemos concluir tanto gráficamente como por medio de la prueba estadística que los residuales cumplen el supuesto de varianza constante.

3. Supuesto de independencia

plot(edad,residuales,xlab="Edad",ylab="Residuales",col="blue")
abline(0,0,col="blue")

dwtest(grasa~edad)
## 
##  Durbin-Watson test
## 
## data:  grasa ~ edad
## DW = 2.2209, p-value = 0.7253
## alternative hypothesis: true autocorrelation is greater than 0

Nuevamente, podemos concluir tanto gráficamente como por medio de la prueba estadística que los residuales cumplen el supuesto de independencia.

Nota importante: al cumplirse todos los supuestos del modelo, toda conclusión que se tenga del modelo ajustado tiene validez.