Clase 4

Carga de librerías

library(tidyverse)
library(GGally)
library(cowplot)

Motivación y objetivo

Nuestro objetivo va a ser crear un modelo estadístico sencillo que nos permita modelar una relación lineal entre dos variables: una será nuestra variable a explicar y la otra será nuestra variable explicativa. Para eso vamos a ver:

1. El concepto de covarianza y correlación, su estimación y los tests asociados
2. Las diferencias entre la correlación y causalidad
3. Modelo de Regresión Lineal Simple: breve introducción teórica, interpretación y evaluación.

Dataset

Vamos a trabajar con el dataset state.x77 que se encuentra en R. Como es una matriz, lo transformamos a un dataframe con as.data.frame()

estados <- state.x77 %>%
  as.data.frame() %>% # Transformo la matriz a dataframe
  rename(life_exp = `Life Exp`, hs_grad=`HS Grad`) # Renombro las variables para no tener problema con los espacios

Para conocer el dataset utilicemos ?state.x77 para abrir la ayuda de R y la funcion glimpse() para ver la estructura de nuestros datos.

?state.x77
glimpse(estados)

Observations: 50
Variables: 8
$ Population <dbl> 3615, 365, 2212, 2110, 21198, 2541, 3100, 579, 8277, 4931, 868, 813, 11197, 5313, 2861, 2280, 3387, 3806, 1058, 41...
$ Income     <dbl> 3624, 6315, 4530, 3378, 5114, 4884, 5348, 4809, 4815, 4091, 4963, 4119, 5107, 4458, 4628, 4669, 3712, 3545, 3694, ...
$ Illiteracy <dbl> 2.1, 1.5, 1.8, 1.9, 1.1, 0.7, 1.1, 0.9, 1.3, 2.0, 1.9, 0.6, 0.9, 0.7, 0.5, 0.6, 1.6, 2.8, 0.7, 0.9, 1.1, 0.9, 0.6,...
$ life_exp   <dbl> 69.05, 69.31, 70.55, 70.66, 71.71, 72.06, 72.48, 70.06, 70.66, 68.54, 73.60, 71.87, 70.14, 70.88, 72.56, 72.58, 70...
$ Murder     <dbl> 15.1, 11.3, 7.8, 10.1, 10.3, 6.8, 3.1, 6.2, 10.7, 13.9, 6.2, 5.3, 10.3, 7.1, 2.3, 4.5, 10.6, 13.2, 2.7, 8.5, 3.3, ...
$ hs_grad    <dbl> 41.3, 66.7, 58.1, 39.9, 62.6, 63.9, 56.0, 54.6, 52.6, 40.6, 61.9, 59.5, 52.6, 52.9, 59.0, 59.9, 38.5, 42.2, 54.7, ...
$ Frost      <dbl> 20, 152, 15, 65, 20, 166, 139, 103, 11, 60, 0, 126, 127, 122, 140, 114, 95, 12, 161, 101, 103, 125, 160, 50, 108, ...
$ Area       <dbl> 50708, 566432, 113417, 51945, 156361, 103766, 4862, 1982, 54090, 58073, 6425, 82677, 55748, 36097, 55941, 81787, 3...

Tenemos 50 observaciones (los 50 estados de Estados Unidos) y 8 variables númericas:

• population: población
• income: ingreso per capita
• illiteracy: porcentaje de la población que es analfabeta
• Life Exp: esperanza de vida en años
• Murder: cantidad de homicidios cada 100.000 habitantes
• HS Grad: porcentaje de la población que terminó la secundaria
• Frost: promedio de días con temperatura mínima por debajo de los 0 grados
• Area: area en millas cuadradas

Veamos algunas estadisticas de resumen para conocer un poco más nuestras variables, usando la funcion summary:

summary(estados)

   Population        Income       Illiteracy       life_exp         Murder          hs_grad          Frost             Area       
 Min.   :  365   Min.   :3098   Min.   :0.500   Min.   :67.96   Min.   : 1.400   Min.   :37.80   Min.   :  0.00   Min.   :  1049  
 1st Qu.: 1080   1st Qu.:3993   1st Qu.:0.625   1st Qu.:70.12   1st Qu.: 4.350   1st Qu.:48.05   1st Qu.: 66.25   1st Qu.: 36985  
 Median : 2838   Median :4519   Median :0.950   Median :70.67   Median : 6.850   Median :53.25   Median :114.50   Median : 54277  
 Mean   : 4246   Mean   :4436   Mean   :1.170   Mean   :70.88   Mean   : 7.378   Mean   :53.11   Mean   :104.46   Mean   : 70736  
 3rd Qu.: 4968   3rd Qu.:4814   3rd Qu.:1.575   3rd Qu.:71.89   3rd Qu.:10.675   3rd Qu.:59.15   3rd Qu.:139.75   3rd Qu.: 81162  
 Max.   :21198   Max.   :6315   Max.   :2.800   Max.   :73.60   Max.   :15.100   Max.   :67.30   Max.   :188.00   Max.   :566432  

