Regreso a la regresión

Tercera sesión

Recapitulación de la regresión lineal. Para esta sesión vamos a repasar la significancia de los coeficientes de regresión, \(\beta\), que nos marcan el efecto de una de las variables independientes sobre la variable dependiente.

Para este primer ejemplo analizaremos variables asociadas con el precio de viviendas seleccionadas para poblaciones de Los Ángeles, California.

# Librerias
library(dplyr)   # Manipulación de datos y uso de %>%
library(ggplot2) # Hacer graficas
options(scipen=999)

# 1. Abrimos las bases 
# Base de precios de casas en comunidades de Los Angeles
LAhomes <- read.csv("BDs/LAhomes.csv")

# Exploramos contenido
head(LAhomes)

##         city type bed bath garage sqft pool spa  price
## 1 Long Beach        0    1         513       NA 119000
## 2 Long Beach        0    1         550       NA 153000
## 3 Long Beach        0    1         550       NA 205000
## 4 Long Beach        0    1      1 1030       NA 300000
## 5 Long Beach        0    1      1 1526       NA 375000
## 6 Long Beach        1    1         552       NA 159900

Recordemos la

Receta para hacer regresiones:

Para hacer regresiones hay que seguir los siguientes pasos:

• 1. Escribir la fórmula de la regresión.

Formato: Escribimos primero la variable dependiente, después una tilde (~) y luego las variables independientes, sumandolas con el signo +.

Ejemplo:

formula <- var_dependiente ~ var_indep_1 + var_indep_2

• 2. Utilizamos la función lm(formula, data) donde formula es la fórmula que escribimos en el paso 1, y data es la base de datos de donde se van a obtener los datos.
• 3. Le damos enter, y listo! Fácil 😉.

Y la utilizaremos para calcular regresiones…

Para nuestra primer regresión, calcularemos el efecto de la superficie de la vivienda sqft medido en pies cuadrados, sobre el precio final, price, medido en USD.

Posteriormente, graficaremos los datos observados para ver esta relación.

#################################
# 1. Precio de las casas en LA. #
#################################

# 0. Exploramos los datos: 
lapply(LAhomes, class)

## $city
## [1] "factor"
## 
## $type
## [1] "factor"
## 
## $bed
## [1] "integer"
## 
## $bath
## [1] "numeric"
## 
## $garage
## [1] "factor"
## 
## $sqft
## [1] "integer"
## 
## $pool
## [1] "factor"
## 
## $spa
## [1] "logical"
## 
## $price
## [1] "numeric"

# 1. Determinamos la formula
fmla <- price ~ sqft

# 2. Aplicamos el modelo
model <- lm(fmla, LAhomes)
print(model)

## 
## Call:
## lm(formula = fmla, data = LAhomes)
## 
## Coefficients:
## (Intercept)         sqft  
##    -1661892         1486

# 3. Checamos mayor informaci'on del modelo
summary(model)

## 
## Call:
## lm(formula = fmla, data = LAhomes)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -7732777  -363898   246242   601763 50004033 
## 
## Coefficients:
##                Estimate  Std. Error t value            Pr(>|t|)    
## (Intercept) -1661892.39    64459.91  -25.78 <0.0000000000000002 ***
## sqft            1485.99       22.71   65.44 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1859000 on 1592 degrees of freedom
## Multiple R-squared:  0.729,  Adjusted R-squared:  0.7288 
## F-statistic:  4282 on 1 and 1592 DF,  p-value: < 0.00000000000000022

# 4. Prediction 
LAhomes$pred1 <- predict(model)

# 5. Graficamos
(plot <- ggplot(LAhomes, 
       aes(x = sqft, y = price)) + 
  geom_point() + 
  geom_line(aes(x = sqft, y = pred1), colour = "blue") + 
  theme_bw()   # Fondo blanco con rayas grises
)

plotly::ggplotly(plot)

# Preguntas... que tan buena es esta regresipn para explicar el precio de la vivienda?

