Regresión Lineal simple con R

Los datos

Los datos del fichero EdadPesoGrasas.txt corresponden a tres variables medidas en 25 individuos: edad, peso y cantidad de grasas en sangre. Para leer el fichero de datos y saber los nombres de las variables:

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.

Para cuantificar el grado de relación lineal, calculamos la matriz de coeficientes de correlación:

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

Cálculo y representacion 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 ~ edad, data = grasas)
summary(regresion)

Call:
lm(formula = grasas ~ edad, data = grasas)

Residuals:
    Min      1Q  Median      3Q     Max 
-63.478 -26.816  -3.854  28.315  90.881 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 102.5751    29.6376   3.461  0.00212 ** 
edad          5.3207     0.7243   7.346 1.79e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 43.46 on 23 degrees of freedom
Multiple R-squared:  0.7012,    Adjusted R-squared:  0.6882 
F-statistic: 53.96 on 1 and 23 DF,  p-value: 1.794e-07

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=102.575+5.321x\)

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):

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.701.

Cálculo 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:

nuevas.edades <- data.frame(edad = seq(30, 50))
predict(regresion, nuevas.edades)

       1        2        3        4        5        6        7        8        9 
262.1954 267.5161 272.8368 278.1575 283.4781 288.7988 294.1195 299.4402 304.7608 
      10       11       12       13       14       15       16       17       18 
310.0815 315.4022 320.7229 326.0435 331.3642 336.6849 342.0056 347.3263 352.6469 
      19       20       21 
357.9676 363.2883 368.6090 

Por ejemplo, para un individuo de 30 años, predecimos una cantidad de grasas de 262.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:

\(yi=β0+β1xi+ϵi, i=1,…,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 \(H_{0}:β_{0}=0\) y \(H_0:β_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) 41.265155 163.885130
edad         3.822367   6.818986

confint(regresion, level = 0.90)

                  5 %       95 %
(Intercept) 51.780153 153.370132
edad         4.079335   6.562018

• 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):

nuevas.edades <- data.frame(edad = seq(20, 60))
# Grafico de dispersion y recta
plot(grasas$edad, grasas$grasas, xlab='Edad', ylab='Grasas')
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, nuevas.edades, interval = 'confidence')
lines(nuevas.edades$edad, ic[, 2], lty = 2)
lines(nuevas.edades$edad, ic[, 3], lty = 2)

# Intervalos de prediccion
ic <- predict(regresion, nuevas.edades, interval = 'prediction')
lines(nuevas.edades$edad, ic[, 2], lty = 2, col = 'red')
lines(nuevas.edades$edad, ic[, 3], lty = 2, col = 'red')

• La tabla de análisis de la varianza de los errores se obtiene con el comando anova:

anova(regresion)

Analysis of Variance Table

Response: grasas
          Df Sum Sq Mean Sq F value    Pr(>F)    
edad       1 101933  101933  53.964 1.794e-07 ***
Residuals 23  43444    1889                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Diagnóstico del modelo

Los valores ajustados \(\hat{y}_i\) y los residuos \(e_i=\hat{y}_i−y_i\) se pueden obtener con los comandos fitted y residuals respectivamente. Los residuos estandarizados se obtienen con rstandard. Por ejemplo, el siguiente código obtiene una representación de los residuos estandarizados frente a los valores ajustados, que resulta útil al llevar a cabo el diagnóstico del modelo:

No se observa ningún patrón especial, por lo que tanto la homocedasticidad como la linealidad resultan hipótesis razonables.

La hipótesis de normalidad se suele comprobar mediante un QQ plot de los residuos. El siguiente código sirve para obtenerlo:

Dado que los puntos están bastante alineados, la normalidad también parece aceptable.

# Variable regresora (dieño fijo) y parámetros
x = seq(1,10)
beta0 <- 0
beta1 <- 1
sigma <- 0.3

# Genera la variable respuesta
y <- beta0 + beta1*x + rnorm(length(x), sd=sigma)

# Ajusta el modelo
reg <- lm(y~x)

# Extrae el valor de la pendiente estimada
coefficients(reg)[2]

       x 
1.031375 

# Resume el ajuste
summary(reg)

Call:
lm(formula = y ~ x)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.54850 -0.30062 -0.00401  0.17426  0.70650 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.06301    0.27941  -0.226    0.827    
x            1.03138    0.04503  22.904  1.4e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.409 on 8 degrees of freedom
Multiple R-squared:  0.985, Adjusted R-squared:  0.9831 
F-statistic: 524.6 on 1 and 8 DF,  p-value: 1.4e-08