Podemos ver el valor mínimo, primer cuartil, mediana, promedio, tercer cuartil y valor máximo de nuestras variables numéricas

Covarianza y correlación

Covarianza

Definición

Es un estadístico que permite medir la variabilidad conjunta de 2 variables.

Puntualmente, nos permite medir la asociacion lineal entre dos variables.

La fórmula de cálculo de esta medida es:

\(q= \frac{1}{N-1} \sum\limits_{i=1}^N (x_i-\bar{x})(y_i-\bar{y})\)

• Cuando valores altos de x corresponden a valores altos de y, y valores bajos de x corresponden a valores bajos de y la COVARIANZA es positiva

• Cuando valores altos de x corresponden a valores bajos de y, y valores bajos de x corresponden a valores altos de y la COVARIANZA es negativa

• Cuando no ocurre ninguna de las dos cosas anteriores, la covarianza será muy cercana a cero

Interpretacion grafica

A continuación presentamos algunos gráficos de dispersión entre variables del dataset de estados.

ASOCIACION LINEAL POSITIVA

• En el eje horizontal se encuentra la variable Asesinatos (Murder) y en el eje vertical se encuentra la variable Analfabetismo (Illiteracy)

• Un punto representa las mediciones de esas dos variables para un estado

• La línea vertical punteada marca el promedio de Asesinatos (Murder)

• La línea horizontal punteada marca el promedio de Analfabetismo (Illiteracy)

ggplot(estados, aes(x=Murder,y=Illiteracy)) + geom_point(size=2) +
  geom_hline(yintercept = mean(estados$Illiteracy),color='steelblue', linetype='dashed', size=1) +
  geom_vline(xintercept = mean(estados$Murder), color='steelblue', linetype='dashed', size=1) +
  annotate(geom='text', x=c(6.5,8.3,15,15), y=c(2.9,2.9,1.5,0.8), label='+', colour='forestgreen', size=19) +
  annotate(geom='text', x=c(6.5,8.3,0,0), y=c(0.3,0.3,1.5,0.8), label='-', colour='firebrick', size=27) +
  annotate(geom='text', x=15, y=2.9, label='1', colour='steelblue', size=10) +
  annotate(geom='text', x=15, y=0.3, label='4', colour='steelblue', size=10) +
  annotate(geom='text', x=0, y=0.3, label='3', colour='steelblue', size=10) +
  annotate(geom='text', x=0, y=2.9, label='2', colour='steelblue', size=10) +
  labs(title='Asociacion lineal positiva', x='Asesinatos (Murder)', y='Analfabetismo (Illiteracy)')+
  theme_bw() + scale_y_continuous(limits = c(0,3)) + scale_x_continuous(limits = c(0,16))

Vemos que los puntos se concentran en:

• En el cuadrante 1) donde el valor de cada variable es mayor a su promedio. Entonces, la fórmula de la covarianza es positiva

• En el cuadrante 3) donde el valor de cada variable es menor a su promedio. Entonces, la fórmula de la covarianza tambien es positiva

Por lo tanto, la covarianza es positiva:

La funcion cov permite calcular la covarianza entre dos variables.

cov(x=estados$Murder, y=estados$Illiteracy)

[1] 1.581776

ASOCIACION LINEAL NEGATIVA

• En el eje horizontal se encuentra la variable Asesinatos (Murder) y en el eje vertical se encuentra la variable Life Exp (Experanza de Vida)

• Un punto representa las mediciones de esas dos variables para un estado

• La línea vertical punteada marca el promedio de Asesinatos (Murder)

• La línea horizontal punteada marca el promedio de Life Exp (Experanza de Vida)

ggplot(estados, aes(x=Murder,y=life_exp)) + geom_point(size=2) +
  geom_hline(yintercept = mean(estados$life_exp), color='steelblue', linetype='dashed', size=1, alpha=0.7) +
  geom_vline(xintercept = mean(estados$Murder), color='steelblue', linetype='dashed', size=1, alpha=0.7) +
  annotate(geom='text', x=c(6.5,8.3,15,15), y=c(74.2,74.2,71.7,70.1), label='+', colour='forestgreen', size=19) +
  annotate(geom='text', x=c(6.5,8.3,0,0), y=c(67.9,67.9,71.7,70.1), label='-', colour='firebrick', size=27) +
  annotate(geom='text', x=15, y=74.2, label='1', colour='steelblue', size=10) +
  annotate(geom='text', x=15, y=67.9, label='4', colour='steelblue', size=10) +
  annotate(geom='text', x=0, y=67.9, label='3', colour='steelblue', size=10) +
  annotate(geom='text', x=0, y=74.2, label='2', colour='steelblue', size=10) +
  labs(title='Asociacion lineal negativa', x='Asesinatos (Murder)', y='Life Exp (Experanza Vida)')+
  theme_bw() +
  scale_x_continuous(limits = c(0,16)) + scale_y_continuous(limits = c(67.5,74.5))