Ahora, para nuestra segunda regresión, utilizaremos una base de datos de salud para estimar la mortalidad por 100,000 habitantes deathrate en base a diversas variables como:

• Disponibilidad de doctores por cada 100,000 habs.

• Disponibilidad de hospitales.

• Ingreso Anual Per Cápita.

• Densidad de Población.

##########################
# 2. Tasa de Mortalidad  #
##########################

# 1. Abrimos la base
health <- readxl::read_xls("BDs/health.xls")

# 2. Exploramos los datos
head(health)

## # A tibble: 6 x 5
##   deathRate doctorAvailabil… hospitalAvailab… annualPerCapita…
##       <dbl>            <dbl>            <dbl>            <dbl>
## 1      8                  78              284             9.10
## 2      9.30               68              433             8.70
## 3      7.5                70              739             7.20
## 4      8.90               96             1792             8.90
## 5     10.2                74              477             8.30
## 6      8.30              111              362            10.9 
## # … with 1 more variable: populationDensity <dbl>

# Exploración de datos. Checar las clases de las variables de salud.
lapply(health, class)

## $deathRate
## [1] "numeric"
## 
## $doctorAvailability
## [1] "numeric"
## 
## $hospitalAvailability
## [1] "numeric"
## 
## $annualPerCapitaIncome
## [1] "numeric"
## 
## $populationDensity
## [1] "numeric"

# 3. Creamos el modelo (los modelos) declarando primero las fórmulas.
fmla1 <- deathRate ~ doctorAvailability
fmla2 <- deathRate ~ hospitalAvailability
fmla3 <- deathRate ~ annualPerCapitaIncome
fmla4 <- deathRate ~ populationDensity
fmla5 <- deathRate ~ populationDensity + annualPerCapitaIncome + hospitalAvailability + doctorAvailability

Una vez declaradas las fórmulas, pasamos a declarar los modelos. Recordemos que para la regresión lineal, utilizamos la función lm.

La función predict sirve para, una vez declarado el modelo, podamos obtener los valores que la variable independiente debería tener si se comportara exactamente como se comporta el modelo.

Recordemos, de la sesión anterior:

¿Cómo interpreto los P-valores en el análisis de regresión lineal?

El p-valor para cada término comprueba la hipótesis nula de que el coeficiente es igual a cero (no tiene efecto). Un p-valor bajo (< 0.05) indica que puedes rechazar la hipótesis nula. En otras palabras, un predictor que tenga un p-valor bajo es probable que tenga una adición significativa a su modelo porque los cambios en el valor del predictor están relacionados con cambios en la variable de respuesta.

Recíprocamente, un p-valor grande (insignificante) sugiere que los cambios en el predictor no están asociados con cambios en la respuesta. Mas info en este enlace.

Entonces, el valor nos sirve para probar la siguiente hipótesis:

\[H_0: \beta = 0 \] vs.

\[H_a: \beta \neq 0\]

¿Cómo interpretamos la R²?

\(R^2\) nos dice qué porcentaje de la variabilidad total en la variable Y puede ser explicada por la variable regresora, en consecuencia es una medida de la capacidad de PREDICCIÓN del modelo.

\(R^2\) también puede verse como es una medida de la fuerza de la ASOCIACIÓN LINEAL entre X e Y. (Hacemos énfasis en la palabra lineal porque fue obtenido bajo un modelo lineal)

# 4. Hacemos las regresiones (modelo lineal).
model1 <- lm(fmla1, health) 
summary(model1)

## 
## Call:
## lm(formula = fmla1, data = health)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.5426 -0.9798  0.2387  0.9751  3.7387 
## 
## Coefficients:
##                    Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)        8.715905   0.744665  11.704 0.000000000000000461 ***
## doctorAvailability 0.005080   0.006103   0.832                0.409    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.667 on 51 degrees of freedom
## Multiple R-squared:  0.0134, Adjusted R-squared:  -0.005944 
## F-statistic: 0.6928 on 1 and 51 DF,  p-value: 0.4091

