Regresion
Badouin
5/2/2022
setwd("~/Documents/pye1pm")
library(xfun)##
## Attaching package: 'xfun'## The following objects are masked from 'package:base':
##
## attr, isFALSERegresión lineal
La correlación lineal y la regresión lineal simple son métodos estadísticos que estudian la relación lineal existente entre dos variables. Antes de profundizar en cada uno de ellos, conviene destacar algunas diferencias:
La correlación cuantifica como de relacionadas están dos variables, mientras que la regresión lineal consiste en generar una ecuación (modelo) que, basándose en la relación existente entre ambas variables, permita predecir el valor de una a partir de la otra.
El cálculo de la correlación entre dos variables es independiente del orden o asignación de cada variable a X e Y , mide únicamente la relación entre ambas sin considerar dependencias. En el caso de la regresión lineal, el modelo varía según qué variable se considere dependiente de la otra (lo cual no implica causa-efecto)
La correlación no necesariamente implica causalidad
A nivel experimental, la correlación se suele emplear cuando ninguna de las variables se ha controlado, simplemente se han medido ambas y se desea saber si están relacionadas. En el caso de estudios de regresión lineal, es más común que una de las variables se controle (tiempo, concentración de reactivo, temperatura…) y se mida la otra.
Por norma general, los estudios de correlación lineal preceden a la generación de modelos de regresión lineal. Primero se analiza si ambas variables están correlacionadas y, en caso de estarlo, se procede a generar el modelo de regresión. ## Correlación lineal
Para estudiar la relación lineal existente entre dos variables continuas es necesario disponer de parámetros que permitan cuantificar dicha relación. Uno de estos parámetros es la covarianza, que indica el grado de variación conjunta de dos variables aleatorias.
\[ \text{Covarianza muestral} = Cov(X,Y)=\dfrac {\sum _{i=1}^{n}\left( x_{i}-\overline {x}\right) \left( y_{i}-\overline {y}\right) }{N-1} \]
siendo x⎯⎯⎯ e y⎯⎯⎯ la media de cada variable y xi e yi el valor de las variables para la observación i.
La covarianza depende de las escalas en que se miden las variables estudiadas, por lo tanto, no es comparable entre distintos pares de variables. Para poder hacer comparaciones se estandariza la covarianza, generando lo que se conoce como coeficientes de correlación. Existen diferentes tipos, de entre los que destacan el coeficiente de Pearson, Rho de Spearman y Tau de Kendall.
Todos ellos varían entre +1 y -1. Siendo +1 una correlación positiva perfecta y -1 una correlación negativa perfecta.
Se emplean como medida de fuerza de asociación (tamaño del efecto):
0: asociación nula.
0.1: asociación pequeña.
0.3: asociación mediana.
0.5: asociación moderada.
0.7: asociación alta.
0.9: asociación muy alta.
Las principales diferencias entre estos tres coeficientes de asociación son:
La correlación de Pearson funciona bien con variables cuantitativas que tienen una distribución normal. En el libro Handbook of Biological Statatistics se menciona que sigue siendo bastante robusto a pesar de la falta de normalidad. Es más sensible a los valores extremos que las otras dos alternativas.
La correlación de Spearman se emplea cuando los datos son ordinales, de intervalo, o bien cuando no se satisface la condición de normalidad para variables continuas y los datos se pueden transformar a rangos. Es un método no paramétrico.
La correlación de Kendall es otra alternativa no paramétrica para el estudio de la correlación que trabaja con rangos. Se emplea cuando se dispone de pocos datos y muchos de ellos ocupan la misma posición en el rango, es decir, cuando hay muchas ligaduras.
Además del valor obtenido para el coeficiente de correlación, es necesario calcular su significancia. Solo si el p-value es significativo se puede aceptar que existe correlación, y esta será de la magnitud que indique el coeficiente. Por muy cercano que sea el valor del coeficiente de correlación a +1 o −1 , si no es significativo, se ha de interpretar que la correlación de ambas variables es 0, ya que el valor observado puede deberse a simple aleatoriedad.
El test paramétrico de significancia estadística empleado para el coeficiente de correlación es el t-test. Al igual que ocurre siempre que se trabaja con muestras, por un lado está el parámetro estimado (en este caso el coeficiente de correlación) y por otro su significancia a la hora de considerar la población entera. Si se calcula el coeficiente de correlación entre X e Y en diferentes muestras de una misma población, el valor va a variar dependiendo de las muestras utilizadas. Por esta razón se tiene que calcular la significancia de la correlación obtenida y su intervalo de confianza.
\[ t=\dfrac {r\sqrt{N-2}}{\sqrt{1-r^{2}}} , \ \ \ df= N-2 \]
ara este test de hipótesis, H0 considera que las variables son independientes (coeficiente de correlación poblacional = 0) mientras que, la Ha , considera que existe relación (coeficiente de correlación poblacional ≠ 0) La correlación lineal entre dos variables, además del valor del coeficiente de correlación y de sus significancia, también tiene un tamaño de efecto asociado. Se conoce como coeficiente de determinación R2 . Se interpreta como la cantidad de varianza de Y explicada por X . En el caso del coeficiente de Pearson y el de Spearman, R2 se obtiene elevando al cuadrado el coeficiente de correlación. En el caso de Kendall no se puede calcular de este modo. (No he encontrado como se calcula).
Coeficiente de Pearson
El coeficiente de correlación de Pearson es la covarianza estandarizada, y su ecuación difiere dependiendo de si se aplica a una muestra, Coeficiente de Pearson muestral (r), o si se aplica la población Coeficiente de Pearson poblacional (ρ ).
Poblacional
\[ \rho=\dfrac {Cov(X,Y)}{\sigma_{x}\sigma_{y}} \]
Muestral
\[ r_{xy}=\dfrac {\sum _{i=1}^{n}\left( x_{i}-\overline {x}\right) \left( y_{i}-\overline {y}\right) } {\sqrt {\sum _{i=1}^{n}\left( x_{i}-\overline {x}\right) ^{2}\sum _{i=1}^{n}}\overline {\left( y_{i}-\overline {y}\right) ^{2}}} \]
Condiciones
La relación que se quiere estudiar entre ambas variables es lineal (de lo contrario, el coeficiente de Pearson no la puede detectar).
Las dos variables deben de ser cuantitativas.
Normalidad: ambas variables se tienen que distribuir de forma normal. Varios textos defienden su robustez cuando las variables se alejan moderadamente de la normal.
Homocedasticidad: La varianza de Y debe ser constante a lo largo de la variable X . Esto se puede identificar si en el scatterplot los puntos mantienen la misma dispersión en las distintas zonas de la variable X . Esta condición no la he encontrado mencionada en todos los libros.
Características
Toma valores entre [-1, +1], siendo +1 una correlación lineal positiva perfecta y -1 una correlación lineal negativa perfecta.
Es una medida independiente de las escalas en las que se midan las variables.
No varía si se aplican transformaciones a las variables.
No tiene en consideración que las variables sean dependientes o independientes.
El coeficiente de correlación de Pearson no equivale a la pendiente de la recta de regresión.
Es sensible a outliers, por lo que se recomienda en caso de poder justificarlos, excluirlos del análisis.
Interpretación
Además del valor obtenido para el coeficiente, es necesario calcular su significancia. Solo si el p-value es significativo se puede aceptar que existe correlación y esta será de la magnitud que indique el coeficiente. Por muy cercano que sea el valor del coeficiente de correlación a +1 o -1, si no es significativo, se ha de interpretar que la correlación de ambas variables es 0 ya que el valor observado se puede deber al azar. (Ver más adelante como calcular la significancia).
Ejemplo de modelo de regresión lineal para contenido de grasa en sangre de personas
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)
grasas## peso edad grasas
## 1 84 46 354
## 2 73 20 190
## 3 65 52 405
## 4 70 30 263
## 5 76 57 451
## 6 69 25 302
## 7 63 28 288
## 8 72 36 385
## 9 79 57 402
## 10 75 44 365
## 11 27 24 209
## 12 89 31 290
## 13 65 52 346
## 14 57 23 254
## 15 59 60 395
## 16 69 48 434
## 17 60 34 220
## 18 79 51 374
## 19 75 50 308
## 20 82 34 220
## 21 59 46 311
## 22 67 23 181
## 23 85 37 274
## 24 55 40 303
## 25 63 30 244Matriz de diagramas de dispersion
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)
Para cuantificar el grado de relación lineal a traves del coeficiente de correlacion de pearson, 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.0000000Cá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-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.3207x \]
Verificando el modelo para la edad de 23 años y con el contenido de grasa en sangre de 254, vamos a ver si el modelo tiene un resultado parecido a este dato real
y = 102.5751 + 5.3207 *23
y## [1] 224.9512El resultado que obtuvimos con el modelo es de: 224.9512 y el dato real de contenido en sangre para una persona de 23 años es de 254, el residuo de esto es:
residuo = 254 - 224.9512
residuo## [1] 29.0488Los 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$edad, grasas$grasas, xlab='Edad', ylab='Grasas')
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.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 de 30 hasta 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
## 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.6090Por 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:
\[ 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) 41.265155 163.885130
## edad 3.822367 6.818986Los 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 ' ' 1Diagnóstico del modelo
Los valores ajustados ŷ i y los residuos ei=ŷ 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)
¿Los datos son normales?
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:
qqnorm(residuos)
qqline(residuos)
Dado que los puntos están bastante alineados, la normalidad también parece aceptable.