Vemos que los puntos se concentran en:

• En el cuadrante 2) donde el valor de Asesinatos es menor a su promedio y el valor de Esperanza de Vida es mayor a su promedio. Entonces, la fórmula de la covarianza tambien es positiva

• En el cuadrante 4) donde el valor de Asesinatos es mayor a su promedio y el valor de Esperanza de Vida es menor a su promedio. Entonces, la fórmula de la covarianza es negativa

Por lo tanto, la covarianza es negativa:

cov(x=estados$Murder, y=estados$life_exp)

[1] -3.86948

SIN ASOCIACION LINEAL

• En el eje horizontal se encuentra la variable Area y en el eje vertical se encuentra la variable Analfabetismo (Illiteracy)

• Un punto representa las mediciones de esas dos variables para un estado

• La línea vertical punteada marca el promedio de Area

• La línea horizontal punteada marca el promedio de Analfabetismo (Illiteracy)

ggplot(estados, aes(x=Area,y=Illiteracy)) + geom_point(size=2) +
  geom_hline(yintercept = mean(estados$Illiteracy),color='steelblue', linetype='dashed', size=1) +
  geom_vline(xintercept = mean(estados$Area), color='steelblue', linetype='dashed', size=1) +
  annotate(geom='text', x=c(-10000,150000,500000,500000), y=c(2.9,2.9,1.5,0.8), label='+', colour='forestgreen', size=19) +
  annotate(geom='text', x=c(0,150000,-400000,-400000), y=c(-0.6,-0.6,1.5,0.8), label='-', colour='firebrick', size=27) +
  annotate(geom='text', x=500000, y=2.9, label='1', colour='steelblue', size=10) +
  annotate(geom='text', x=500000, y=-0.6, label='4', colour='steelblue', size=10) +
  annotate(geom='text', x=-400000, y=-0.6, label='3', colour='steelblue', size=10) +
  annotate(geom='text', x=-400000, y=2.9, label='2', colour='steelblue', size=10) +
  labs(title='Sin asociacion lineal',y='Analfabetismo (Illiteracy)')+
  theme_bw() + scale_y_continuous(limits = c(-1,3)) + scale_x_continuous(limits = c(-400000,600000))

Vemos que los puntos se distribuyen entre el cuadrante 2) y 3) sin una estructura muy clara. Entonces decimos que no hay asociación lineal entre las variables:

cov(x = estados$Area, y=estados$Illiteracy)

[1] 4018.337

¿Qué pasó? Habíamos dicho que la covarianza entre Asesinatos y Analfabetismo era positiva y la covarianza entre Area y Analfabetismo era inexsistente. Sin embargo, la segunda es mucho más grande que la primera.

La covarianza tiene una CARACTERÍSTICA que puede ser un PROBLEMA importante: se ve afectada por la unidad de medida de las variables.

Por ejemplo:

# Covarianza entre area (medida en millas cuadradas) y analfabetismo
cov_1 = cov(x = estados$Area, y=estados$Illiteracy)
cov_1

[1] 4018.337

# 1 milla cuadrada = 2.59 kilometros cuadrados
area_kilometros = estados$Area*2.59
# Covarianza entre area (medida en kilometros cuadradas) y analfabetismo
cov_2=cov(x = area_kilometros, y=estados$Illiteracy)
cov_2

[1] 10407.49

La covarianza entre el área y analfabetismo es distinta dependiendo de como midamos el área. Si la medimos en millas cuadras la covarianza es 4018.3 y si la medimos en kilómetros cuadrados es igual a 10407.5.

Como la covarianza se ve afectada por la unidad de medida de las variables puede ser muy difícil (o imposible) realizar una comparación entre las covarianzas de distintas variables. Para eso podemos utilizar otra medida

Correlación

Definición

Al igual que la covarianza, la correlación nos permite medir la asociación lineal entre dos variables.

La fórmula de calculo para el coeficiente de correlación es:

\(r= \frac{1}{N-1} \sum\limits_{i=1}^N \frac{(x_i-\bar{x})(y_i-\bar{y})}{S_x \cdot S_y}\)

\(S_x\) es el desvío estándar estimado de X y \(S_y\) es el desvío estándar estimado de Y.

