Correlación y Regresión Lineal

Introducción
ANALISIS DE REGRESION
- SUPUESTOS DEL MODELO DE REGRESION LINEAL
- Regresion Multiple

Introducción

En ocasiones nos puede interesar estudiar si existe o no algún tipo de relación entre dos variables aleatorias. Las relaciones o asociaciones de este tipo entre variables se denominan correlaciones; y están medidas en escalas ordinales o de intervalos. En particular, nos interesa cuantificar la intensidad de la relación lineal entre dos variables.

Cuando el aumento de una de las variables viene acompañado del aumento de la otra, se trata de una correlación positiva o directa. La superficie de un bosque y el numero de nidos en el están correlacionados positivamente. Si el aumento de una variable viene acompañado de una disminución en la otra, la correlación es negativa o inversa.

El hecho de que dos variables estén correlacionadas no quiere decir que una sea la causante de la otra dos variables pueden estar independientemente relacionadas con una tercera (quizás no identificada).

Es importante notar que la existencia de correlación entre variables no implica causalidad.

Una Correlación perfecta se da cuando todos los puntos de un diagrama de dispersión forman una linea recta perfecta. Estas correlaciones perfectas (positivas o negativas) prácticamente no existen en biología; suelen ser privilegio exclusivo de los físicos.

En lenguaje llano, la palabra “correlación” describe cualquier tipo de relación entre objetos y hechos. No obstante, en estadística tiene un significado mas preciso: se refiere a la relación cuantitativa entre dos variables medidas en escalas ordinales o de intervalos.

Los coeficientes de correlación se pueden calcular mediante métodos paramétricos y no paramétricos. Un coeficiente paramétrico es el Coeficiente de Correlación de Pearson, que se usa para observaciones obtenidas sobre una escala de intervalos y esta sujeto a condiciones mas restrictivas que las alternativas no paramétricas. De estas, una de las mas ampliamente utilizadas es el Coeficiente de Correlación por Rangos de Spearman.

Si no hay correlación de ningún tipo entre dos variable, entonces tampoco habrá correlación lineal, por lo que r = 0. Sin embargo, el que ocurra r = 0 solo nos dice que no hay correlación lineal, pero puede que la haya de otro tipo. Mientras mas cercano a cero la correlación es mas débil.

alt text

En ocasiones, una correlación fuerte puede ser considerada como sistemáticamente no significativa, mientras que una débil puede ser sistemáticamente significativa, por lo que debemos resolver esta aparente paradoja (revisar, supuestos, tamaños de muestras, diseño estadístico y poder de la prueba, entre otras).

EL COEFICIENTE DE CORRELACION DE PEARSON

El Coeficiente de Correlación de Pearson es un estadístico paramétrico cuya aplicación es adecuada cuando las observaciones, de unidades maestreadas aleatoriamente, están medidas en escalas de intervalos. Se asume que ambas variables tienen una distribución aproximadamente normal, o sea, distribución normal bivariante. Esto puede comprobarse mediante un diagrama de dispersión de los datos, puesto que un diagrama de este tipo para datos normales bivariantes presenta un contorno aproximadamente circular o elíptico. El circulo se acerca mas a una elipse en tanto r aumenta su valor

peso <- c(51, 59, 49, 54, 50, 55, 48, 53, 52, 57)
long <- c(33.5, 38, 32, 37.5, 31.5, 33, 31, 36.5, 34, 35)
pairs(long ~ peso)  #permite elaborar un plot de correlacion

Otra forma de hacerlo de una forma elegante es a través del paquete “PerformanceAnalytics”

library(PerformanceAnalytics)
# Nuestros datos es mejor tenerlos en un data.frame
dat1 <- data.frame(peso, long)
chart.Correlation(dat1)

Calculo de la correlación.

cor(peso, long)  #cor sirve para analizar correlaciones tipo matriz

[1] 0.7794691

