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 ~ 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-07Los 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.5751 + 5.3207 x \]
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:
## 1 2 3 4 5 6 7 8
## 262.1954 267.5161 272.8368 278.1575 283.4781 288.7988 294.1195 299.4402
## 9 10 11 12 13 14 15 16
## 304.7608 310.0815 315.4022 320.7229 326.0435 331.3642 336.6849 342.0056
## 17 18 19 20 21
## 347.3263 352.6469 357.9676 363.2883 368.6090Inferencia 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,
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):
## 2.5 % 97.5 %
## (Intercept) 41.265155 163.885130
## edad 3.822367 6.818986## 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))
#Gráfica de disperión y recta
plot(grasas$edad, grasas$grasas, xlab="Edad", 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.edades, interval="confidence" )
lines(nuevas.edades$edad, ic[, 2], lty =2, col="red" )
lines(nuevas.edades$edad, ic[, 3], lty =2, col="red" )
#Intervalo de predicción
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" )
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)
## edad 1 101933 101933 53.964 1.794e-07 ***
## Residuals 23 43444 1889
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Diagnó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.
Métodos para realizar análisis de correlación
- Existen diferentes métodos para realizar un análisis de correlación
- Correlación de Pearson(r):
Mide una dependencia lineal entre dos variables (x e y). También se conoce como prueba de correlación paramétrica porque depende de la distribución de los datos. Se puede usar solo cuando x e y son de distribución normal. La gráfica de \(y = f (x)\) se denomina curva de regresión lineal.
- Kendall \(\tau\) y Spearman \(\rho\): Son coeficientes de correlación basados en rangos (no paramétricos)
El Método que por lo regular se usa más es el de Pearson
Funciones en R para análisis de correlación
cor () calcula el coeficiente de correlación prueba
cor.test () para asociación / correlación entre muestras emparejadas. Devuelve tanto el coeficiente de correlación como el nivel de significancia (o valor p) de la correlación.
- Los formatos simples para implementar esto en R son:
#cor(x, y, method = c("pearson", "kendall", "spearman"))
#cor.test(x, y, method=c("pearson", "kendall", "spearman"))- Consideraciones preeliminares Usaremos la base de datos “mtcars”
## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4
## Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4
## Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
## Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
## Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
## Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1- Lo que queremos hacer sería calcular la correlación entre mpg y wt
Primero se visualiza la relación en una gráfica de dispersión
## 'data.frame': 32 obs. of 11 variables:
## $ mpg : num 21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
## $ cyl : Factor w/ 3 levels "4","6","8": 2 2 1 2 3 2 3 1 1 2 ...
## $ disp: num 160 160 108 258 360 ...
## $ hp : num 110 110 93 110 175 105 245 62 95 123 ...
## $ drat: num 3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
## $ wt : num 2.62 2.88 2.32 3.21 3.44 ...
## $ qsec: num 16.5 17 18.6 19.4 17 ...
## $ vs : num 0 0 1 1 0 1 0 1 1 1 ...
## $ am : num 1 1 1 0 0 0 0 0 0 0 ...
## $ gear: num 4 4 4 3 3 3 3 4 4 4 ...
## $ carb: num 4 4 1 1 2 1 4 2 2 4 ...- Ver los datos con ggscatter para dispersión
ggscatter (my_data, x="wt", y= "mpg",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Peso (1000 libras)", ylab= "Millas por galón")## `geom_smooth()` using formula 'y ~ x'
* La relación es negativa, por tanto es inversamente propocional lo que significa que a mayor peso del motor, menor eficiencia en mpg