Al agregar los desvíos, la correlación sólo puede tomar valores entre -1 y 1 y así se resuelve el problema que tiene la covarianza de verse afectada por la unidad de medida de las variables

Características

1. Puede tomar valores entre -1 y 1

2. No se ve afectada por las unidades de medida de las variables

3. El valor absoluto (módulo) del coeficiente mide la fuerza de la relacion lineal entre X e y. Cuanto mayor sea el valor absoluto, más fuerte es la relación LINEAL entre X e Y

4. 1. \(r=1\) indica una relación lineal perfecta positiva (los puntos se encuentra una recta diagonal de pendiente positiva)
   2. \(r=0\) indica que no existe relación lineal entre X e Y
   3. \(r=-1\) indica una relación lineal perfecta negativa (los puntos se encuentra una recta diagonal de pendiente negativa)

Veamos la caracteristica 2) con nuestro ejemplo de la relacion entre analfabetismo y area

La funcion cor permite calcular la correlacion entre dos variables.

# Correlacion entre area (medida en millas cuadradas) y analfabetismo
corr_1 = cor(x = estados$Area, y=estados$Illiteracy)
corr_1

[1] 0.07726113

# 1 milla cuadrada = 2.59 kilometros cuadrados
area_kilometros = estados$Area*2.59
# Correlacion entre area (medida en kilometros cuadradas) y analfabetismo
corr_2=cor(x = area_kilometros, y=estados$Illiteracy)
corr_2

[1] 0.07726113

Vemos que ambos coeficientes de correlación son iguales.

Ejemplos

Veamos el coeficiente de correlación para las relaciones que graficamos antes:

# Correlacion entre homicidios y analfabetismo
corr_positiva = cor(x = estados$Murder, y=estados$Illiteracy)
corr_positiva

[1] 0.7029752

Existe una correlación positiva fuerte entre la cantidad de homicidios y el porcentaje de analfabetas

# Correlacion entre homicidios y esperanza de vida
corr_negativa = cor(x = estados$Murder, y=estados$life_exp)
corr_negativa

[1] -0.7808458

Existe una correlación negativa fuerte entre la cantidad de homicidios y la esperanza de vida

# Covarianza entre area (medida en kilometros cuadradas) y analfabetismo
corr_nula = cor(x = estados$Area, y=estados$Illiteracy)
corr_nula

[1] 0.07726113

Existe una correlación positiva muy débil entre área y analfabetismo

GGAlly

La libreria GGally es muy útil porque nos permite resumir información sobre las relaciones de las variables de manera gráfica.

Con el comando ggpairs vamos a poder ver:

1. En la diagonal: los gráficos de distribucion de la variable

2. Por debajo de la diagonal: los gráficos de dispersion para cada par de variables

3. Por encima de la diagonal: los coeficientes de correlacion para cada par de variables

ggpairs(estados) + theme_bw() + theme(axis.text.x = element_text(angle = 90, hjust = 1))

Por ejemplo:

• En la fila 2, columna 2 vemos la distribución de la variable ingreso

• En la fila 2, columna 1 vemos el diagrama de dispersión entre población e ingreso

• En la fila 2, columna 3 vemos el coeficiente de correlación entre el ingreso y el porcentaje de analfabetas

Test de Hipótesis

Motivación

El coeficiente de correlación es el resultado de una estimación sobre nuestro conjunto de datos. Por lo tanto, el número que se obtiene va a depender de que individuos fueron seleccionados en la muestra.

Veamos esto con un ejemplo

# Sacamos el promedio de días con temperatura mínima por debajo de los 0 grados para los 50 estados
mean(estados$Frost)

[1] 104.46

# Sacamos la correlacion entre dias debajo de los 0 grados (Frost) y el analfabetismo 
cor(estados$Frost, estados$Illiteracy)

[1] -0.671947

#Seleccionamos los 10 estados mas frios 
frios= estados %>% arrange(desc(Frost)) %>% slice(1:10)
# Sacamos el promedio de días con temperatura mínima por debajo de los 0 grados para esos 10 estados
mean(frios$Frost)

[1] 170.3

# Sacamos la correlacion entre dias debajo de los 0 grados (Frost) y el analfabetismo para esos 10 estados
cor(frios$Frost, frios$Illiteracy)

[1] 0.02294539

# Ahora seleccionamos los 10 estados mas calidos 
calidos= estados %>% arrange(desc(Frost)) %>% slice(41:50)
# Sacamos el promedio de días con temperatura mínima por debajo de los 0 grados para esos 10 estados
mean(calidos$Frost)

[1] 23.9

cor(calidos$Frost, calidos$Illiteracy)

[1] -0.1864304