cor.test(peso, long)  # cor.test analiza la significancia de la correlacion, utilizando el contraste t-student.


    Pearson's product-moment correlation

data:  peso and long
t = 3.5194, df = 8, p-value = 0.007853
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.2942560 0.9452104
sample estimates:
      cor 
0.7794691

Antes de continuar, para saber que tipo correlación aplicar es importante revisar los supuestos de las pruebas.

shapiro.test(peso)


    Shapiro-Wilk normality test

data:  peso
W = 0.97343, p-value = 0.9207

shapiro.test(long)


    Shapiro-Wilk normality test

data:  long
W = 0.93757, p-value = 0.5263

Ambas variables presentan una distribución normal, por lo que Pearson es la prueba adecuada para nuestros datos.

EL COEFICIENTE DE DETERMINACION R2

Una pregunta importante que se plantea en el análisis de regresión es la siguiente: Que porcentaje de la variación total en Y se debe a la variación en X? En otras palabras, cual es la proporción de la variación total en Y que puede ser explicada por la variación en X?

El cuadrado del coeficiente de correlación de Pearson es, en si mismo, un estadístico de gran utilidad denominado coeficiente de determinación. Es una medida de la proporción de la variabilidad en una variable, atribuible a la variabilidad de la otra. En una correlación perfecta donde r = +1 o -1, la variación de variable se corresponde con una variación exacta en el valor de la otra. Esta situación no se suele dar en biología, puesto que son muchos los factores que, por lo general, regulan las relaciones entre variables entre organismos.

En el ejemplo anterior el coeficiente de determinación es 0.6075721 = 0,60.75 o 60.75%, de lo que se deduce que aproximadamente el 40% del peso de las aves no esta relacionado con la longitud del pico.

0.7794691^2

[1] 0.6075721

EL COEFICIENTE DE CORRELACION POR RANGOS DE SPEARMAN

Cuando exista alguna duda acerca de si se cumplen las relativamente rigurosas premisas para la aplicación del coeficiente de correlación de Pearson, entonces debemos considerar el uso de alguna alternativa no paramétrica.

Esta prueba estadística permite medir la correlación o asociación de dos variables y es aplicable cuando las mediciones se realizan en una escala ordinal, aprovechando la clasificación por rangos.

El coeficiente de correlación de Spearman se rige por las reglas de la correlación simple de Pearson, y las mediciones de este indice corresponden de + 1 a - 1, pasando por el cero, donde este ultimo significa no correlación entre las variables estudiadas, mientras que los dos primeros denotan la correlación máxima.

x <- c(44.4, 45.9, 41.9, 53.3, 44.7, 44.1, 40.7, 45.2, 60.1)
y <- c(2.6, 3.1, 2.5, 3, 3.6, 4, 3.2, 2.8, 3.8)
shapiro.test(x)


    Shapiro-Wilk normality test

data:  x
W = 0.80304, p-value = 0.02212

shapiro.test(y)


    Shapiro-Wilk normality test

data:  y
W = 0.94856, p-value = 0.6743

Una de las dos variables es asimétrica, por lo que tendremos que utilizar una correlación no paramétrica de Spearman

cor(x, y, method = "spearman")

[1] 0.1666667

cor.test(x, y, method = "spearman")


    Spearman's rank correlation rho

data:  x and y
S = 100, p-value = 0.6777
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.1666667

revise los argumentos de esas pruebas , ?cor.test.

Utilice la base de datos “women”. Que tipo de correlacion debe aplicar para mostrar si existe una asociacion entre la altura (height) y el peso (weight)

data(women)
head(women)

  height weight
1     58    115
2     59    117
3     60    120
4     61    123
5     62    126
6     63    129

Se recomienda dar un vista al paquete “corrplot”, dado que permite ejecutar gráficas y correlaciones muy elegantes para ser usadas en publicación o presentaciones de sus trabajos. Además de que permite al lector hacer una lectura rápida del comportamiento de las variables en cuanto a su asociación.

