Análisis de correlacion
Importar paquetes y definir folder de trabajo
Importar datos
grasas <- read.table("http://verso.mat.uam.es/~joser.berrendero/datos/EdadPesoGrasas.txt", header= TRUE)
names(grasas)## [1] "peso" "edad" "grasas"
Cuantificar el grado de relación lineal (coef. de correlación)
## peso edad grasas
## peso 1.0000000 0.2400133 0.2652935
## edad 0.2400133 1.0000000 0.8373534
## grasas 0.2652935 0.8373534 1.0000000Estimación 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:
##
## 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.2Los 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.622x \]
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):
Calcula y representa gráficamente la recta de regresión, junto con la correspondiente nube de puntos.
¿Cuánto vale el coeficiente de correlación al cuadrado en este caso? 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.02996.
Cálculo de predicciones
Supongamos que queremos utilizar la recta de mínimos cuadrados para predecir la cantidad de grasas para individuos de pesos 50, 51, 53 … 80. Basta crear un fichero de datos que contenga las nuevas variables regresoras y usar el comando predict:
## 1 2 3 4 5 6 7 8
## 280.4146 282.0370 283.6593 285.2817 286.9040 288.5264 290.1487 291.7710
## 9 10 11 12 13 14 15 16
## 293.3934 295.0157 296.6381 298.2604 299.8828 301.5051 303.1274 304.7498
## 17 18 19 20 21 22 23 24
## 306.3721 307.9945 309.6168 311.2391 312.8615 314.4838 316.1062 317.7285
## 25 26 27 28 29 30 31
## 319.3509 320.9732 322.5955 324.2179 325.8402 327.4626 329.0849Inferencia 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, \]
en el cual los errores aleatorios \(\epsilon_i\) son independientes normal de media 0 y varianza σ2.
Bajo este modelo,
¿Cuánto valen los estimadores de todos los parámetros del modelo? 1. 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 85.818 y 1.229 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.0294 y H0:β1=0.2000. Los correspondientes p-valores aparecen en la columna Pr(>|t|).
El estimador de la desviación típica de los errores σ aparece como Residual standard error y su valor en el ejemplo es 76.65.
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):
Contrasta la hipótesis de que la pendiente de la recta es cero a nivel 0.05.
## 2.5 % 97.5 %
## (Intercept) 21.7696109 376.825392
## peso -0.9209324 4.165618Calula un intervalo de confianza para la pendiente de la recta de nivel 90%.
## 5 % 95 %
## (Intercept) 52.2166142 346.378389
## peso -0.4847468 3.729432- 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 pesos que va de 30 a 90 kg (los de predicción en rojo):
Calcula y representa los intervalos de confianza al 95% de la cantidad de grasas media para los individuos entre 30 y 90 kg.
nuevas.peso <- data.frame(peso=seq(30,90))
#Gráfica de disperión y recta
plot(grasas$peso, grasas$grasas, xlab="Peso", ylab="Grasas")
abline(regresion)
#Intervalos de confianza de la respuesta media
#ic es una matriz de tres columnas: la primera es la predicción, las otras son los extremos del intervalo
ic <- predict(regresion, nuevas.peso, interval="confidence" )
lines(nuevas.peso$peso, ic[, 2], lty =2, col="red" )
lines(nuevas.peso$peso, ic[, 3], lty =2, col="red" )
#Intervalo de predicción
ic <- predict(regresion, nuevas.peso, interval="prediction" )
lines(nuevas.peso$peso, ic[, 2], lty =2, col="red" )
lines(nuevas.peso$peso, ic[, 3], lty =2, col="red" )
Análisis ANOVA
- La tabla de análisis de la varianza de los errores se obtiene con el comando anova:
## Analysis of Variance Table
##
## Response: grasas
## Df Sum Sq Mean Sq F value Pr(>F)
## peso 1 10232 10231.7 1.7413 0.2
## Residuals 23 135145 5875.9Diagnóstico del modelo
Los valores ajustados y^i y los residuos ei=y^i−yi 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:
residuos <- rstandard(regresion)
valores.ajustados <- fitted(regresion)
plot(valores.ajustados, residuos)
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.