Considerando los 50 estados vemos que en promedio tienen 104 días por debajo de 0 grados y la correlación entre días fríos y analfabetismo es negativa y bastante fuerte (igual a -0.67)

Los 10 estados más fríos tienen en promedio 170 días por debajo de los 0 grados (casi la mitad del año) mientras los 10 estados más cálidos tienen en promedio 24 días por debajo de los 0 grados (menos de un mes).

Por su parte, vemos que la correlación entre días fríos y analfabetismo es muy cercana a 0 (es igual a 0.02) para los estados más fríos, mientras que para los estados más cálidos la correlación es negativa y baja (es igual a -0.19).

Vemos que el coeficiente de correlación puede cambiar mucho dependiendo de que individuos hayamos incluido en nuestro dataset (muestra). Esto va suceder con cualquier medida estadística que calculemos sobre nuestro conjunto de datos.

En conclusión, el coeficiente de correlación (y otras medidas estadísticas como el promedio) va a tener cierta VARIABILIDAD O INCERTIDUMBRE en su estimación.

Test de Hipótesis

Puntualmente, nos interesa saber:

1. ¿Puede ser que el coeficiente de correlación sea igual a cero?

2. ¿Qué tan probable es que sea igual a cero?

Porque, si es probable que el coeficiente de correlacion sea igual a cero, entonces no va a existir una relacion lineal importante entre esas dos variables

En la siguiente imagen vemos:

• La curva ROJA es la distribución del coeficiente de correlación entre Area y Analfabetismo.

• La curva VERDE es es la distribución del coeficiente de correlación entre Asesinatos y Analfabetismo

• La linea punteada nos indica el punto en el cual el coeficiente de correlación es cero.

¿Cuál de estas dos curvas está más cerca del cero? ¿Cuál creen que es el coeficiente de correlación que es más probable que sea igual a cero?

ggplot(data = data.frame(x = c(-1, 1)), aes(x)) +
  stat_function(fun = dnorm, n = 100, args = list(mean = 0.07, sd = .1), color="firebrick", size=1.5) +
  stat_function(fun = dnorm, n = 100, args = list(mean = 0.71, sd = .1), color="forestgreen",size=1.5) +
  geom_vline(xintercept = 0, linetype='dashed', size=1) +
  theme_bw() +
  labs(title="Distribucion coeficientes de correlacion", x="Coeficiente de Correlacion", y="")+
  scale_y_continuous(breaks = NULL)

A simple vista, parece ser que el coeficiente de correlacion entre Area y Analfabetismo tiene una probabilidad más alta de ser igual a cero ya que la curva roja está centrada más cerca del cero. Mientras que la curva verde se encuentra bastante alejada del cero, por lo tanto la probabilidad que el coeficiente sea igual a cero es muy baja.

Podemos usar un test estadístico para responder correctamente a las dos preguntas que hicimos arriba. Vamos a TESTEAR la HIPÓTESIS: El coeficiente de correlacion es igual a cero

La función que vamos a utilizar se llama cor.test

# Veamos la documentacion
?cor.test
# Al igual que para calcular el coeficiente de correlacion pasamos nuestras dos variables: X e Y
cor.test(x=estados$Murder, y=estados$Illiteracy)

    Pearson's product-moment correlation

data:  estados$Murder and estados$Illiteracy
t = 6.8479, df = 48, p-value = 1.258e-08
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.5279280 0.8207295
sample estimates:
      cor 
0.7029752 

Observamos:

1. Debajo de la parte que dice 95 percent confidence interval: el intervalo de confianza para el coeficiente de correlación.

• El primer valor es el valor más chico que esperamos ver del coeficiente de correlación
• El segundo valor es el valor más grande que esperamos ver del coeficiente de correlación

IMPORTANTE: nos interesa ver si dentro del intervalo de confianza se encuentra el cero:

• Si esta DENTRO del INTERVALO, entonces es PROBABLE que el COEFICIENTE DE CORRELACION sea igual a cero y NO EXISTA UNA RELACION LINEAL IMPORTANTE ENTRE LAS VARIABLES
• Si esta FUERA del INTERVALO, entonces es MUY POCO PROBABLE que el COEFICIENTE DE CORRELACION sea igual a cero y EXISTA UNA RELACION LINEAL IMPORTANTE ENTRE LAS VARIABLES

2. p-value: es número que va entre 0 y 1:

• Cuanto más cerca de 0 este, es MUY POCO PROBABLE que el COEFICIENTE DE CORRELACION sea igual a cero
• Cuanto más cerca de 1 este, es MUY PROBABLE que el COEFICIENTE DE CORRELACION sea igual a cero