Ejemplo de como funciona

library(corrplot)

data(mtcars)
head(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

M <- cor(mtcars)  #permite ejecutar una matriz de correlación

Aplicación gráfica

corrplot(M, method = "ellipse")

# Visualice otros metodos gráficos como, 'circle', 'square', 'ellipse',
# 'number', 'shade', 'color', 'pie',

ANALISIS DE REGRESION

El Análisis de Regresión se usa cuando el investigador sabe que existe una relación entre las variables porque hay una teoría o investigaciones previas que la han descubierto. Por ejemplo, la relación entre espacio y tiempo ya se sabe que es la velocidad, o como la relación entre voltaje e intensidad de corriente eléctrica. En estos casos, el investigador suele estar interesado en verificar experimentalmente tal relación y el objeto de la regresión es encontrar la curva que mejor ajuste a sus datos experimentales.

Utilizando un diagrama de dispersión de puntos podemos resaltar las características de dos variables relacionadas, la linea que se genera denomina linea de mejor ajuste.

SUPUESTOS DEL MODELO DE REGRESION LINEAL

alt text

Ejemplo 1

peso <- c(61, 60, 78, 62, 66, 60, 54, 84, 68)
altura <- c(162, 154, 180, 158, 171, 169, 166, 176, 163)
ej1 <- data.frame(peso, altura)
reg1 <- lm(peso ~ altura, data = ej1)
summary(reg1)


Call:
lm(formula = peso ~ altura, data = ej1)

Residuals:
    Min      1Q  Median      3Q     Max 
-11.444  -3.447   1.347   4.164  10.549 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) -67.4679    51.1995  -1.318    0.229  
altura        0.8007     0.3071   2.608    0.035 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.268 on 7 degrees of freedom
Multiple R-squared:  0.4927,    Adjusted R-squared:  0.4203 
F-statistic: 6.799 on 1 and 7 DF,  p-value: 0.03504

anova(reg1)

Analysis of Variance Table

Response: peso
          Df Sum Sq Mean Sq F value  Pr(>F)  
altura     1 359.15  359.15  6.7994 0.03504 *
Residuals  7 369.74   52.82                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Observar que el anova es significativo. Muy importante tener este resultado en las pruebas de la regresión, dado que esto se ajusta a lo que buscamos en la hipótesis de la regresión.

Ho = pendientes es igual a cero

H1 = pendientes es diferente a cero

Verificacion de Supuestos

par(mfrow = c(2, 2))
plot(reg1)

# Residual vs Fitted: Deberia estar distribuidos aleatoriamente alrededor de
# la linea horizontal que representa un error residual de cero

# Normal Q-Q: deberia sugerir que los errores residuales se distribuyen
# normalmente.

# Scale-Location Muestra la raiz cuadrada de los residuos estandarizados,
# como una funcion de los valores ajustados. No deberia existir una
# tendencia clara en ese trama.

# Residual vs Leverage Las distancias mas grandes que 1 son sospechosos y
# sugieren la presencia de un valor atipico posible. y su eliminacion podria
# tener efectos sobre la regresion.

Forma grafica de la regresion

Exploraremos el uso de dos paquetes visreg y ggplot2 para traficar. Usted decide cual es su preferido.

library(visreg)
# Notese que solo requiere incluir la variable independiente, para graficar.
visreg(reg1, "altura", partial = F)  #reg1 es nuestro vector que contiene la ecuacion de regresion.

# Nos brinda una resumen de la regresion
summary(reg1)


Call:
lm(formula = peso ~ altura, data = ej1)

Residuals:
    Min      1Q  Median      3Q     Max 
-11.444  -3.447   1.347   4.164  10.549 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) -67.4679    51.1995  -1.318    0.229  
altura        0.8007     0.3071   2.608    0.035 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.268 on 7 degrees of freedom
Multiple R-squared:  0.4927,    Adjusted R-squared:  0.4203 
F-statistic: 6.799 on 1 and 7 DF,  p-value: 0.03504