model2 <- lm(fmla2, health) 
summary(model2)

## 
## Call:
## lm(formula = fmla2, data = health)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.5388 -0.9366  0.0217  0.9567  3.6086 
## 
## Coefficients:
##                       Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)          8.9796442  0.4699376  19.108 <0.0000000000000002 ***
## hospitalAvailability 0.0005528  0.0006956   0.795               0.431    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.668 on 51 degrees of freedom
## Multiple R-squared:  0.01223,    Adjusted R-squared:  -0.007138 
## F-statistic: 0.6315 on 1 and 51 DF,  p-value: 0.4305

model3 <- lm(fmla3, health) 
summary(model3)

## 
## Call:
## lm(formula = fmla3, data = health)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.9811 -0.8620  0.2975  1.1987  2.9530 
## 
## Coefficients:
##                       Estimate Std. Error t value    Pr(>|t|)    
## (Intercept)            11.8145     2.0249   5.835 0.000000372 ***
## annualPerCapitaIncome  -0.2659     0.2132  -1.247       0.218    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.654 on 51 degrees of freedom
## Multiple R-squared:  0.02958,    Adjusted R-squared:  0.01055 
## F-statistic: 1.555 on 1 and 51 DF,  p-value: 0.2181

model4 <- lm(fmla4, health) 
summary(model4)

## 
## Call:
## lm(formula = fmla4, data = health)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.7119 -1.0315  0.1011  1.0413  2.9696 
## 
## Coefficients:
##                    Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)       10.387997   0.569350  18.245 <0.0000000000000002 ***
## populationDensity -0.009782   0.004740  -2.064              0.0442 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.613 on 51 degrees of freedom
## Multiple R-squared:  0.07707,    Adjusted R-squared:  0.05897 
## F-statistic: 4.259 on 1 and 51 DF,  p-value: 0.04416

model5 <- lm(fmla5, health)
summary(model5)

## 
## Call:
## lm(formula = fmla5, data = health)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.6404 -0.7904  0.3053  0.9164  2.7906 
## 
## Coefficients:
##                         Estimate Std. Error t value    Pr(>|t|)    
## (Intercept)           12.2662552  2.0201467   6.072 0.000000195 ***
## populationDensity     -0.0094629  0.0048868  -1.936      0.0587 .  
## annualPerCapitaIncome -0.3302302  0.2345518  -1.408      0.1656    
## hospitalAvailability   0.0005837  0.0007219   0.809      0.4228    
## doctorAvailability     0.0073916  0.0069336   1.066      0.2917    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.601 on 48 degrees of freedom
## Multiple R-squared:  0.1437, Adjusted R-squared:  0.07235 
## F-statistic: 2.014 on 4 and 48 DF,  p-value: 0.1075

# 5. Prediccion 
health$pred1 <- predict(model1)
health$pred2 <- predict(model2)
health$pred3 <- predict(model3)
health$pred4 <- predict(model4)
health$pred5 <- predict(model5)

# 6. BONUS. Cálculo del error!
health$error1 <- health$deathRate - health$pred1

Ahora, pasamos a dibujar las gráficas de puntos y líneas que describen el fenómeno y el modelo predictivo.

# 6. Grafica

# Sobre popDensity, la mas significativa
ggplot(health, 
       aes(x = populationDensity, y = deathRate)) + 
  geom_point() + 
  geom_line(aes(x = populationDensity, y = pred4), colour = "blue") + 
  theme_bw()   # Fondo blanco con rayas grises

# Sobre DoctorAvailability
ggplot(health, 
       aes(x = doctorAvailability, y = deathRate)) + 
  geom_point() + 
  geom_line(aes(x = doctorAvailability, y = pred1), colour = "blue") + 
  theme_bw()   # Fondo blanco con rayas grises

# Como podemos ver, en este caso el modelo casi no predice nada. 

# Y es por esto que esta variable no es significante.