IMPORTANTE: Cuando el p-valor sea menor a 0.05 podemos estar bastante seguros que el COEFICIENTE DE CORRELACIÓN es distinto de cero y EXISTE UNA RELACIÓN LINEAL IMPORTANTE ENTRE LAS VARIABLES

Para el caso de la correlación entre Asesinatos y Analfabetismo vemos que:

1. El intervalo de confianza NO contiene al cero

2. El p-valor es mas chico que 0.05

CONCLUSIÓN: rechazamos que el COEFICIENTE DE CORRELACIÓN es distinto de cero

Veamos el test de la correlacion entre Area y Analfabetismo:

cor.test(x=estados$Area, y=estados$Illiteracy)

    Pearson's product-moment correlation

data:  estados$Area and estados$Illiteracy
t = 0.53689, df = 48, p-value = 0.5938
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.2055061  0.3481223
sample estimates:
       cor 
0.07726113 

Vemos que:

1. El intervalo de confianza SI contiene al cero

2. El p-valor es mucho mas grande que 0.05

CONCLUSIÓN: NO rechazamos que el COEFICIENTE DE CORRELACIÓN es distinto de cero

Correlación vs causalidad

La correlación mide el grado de asociación lineal entre dos variables.

La CAUSALIDAD indica que una o más variables (variables predictoras/exogenas/explicativas) sirven para explicar el comportamiento de otra variable (variable a predecir/endogena/ a explicar)

Puede haber variables con una baja (o nula) correlación, y sin embargo, que exista otro tipo de relación entre ellas. Como tambien pueden existir variables con una alta correlación pero que no podamos explicar el vinculo entre ellas.

Para este segundo caso, podemos clasificarlo en dos situaciones:

1. Correlación espuria: dos variables tienen una correlación elevada por motivos puramente azarosos

2. Variables ocultas: dos variables tienen una correlación elevada porque hay información “escondida” en alguna de ellas (o en ambas)

Veamos algunos ejemplos para aclarar la situacion

Fuente: http://www.tylervigen.com/spurious-correlations

Fuente: http://www.tylervigen.com/spurious-correlations

Ejemplo de estados con Frost vs Illiteracy

ggplot(estados, aes(x=Frost,y=Illiteracy)) + geom_point(size=2) +
  geom_hline(yintercept = mean(estados$Illiteracy),color='steelblue', linetype='dashed', size=1) +
  geom_vline(xintercept = mean(estados$Frost), color='steelblue', linetype='dashed', size=1) +
  labs(title='¿Correlacion espuria?')+
  theme_bw() + scale_y_continuous(limits = c(-1,3)) + scale_x_continuous(limits = c(-10,200))

cor(x=estados$Frost, y=estados$Illiteracy)

[1] -0.671947

Existe una correlación negativa fuerte entre Analfabetismo y Cantidad de Dias debajo de los 0 grados.

No tiene sentido decir: Para que haga menos frío, tenemos que reducir el analfabetismo en la población. Sin embargo, pueden existir factores que expliquen que los estados que más días fríos tienen, también tienen las menores tasas de analfabetismo. Por ejemplo, puede ser que en los estados donde hace más frío hay más escuelas porque la gente no viaja tanto por la nieve y eso lleva a que haya una menor tasa de analfabetismo.

Poder encontrar estas variables ocultas es una parte muy importante del modelado estadístico.

Modelo

Realizar un modelo implica contar con una variable que deseamos explicar (variable a predecir) y una o más variables que sirven para explicarla (variable/s predictora/s). A su vez, necesitamos SUPONER cual es el tipo de relación entre las variables.

Modelo lineal simple

Definicion

El modelo que vamos a usar es:

\(Y= \beta_0 + \beta_1 \cdot X + \varepsilon\)

1. Y es la variable a predecir

2. X es la variable predictora

3. La relación que suponemos entre ambas variables es una recta, por lo tanto la relación es LINEAL

4. \(\varepsilon\) es un termino de error. La relación entre las variables nunca va a ser exactamente igual a una recta, el termino de error da cuenta de cuán lejos esta nuestro modelo (nuestra recta) de los datos que vemos.

Tarea de ejemplo

Supongamos que nos interesa modelar la relación entre asesinatos y la esperanza de vida en los estados de Estados Unidos.

ggplot(estados, aes(x=Murder,y=life_exp)) +
  geom_point(size=2) +
  labs(title='Asesinatos y Esperanza de Vida', x='Asesinatos (Murder)', y='Life Exp (Experanza Vida)')+
  theme_bw()

Entonces podemos usar el modelo lineal simple para esta tarea:

1. Y sera Esperanza de Vida

2. X sera Asesinatos

Nuestro modelo nos queda igual a:

\(EsperanzaVida= \beta_0 + \beta_1 \cdot Asesinatos + \varepsilon\)