Los valores que debemos interpretar del summary son:

alt text

La ecuacion de prediccion obtenida por la regresion seria:

==========================

peso = -67.4679 + 0.8007 * altura

==========================

La interpretacion de los resultados es:

En promedio, cada incremento en una unidad de la altura, corresponde a un incremento del peso de 0.80 kg, asi mismo que una persona cuyo peso sea de cero (origen de y), tenga una altura de -67.46, esto obviamente no parece posible. Dado que la regresion estima como logramos ver, estimaciones mas alla de los valores observados.

#Otra forma grafica con ggplot2
library(ggplot2)#ggplot2 es una extension poderosa para graficar
ggplot(ej1, aes(x=altura, y=peso)) +
geom_point(shape=1) +    # genera circulos en el grafico
geom_smooth(method=lm)   # adjunta la linea de regresion por defecto es al 95% de confianza

Prueba sobre los residuos

shapiro.test(resid(reg1))  #Normalidad


    Shapiro-Wilk normality test

data:  resid(reg1)
W = 0.97706, p-value = 0.9475

library(car)
ncvTest(reg1)  #Prueba de heterocedasticidad, Ho=varianzas son constantes

Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 0.3545291    Df = 1     p = 0.5515605

Revisar el vídeo tutorial de, Regresión Lineal Simple

Regresion Multiple

# Usaremos los datos de Prestige
library(car)
data(Prestige)
head(Prestige)

                    education income women prestige census type
gov.administrators      13.11  12351 11.16     68.8   1113 prof
general.managers        12.26  25879  4.02     69.1   1130 prof
accountants             12.77   9271 15.70     63.4   1171 prof
purchasing.officers     11.42   8865  9.11     56.8   1175 prof
chemists                14.62   8403 11.68     73.5   2111 prof
physicists              15.64  11030  5.13     77.6   2113 prof

names(Prestige)

[1] "education" "income"    "women"     "prestige"  "census"    "type"

# Generamos un modelo de regresion multiple (variables independientes)
reg2 <- lm(prestige ~ education + women, data = Prestige)
summary(reg2)


Call:
lm(formula = prestige ~ education + women, data = Prestige)

Residuals:
    Min      1Q  Median      3Q     Max 
-28.010  -4.069   1.050   5.027  18.942 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -8.75416    3.54213  -2.471 0.015164 *  
education    5.42780    0.31612  17.170  < 2e-16 ***
women       -0.09305    0.02719  -3.422 0.000904 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.652 on 99 degrees of freedom
Multiple R-squared:  0.7521,    Adjusted R-squared:  0.7471 
F-statistic: 150.2 on 2 and 99 DF,  p-value: < 2.2e-16

# Revisamos el modelo
anova(reg2)

Analysis of Variance Table

Response: prestige
          Df  Sum Sq Mean Sq F value    Pr(>F)    
education  1 21608.4 21608.4 288.685 < 2.2e-16 ***
women      1   876.7   876.7  11.713 0.0009039 ***
Residuals 99  7410.3    74.9                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Verificamos los supuestos
par(mfrow = c(2, 2))
plot(reg2)

# Revisamos la distribucion de los residuos
shapiro.test(residuals(reg2))


    Shapiro-Wilk normality test

data:  residuals(reg2)
W = 0.98367, p-value = 0.2419

# Forma grafica con visreg
library(visreg)
visreg(reg2, "education", partial = F)

visreg(reg2, "women", partial = F)

#Forma grafica con con ggplot2, para la variable women
library(ggplot2)
ggplot(Prestige, aes(x=women, y=prestige)) +
geom_point(shape=1) +    # genera circulos en el grafico
geom_smooth(method=lm)   # adjunta la linea de regresion por defecto es al 95% de confianza