Nosotros ya contamos con los datos de la esperanza de vida y los asesinatos. Entonces la tarea que tenemos por delante es encontrar los valores de \(\beta_0\) y \(\beta_1\).

Rectas

El valor de \(\beta_0\) es la ordenada al origen de mi recta mientras que \(\beta_1\) es la pendiente

Consideremos la recta \(Y= 6 - 2X\)

# Creo la funcion
funcion <- function(x) {6 - 2*x}
# Grafico
ggplot(data.frame(x = c(-1, 3)), aes(x)) +
  stat_function(fun=funcion, colour='forestgreen', size=2)+
  geom_vline(xintercept = 0) +
  geom_hline(yintercept = 0) +
  theme_bw()

Vemos que:

Cuando \(x=0\) la variable Y es igual a la ordenada al origen. La ordenada al origen indica cual es el valor de la variable Y cuando la variable X es igual a cero

Cuando \(x=1\), \(y= 6-2.1=4\) . Cuando X aumento en 1 unidad Y se redujo en 2 unidades. La pendiente nos indica justamente eso: cuanto cambia la variable Y, cuando la variable X cambia en una unidad.

Parámetros del modelo

Volviendo a nuestro modelo de Asesinatos y Esperanza de Vida, surge la pregunta:

¿Hay valores de la ordenada al origen y la pendiente que sean mejores para explicar la relación entre ellas?

Probemos algunas rectas y veamos cual se ajusta mejor a nuestros datos

recta1 <- ggplot(estados, aes(x=Murder,y=life_exp)) +
  geom_point(size=2) +
  geom_abline(intercept = 73, slope = -0.3, colour='forestgreen', size=1.5) +
  labs(title='Recta 1', x='Asesinatos (Murder)', y='Life Exp (Experanza Vida)')+
  theme_bw()
recta2 <- ggplot(estados, aes(x=Murder,y=life_exp)) +
  geom_point(size=2) +
  geom_abline(intercept = 69, slope = 0.2, colour='firebrick', size=1.5) +
  labs(title='Recta 2', x='Asesinatos (Murder)', y='Life Exp (Experanza Vida)')+
  theme_bw()
recta3 <- ggplot(estados, aes(x=Murder,y=life_exp)) +
  geom_point(size=2) +
  geom_abline(intercept = 71, slope = -0.3, colour='steelblue', size=1.5) +
  labs(title='Recta 3', x='Asesinatos (Murder)', y='Life Exp (Experanza Vida)')+
  theme_bw()
plot_grid(recta1,recta2,recta3)

• La recta 2 atraviesa los puntos por la mitad. No se ajusta para nada a los datos

• La recta 3 deja todos los puntos por encima de ella. Sería un mal modelo

• La recta 1 parece ser que se ajusta mucho mejor a la estructura de los datos

Entonces, ya sabemos que puede haber valores de \(\beta_0\) (ordenada al origen) y \(\beta_1\) (pendiente) que son mejores que otros.

Cuando creemos un modelo lineal, la regresión va a obtener los MEJORES valores \(\beta_0\) y \(\beta_1\). Son los MEJORES porque son los que producen la RECTA que MEJOR se AJUSTA a los DATOS.

Modelo en R

La función para crear un modelo lineal en R es lm(). En ella le tenemos que pasar:

• formula: es la formula de nuestro modelo y se escribe Y~X

• data: es la tabla/dataset que estamos usando

#Creo el modelo con life_exp (variable a predecir) ~ Murder (variable predictora)
modelo_asesinatos<-lm(formula = life_exp~Murder, data = estados)
modelo_asesinatos

Call:
lm(formula = life_exp ~ Murder, data = estados)

Coefficients:
(Intercept)       Murder  
    72.9736      -0.2839  

Cuando creamos un modelo vamos a poder acceder a mucha informacion contenida en él. Podemos ver que cosas tiene dentro con modelo_asesinatos$

Interpretación

Recordemos que nuestro modelo era:

\(EsperanzaVida= \beta_0 + \beta_1 \cdot Asesinatos + \varepsilon\)

Y queríamos encontrar los valores de \(\beta_0\) y \(\beta_1\).

Si realizamos modelo_asesinatos$coefficients vamos a poder acceder a los valores de \(\beta_0\) y \(\beta_1\) estimados por la regresión simple:

modelo_asesinatos$coefficients

(Intercept)      Murder 
 72.9735623  -0.2839472 

beta_0 <- modelo_asesinatos$coefficients[[1]] # Guardamos el valor del intercepto
beta_1 <- modelo_asesinatos$coefficients[[2]] # Guardamos el valor de la pendiente

¿Cómo interpretamos los coeficientes?

• El intercepto nos indica que si no hay asesinatos, la esperanza de vida en promedio para los estados es de 73 años