#Forma grafica con con ggplot2, para la variable education
library(ggplot2)
ggplot(Prestige, aes(x=education, y=prestige)) +
geom_point(shape=1) +    # genera circulos en el grafico
geom_smooth(method=lm)   # adjunta la linea de regresion por defecto es al 95% de confianza

#Figura con lineas suavizadas
library(ggplot2)
par(mfrow=c(1,2))
ggplot(Prestige, aes(x=women, y=prestige)) +
geom_point(shape=1) +    # genera circulos en el grafico
geom_smooth()  # adjunta la linea de regresion por defecto es al 95% de confianza, de forma suavizada.

ggplot(Prestige, aes(x=income, y=prestige)) +
geom_point(shape=1) +    # genera circulos en el grafico
geom_smooth()  # adjunta la linea de regresion por defecto es al 95% de confianza, de forma suavizada.

# Generamos un segundo modelo de regresion multiple, incluyendo una variable
# adicional (income) al modelo
reg3 <- lm(prestige ~ education + women + income, data = Prestige)
# Revisamos todo lo aprendido
summary(reg3)


Call:
lm(formula = prestige ~ education + women + income, data = Prestige)

Residuals:
     Min       1Q   Median       3Q      Max 
-19.8246  -5.3332  -0.1364   5.1587  17.5045 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -6.7943342  3.2390886  -2.098   0.0385 *  
education    4.1866373  0.3887013  10.771  < 2e-16 ***
women       -0.0089052  0.0304071  -0.293   0.7702    
income       0.0013136  0.0002778   4.729 7.58e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.846 on 98 degrees of freedom
Multiple R-squared:  0.7982,    Adjusted R-squared:  0.792 
F-statistic: 129.2 on 3 and 98 DF,  p-value: < 2.2e-16

anova(reg3)

Analysis of Variance Table

Response: prestige
          Df  Sum Sq Mean Sq F value    Pr(>F)    
education  1 21608.4 21608.4 350.974 < 2.2e-16 ***
women      1   876.7   876.7  14.240 0.0002757 ***
income     1  1376.7  1376.7  22.361 7.579e-06 ***
Residuals 98  6033.6    61.6                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

par(mfrow = c(2, 2))
plot(reg3)

shapiro.test(residuals(reg3))


    Shapiro-Wilk normality test

data:  residuals(reg3)
W = 0.99363, p-value = 0.918

library(visreg)
visreg(reg3, "income", partial = F)

# Cual modelos es el mejor?
AIC(reg2, reg3)

     df      AIC
reg2  4 734.5998
reg3  5 715.6358

# Buscamos las variables influyentes en el modelo. Esto nos premite ver las
# variables que afectan el modelo de regresion.
avPlots(reg2, id.n = 2, id.cex = 0.7)

avPlots(reg3, id.n = 2, id.cex = 0.7)

Revisar el vídeo tutorial de, Como eliminar datos de un data.frame (matriz de datos) en Rstudio

# Outlier, nos permite identificar los valores extremos.
outlierTest(reg2)


No Studentized residuals with Bonferonni p < 0.05
Largest |rstudent|:
          rstudent unadjusted p-value Bonferonni p
newsboys -3.437184         0.00086402      0.08813

qqPlot(reg2, id.n = 3)

                 newsboys service.station.attendant 
                        1                         2 
       electronic.workers 
                      102

outlierTest(reg3)


No Studentized residuals with Bonferonni p < 0.05
Largest |rstudent|:
          rstudent unadjusted p-value Bonferonni p
newsboys -2.694442          0.0083101      0.84763

qqPlot(reg3, id.n = 3)

          newsboys electronic.workers            farmers 
                 1                101                102

ncvTest(reg3)  # Test de heretocesdasticidad. H0= varianza constante.

Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 0.2490143    Df = 1     p = 0.61777

vif(reg3)  #Test de multicolineariedad, valores >4 sugieren multicolineariedad

education     women    income 
 1.845165  1.526593  2.282038

###FIN