• La pendiente nos indica que cada asesinato adicional reduce la esperanza de vida en promedio en 0.3 años. Es decir, que más o menos, cada 3 homicidios la esperanza de vida promedio se reduce en 1 año.

Veamos como luce nuestro modelo lineal gráficamente:

ggplot(estados, aes(x=Murder,y=life_exp)) +
  geom_point(size=2) +
  geom_abline(intercept = beta_0, slope = beta_1, colour='forestgreen', size=1.5, alpha=0.6) +
  geom_vline(xintercept = 0) +
  labs(title='Modelo lineal', x='Asesinatos (Murder)', y='Life Exp (Experanza Vida)')+
  theme_bw()

Evaluación del modelo

Hasta ahora no sabemos que tan bueno es nuestro modelo para explicar la esperanza de vida.

Para evaluar que tan bueno es nuestro modelo vamos a ver dos cosas:

1. El p-valor de la ordenada al origen y la pendiente

2. El valor del R cuadrado

p-valor

Al igual que con el coeficiente de correlación, las estimaciones de \(\beta_0\) y \(\beta_1\) van a tener un p-valor asociado.

En este caso:

• Si el p-valor es inferior a 0.05, estamos bastante seguros que la variable es buena para explicar a nuestra variable Y

• Si el p-valor es superior a 0.05, estamos bastante seguros que la variable no sirve para explicar a nuestra variable Y

R cuadrado

Es un valor que va entre 0 y 1, e indica que tan bueno es mi modelo en su totalidad para explicar mi variable Y.

IMPORTANTE: podemos usar el R cuadrado para comparar modelos distintos entre sí

• Si se encuentra muy cerca de 0, el modelo es malo para explicar a la variable Y

• Si se encuentra muy cerca de 1, el modelo es bueno para explicar a la variable Y

Con la función summary accedemos a un resumen del modelo que va a tener todas las cosas mencionadas anteriormente

# Resumen del modelo
summary(modelo_asesinatos)

Call:
lm(formula = life_exp ~ Murder, data = estados)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.81690 -0.48139  0.09591  0.39769  2.38691 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 72.97356    0.26997  270.30  < 2e-16 ***
Murder      -0.28395    0.03279   -8.66 2.26e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8473 on 48 degrees of freedom
Multiple R-squared:  0.6097,    Adjusted R-squared:  0.6016 
F-statistic: 74.99 on 1 and 48 DF,  p-value: 2.26e-11

Notemos que:

• Ambos p-valores son mas chicos que 0.05. Entonces la variable Asesinatos es buena para explicar la Esperanza de Vida

• El R cuadrado (Multiple R-squared) es igual a 0.61. Por lo que podemos decir que es un buen modelo

Ahora creemos otro modelo, para ver si su R-cuadrado es mayor o menor.

Nuestro nuevo modelo va a tratar de explicar la Esperanza de Vida a partir del Ingreso per capita

#Creo el modelo con life_exp (variable a predecir) ~ Income (variable predictora)
modelo_ingreso<-lm(formula = life_exp~Income, data = estados)
# Resumen del modelo
summary(modelo_ingreso)

Call:
lm(formula = life_exp ~ Income, data = estados)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.96547 -0.76381 -0.03428  0.92876  2.32951 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 6.758e+01  1.328e+00  50.906   <2e-16 ***
Income      7.433e-04  2.965e-04   2.507   0.0156 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.275 on 48 degrees of freedom
Multiple R-squared:  0.1158,    Adjusted R-squared:  0.09735 
F-statistic: 6.285 on 1 and 48 DF,  p-value: 0.01562

Notemos que:

• Ambos p-valores son mas chicos que 0.05. Entonces la variable Ingreso es buena para explicar la Esperanza de Vida

• El R cuadrado (Multiple R-squared) es igual a 0.12

En conclusión, nuestro primer modelo tiene un R cuadrado mucho más alto que el segundo modelo. Por lo tanto, el primer modelo es mucho mejor para predecir la esperanza de vida de los estados de Estados Unidos que el segundo.

Ejercicios de tarea

1. Realizar el gráfico de dispersión entre Income y Illiteracy, y entre Income y hs_grad

2. Calcular la covarianza y correlación entre Income y Illiteracy, y entre Income y hs_grad

3. 1. Realizar un modelo lineal que explique Income en función de Illiteracy

   2. Realizar un modelo lineal que explique Income en función de hs_grad

4. Obtener los coefientes de los 2 modelos previos. Interpretarlos brevemente

5. Obtener el R-cuadrado de los 2 modelos anteriores. Elegir el que mejor explique la variable Income

Bonus

6. Graficar la recta de cada modelo sobre el diagrama de dispersión correspondiente