Regresión Lineal multiple

Librerías

setwd("~/EALMV9") # Directorio de trabajo.

library("pacman") # Importa biblioteca "pacman". Se utiliza para hacer una mejor gestión de paquetes.

p_load("base64enc", "htmltools", "mime", "xfun", "prettydoc","readr", "knitr","DT","dplyr", "ggplot2","plotly", "gganimate","gifski","scales", "readxl", "tidyverse","cluster", "factoextra","NbClust","tidyr", "hpackedbubble", "gridExtra", "psych", "GGally", "gridExtra", "lmtest", "corrplot", "car") # Paquetes necesarios para la elaboración.

Descarga de este código.

Para fines de reproducibilidad se incluye el código para su descarga.

xfun::embed_file("U2A4.Rmd")

Download U2A4.Rmd

Descarga los datos

xfun::embed_file("confirmadosEstados.csv")

Download confirmadosEstados.csv

xfun::embed_file("defuncionesEstados.csv")

Download defuncionesEstados.csv

xfun::embed_file("negativosEstados.csv")

Download negativosEstados.csv

xfun::embed_file("sospechososEstados.csv")

Download sospechososEstados.csv

xfun::embed_file("movilidadSonora.csv")

Download movilidadSonora.csv

xfun::embed_file("datosRegresion.csv")

Download datosRegresion.csv

Confirmados, sospechosos, decesos y negativos en los municipios de Sonora

Datos del CONACYT

Importación de datos

# Leer los archivos .csv 

datos_confirmados <- read.csv("confirmadosEstados.csv")
datos_decesos <- read.csv("defuncionesEstados.csv")
datos_negativos <- read.csv("negativosEstados.csv")
datos_sospechosos <- read.csv("sospechososEstados.csv")

# Definir variables

conf_son <- t(datos_confirmados[datos_confirmados$nombre == "SONORA",]) # Confirmados SONORA

dec_son <- t(datos_decesos[datos_decesos$nombre == "SONORA",]) # Decesos SONORA

neg_son <- t(datos_negativos[datos_negativos$nombre == "SONORA",]) # negativos SONORA

sos_son <- t(datos_sospechosos[datos_sospechosos$nombre == "SONORA",]) # sospechosos SONORA

# Vector de fecha (toma desde el 22 de enero de 2020, hasta el 14 de febrero de 2021).

Fecha <- seq(from = as.Date("2020-01-20"), to = as.Date("2021-05-04"), by = "day")

# Casos confirmados
vec1 <- as.vector(conf_son)
vec2 <- vec1[4:474]
num1 <- as.numeric(vec2)
Confirmados <- as.vector(num1)

# Decesos
vec1 <- as.vector(dec_son)
vec2 <- vec1[4:474]
num1 <- as.numeric(vec2)
Decesos <- as.vector(num1)

# negativos
vec1 <- as.vector(neg_son)
vec2 <- vec1[4:474]
num1 <- as.numeric(vec2)
Negativos <- as.vector(num1)

# sospechosos
vec1 <- as.vector(sos_son)
vec2 <- vec1[4:474]
num1 <- as.numeric(vec2)
Sospechosos <- as.vector(num1)

# Generación de un marco de datos ("data frame")
datos1 <- data.frame(Fecha, Confirmados, Decesos, Negativos, Sospechosos)

Graficación de datos

A continuación, se presenta una visualización de los datos de: casos confirmados, decesos, negativos y sospechosos en Sonora

Gráfica

gcov <- ggplot(data = datos1) + 
  geom_line(aes(Fecha, Confirmados, colour = "Confirmados")) +
  geom_line(aes(Fecha, Decesos, colour = "Decesos")) +
  geom_line(aes(Fecha, Negativos, colour = "Negativos")) +
  geom_line(aes(Fecha, Sospechosos, colour = "Sospechosos")) +
  xlab("Fecha") +
  ylab("COVID-19 en México") +
  labs(colour = "casos") +
  ggtitle("Casos diarios de COVID-19 en Sonora (Fuente: CONACYT)") +
  scale_y_continuous(labels = comma)
ggplotly(gcov)

Datos de enfermedades crónicas en Sonora.

Fuente de los datos Obtenida de: https://datos.gob.mx/busca/dataset/expediente-clinico-electronico-unemes-enfermedades-cronicas--2018
URL de datos CSV crudos: http://www.cenaprece.salud.gob.mx/descargas/Excel/Diagnosticos18.csv

Expediente Clínico Electrónico UNEMES (unidades de especialidades médicas) Enfermedades Crónicas 2018

El Centro Nacional de Programas Preventivos y Control de Enfermedades (CENAPRECE), es el órgano desconcentrado de la Secretaría de Salud responsable de conducir e implementar los programas sustantivos para la prevención y control de enfermedades, para reducir la morbilidad y mortalidad en la población mexicana.

Importar datos de enfermedades crónicas en Sonora

datos2 <- read.csv("Diagnosticos18.csv", encoding = "latin1")
class(datos2)

## [1] "data.frame"

head(datos2)

##       Estado Jurisdicción             Uneme        CLUES Cve.Persona Genero
## 1 Guanajuato     IRAPUATO UNEME EC IRAPUATO GTSSA017250   GuIRUN5910  Mujer
## 2 Guanajuato     IRAPUATO UNEME EC IRAPUATO GTSSA017250   GuIRUN5953  Mujer
## 3 Guanajuato     IRAPUATO UNEME EC IRAPUATO GTSSA017250   GuIRUN6045  Mujer
## 4 Guanajuato     IRAPUATO UNEME EC IRAPUATO GTSSA017250   GuIRUN6091  Mujer
## 5 Guanajuato     IRAPUATO UNEME EC IRAPUATO GTSSA017250   GuIRUN6138  Mujer
## 6 Guanajuato     IRAPUATO UNEME EC IRAPUATO GTSSA017250   GuIRUN5706 Hombre
##   Cve.Diagnóstico   Diagnóstico Fecha.Diagnótico
## 1           E78.2 Dislipidemias       02/01/2018
## 2           E78.2 Dislipidemias       02/01/2018
## 3           E78.2 Dislipidemias       02/01/2018
## 4           E78.2 Dislipidemias       02/01/2018
## 5           E78.2 Dislipidemias       02/01/2018
## 6           E78.2 Dislipidemias       02/01/2018

Formateo de datos

SonoraS <- t(datos2[datos2$Estado == "Sonora" ,])
SonoraS <- (datos2[datos2$Estado == "Sonora" ,])

Graficación

# Gráfico agrupado
  ggplot(SonoraS, aes(fill = Jurisdicción, y = Diagnóstico)) +
  geom_bar(position = "dodge", stat = "count") +
  xlab("Número de casos") +
  ylab("Diagnóstico") +
  ggtitle("Enfermedades crónicas en Sonora, 2018 (CENAPRECE)")

Tabla interactiva de datos de enfermedades crónicas

datatable(datos2)

## Warning in instance$preRenderHook(instance): It seems your data is too big
## for client-side DataTables. You may consider server-side processing: https://
## rstudio.github.io/DT/server.html

Movilidada en Sonora

Cambio en la linea base de movilidad en el Estado de Sonora

Datos obtenidos por Google

Importación de datos

movilidad <- read.csv("movilidadSonora.csv")

FechaMR <- seq(from=as.Date("2020-02-15"), to= as.Date("2020-12-31"), by="day")


recreacion <- movilidad$retail_and_recreation_percent_change_from_baseline

farmacias_abarrotes <- movilidad$grocery_and_pharmacy_percent_change_from_baseline

parques <- movilidad$parks_percent_change_from_baseline

estaciones <- movilidad$transit_stations_percent_change_from_baseline

trabajos<- movilidad$workplaces_percent_change_from_baseline

hogares <- movilidad$residential_percent_change_from_baseline


Porcentaje_VR <- movilidad$retail_and_recreation_percent_change_from_baseline

Porcentaje_F <-  movilidad$grocery_and_pharmacy_percent_change_from_baseline

Porcentaje_P <- movilidad$parks_percent_change_from_baseline

Porcentaje_ET <- movilidad$transit_stations_percent_change_from_baseline

Porcentaje_AT <- movilidad$workplaces_percent_change_from_baseline

Porcentaje_H <- movilidad$residential_percent_change_from_baseline


#Generacion de un data frame 
datosMRRR <- data.frame(recreacion, farmacias_abarrotes, parques, estaciones, trabajos, hogares)

# Create data MR
dataMR <- data.frame(x=FechaMR,y=recreacion)

Gráfica movilidad

gMR1<- ggplot(data=dataMR) +
  geom_line(aes(x=FechaMR, y=Porcentaje_VR), size=1, colour="blue") +
  geom_hline(yintercept = 0) +
  theme_light() +
  xlab("Fecha") +
  ylab("Cambio línea base") +
  ggtitle("Comercio y recreación")

gMR2 <- ggplot(data=dataMR) +
  geom_line(aes(x=FechaMR, y=Porcentaje_F), size=1, colour="green") +
  geom_hline(yintercept = 0) +
  theme_light() +
  xlab("Fecha") +
  ylab("Cambio línea base") +
  ggtitle("Farmacias y abarrotes")

gMR3 <- ggplot(data=dataMR) +
  geom_line(aes(x=FechaMR, y=Porcentaje_P), size=1, colour="purple") +
  geom_hline(yintercept = 0) +
  theme_light() +
  xlab("Fecha") +
  ylab("Cambio línea base") +
  ggtitle("Parques")

gMR4 <- ggplot(data=dataMR) +
  geom_line(aes(x=FechaMR, y=Porcentaje_ET), size=1, colour="brown") +
  geom_hline(yintercept = 0) +
  theme_light() +
  xlab("Fecha") +
  ylab("Cambio línea base") +
  ggtitle("Estaciones de tránsito")

gMR5 <- ggplot(data=dataMR) +
  geom_line(aes(x=FechaMR, y=Porcentaje_AT), size=1, colour="red") +
  geom_hline(yintercept = 0) +
  theme_light() +
  xlab("Fecha") +
  ylab("Cambio línea base") +
  ggtitle("Espacios de trabajo")

gMR6 <- ggplot(data=dataMR) +
  geom_line(aes(x=FechaMR, y=Porcentaje_H), size=1, colour="pink") +
  geom_hline(yintercept = 0) +
  theme_light() +
  xlab("Fecha") +
  ylab("Cambio línea base") +
  ggtitle("Hogares")
grid.arrange(gMR1,gMR2,gMR3,gMR4,gMR5,gMR6)

Regresión Lineal

Importación de datos

# Los datos representan los días comprendidos del 15/02/2020 al 31/12/2020

datos <- as.data.frame(read.csv("datosRegresion.csv"))

Analizar la relación entre variables

A continuación, se estudia la relación existente entre las variables estudiadas. Es necesario conocer su relación para poder identificar los mejores predictores para nuestro modelo.

Se calcula el coeficiente de correlación de cada par de variables.
Se relaliza una representación gráficas mediante gráficos de dispersión.

round(cor(x = datos, method = "pearson"), 3)

##             comeYrecre farmaYabarr parques estTransito trabajos hogares
## comeYrecre       1.000       0.867   0.873       0.957    0.616  -0.846
## farmaYabarr      0.867       1.000   0.684       0.885    0.523  -0.708
## parques          0.873       0.684   1.000       0.819    0.458  -0.712
## estTransito      0.957       0.885   0.819       1.000    0.542  -0.785
## trabajos         0.616       0.523   0.458       0.542    1.000  -0.896
## hogares         -0.846      -0.708  -0.712      -0.785   -0.896   1.000
## confirmados     -0.051       0.090  -0.144       0.027   -0.291   0.268
##             confirmados
## comeYrecre       -0.051
## farmaYabarr       0.090
## parques          -0.144
## estTransito       0.027
## trabajos         -0.291
## hogares           0.268
## confirmados       1.000

multi.hist(x = datos, dcol = c("blue", "red"), dlty = c("dotted", "solid"),
           main = "")

ggpairs(datos, lower = list(continuous = "smooth"),
        diag = list(continuous = "barDiag"), axisLabels = "none")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Del análisis preliminar se pueden extraer las siguientes conclusiones:

Las variables de movilidad que tienen una mayor relación lineal con los casos diarios confirmados por COVID-19 son: trabajos (r= -0.291) y hogares (r= -0.268).

Generar el modelo

Se utilizará el método mixto iniciando el modelo con todas las variables como predictores y realizando la selección de los mejores predictores con la medición Akaike (AIC).

modelo <- lm(confirmados ~ datos$comeYrecre + datos$farmaYabarr + datos$parques 
             + datos$estTransito + datos$trabajos + datos$hogares, data = datos )
summary(modelo)

## 
## Call:
## lm(formula = confirmados ~ datos$comeYrecre + datos$farmaYabarr + 
##     datos$parques + datos$estTransito + datos$trabajos + datos$hogares, 
##     data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -324.74  -69.69    1.59   75.10  430.72 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         56.692     18.842   3.009  0.00284 ** 
## datos$comeYrecre     2.740      1.927   1.422  0.15599    
## datos$farmaYabarr    3.268      1.550   2.108  0.03580 *  
## datos$parques       -3.725      1.085  -3.432  0.00068 ***
## datos$estTransito    3.768      1.415   2.662  0.00817 ** 
## datos$trabajos       1.468      1.149   1.277  0.20249    
## datos$hogares       23.880      4.871   4.902 1.52e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 125.5 on 314 degrees of freedom
## Multiple R-squared:  0.2855, Adjusted R-squared:  0.2719 
## F-statistic: 20.92 on 6 and 314 DF,  p-value: < 2.2e-16

El modelo con todas las variables introducidas como predictores tiene un R2 baja (0.2855), es capaz de explicar el 28,55% de la variabilidad observada en los casos diarios confirmados de COVID-19. El p-value del modelo es 2.2e-16 por lo que no se puede aceptar que el modelo no es por azar, al menos uno de los coeficientes parciales de regresión es distinto de 0. Muchos de ellos no son significativos, lo que es un indicativo de que podrían no contribuir al modelo.

Selección de los mejores predictores

step(object = modelo, direction = "both", trace = 1)

## Start:  AIC=3109.12
## confirmados ~ datos$comeYrecre + datos$farmaYabarr + datos$parques + 
##     datos$estTransito + datos$trabajos + datos$hogares
## 
##                     Df Sum of Sq     RSS    AIC
## - datos$trabajos     1     25679 4969038 3108.8
## <none>                           4943359 3109.1
## - datos$comeYrecre   1     31838 4975197 3109.2
## - datos$farmaYabarr  1     69978 5013337 3111.6
## - datos$estTransito  1    111565 5054924 3114.3
## - datos$parques      1    185421 5128779 3118.9
## - datos$hogares      1    378326 5321685 3130.8
## 
## Step:  AIC=3108.78
## confirmados ~ datos$comeYrecre + datos$farmaYabarr + datos$parques + 
##     datos$estTransito + datos$hogares
## 
##                     Df Sum of Sq     RSS    AIC
## - datos$comeYrecre   1     27826 4996864 3108.6
## <none>                           4969038 3108.8
## + datos$trabajos     1     25679 4943359 3109.1
## - datos$farmaYabarr  1     83187 5052226 3112.1
## - datos$estTransito  1     93517 5062555 3112.8
## - datos$parques      1    232996 5202034 3121.5
## - datos$hogares      1   1111000 6080038 3171.6
## 
## Step:  AIC=3108.57
## confirmados ~ datos$farmaYabarr + datos$parques + datos$estTransito + 
##     datos$hogares
## 
##                     Df Sum of Sq     RSS    AIC
## <none>                           4996864 3108.6
## + datos$comeYrecre   1     27826 4969038 3108.8
## + datos$trabajos     1     21667 4975197 3109.2
## - datos$farmaYabarr  1    121135 5117999 3114.3
## - datos$parques      1    220284 5217148 3120.4
## - datos$estTransito  1    287598 5284463 3124.5
## - datos$hogares      1   1267057 6263921 3179.1

## 
## Call:
## lm(formula = confirmados ~ datos$farmaYabarr + datos$parques + 
##     datos$estTransito + datos$hogares, data = datos)
## 
## Coefficients:
##       (Intercept)  datos$farmaYabarr      datos$parques  datos$estTransito  
##            56.855              4.096             -3.275              4.533  
##     datos$hogares  
##            16.843

El mejor modelo resultante del proceso de selección ha sido:

modelo <- (lm(formula = confirmados ~ datos$farmaYabarr + datos$parques + 
    datos$estTransito + datos$hogares, data = datos))
summary(modelo)

## 
## Call:
## lm(formula = confirmados ~ datos$farmaYabarr + datos$parques + 
##     datos$estTransito + datos$hogares, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -339.92  -72.91    0.77   73.01  432.15 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        56.8553    18.8023   3.024 0.002700 ** 
## datos$farmaYabarr   4.0957     1.4798   2.768 0.005977 ** 
## datos$parques      -3.2747     0.8774  -3.732 0.000225 ***
## datos$estTransito   4.5329     1.0629   4.265 2.65e-05 ***
## datos$hogares      16.8427     1.8816   8.951  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 125.7 on 316 degrees of freedom
## Multiple R-squared:  0.2778, Adjusted R-squared:  0.2687 
## F-statistic: 30.39 on 4 and 316 DF,  p-value: < 2.2e-16

Es recomendable mostrar el intervalo de confianza para cada uno de los coeficientes parciales de regresión:

confint(lm(formula = confirmados ~ datos$farmaYabarr + datos$parques + 
    datos$estTransito + datos$hogares, data = datos))

##                       2.5 %    97.5 %
## (Intercept)       19.861819 93.848816
## datos$farmaYabarr  1.184227  7.007206
## datos$parques     -5.000866 -1.548442
## datos$estTransito  2.441669  6.624150
## datos$hogares     13.140698 20.544635

Validación de condiciones para la regresión múltiple lineal

plot1 <- ggplot(data = datos, aes(farmaYabarr, modelo$residuals)) +
    geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
    theme_bw()
plot2 <- ggplot(data = datos, aes(parques, modelo$residuals)) +
    geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
    theme_bw()
plot3 <- ggplot(data = datos, aes(estTransito, modelo$residuals)) +
    geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
    theme_bw()
plot4 <- ggplot(data = datos, aes(hogares, modelo$residuals)) +
    geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
    theme_bw()
grid.arrange(plot1, plot2, plot3, plot4)

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Distribución normal:

qqnorm(modelo$residuals)
qqline(modelo$residuals)

Shapiro test:

shapiro.test(modelo$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  modelo$residuals
## W = 0.99118, p-value = 0.05186

Tanto el análisis gráfico como es test de hipótesis (Shapiro test) confirman la normalidad.

Variabilidad constante de los residuos (homocedasticidad):

Al representar los residuos frente a los valores ajustados por el modelo, los primeros se tienen que distribuir de forma aleatoria en torno a cero, manteniendo aproximadamente la misma variabilidad a lo largo del eje X. Si se observa algún patrón específico, por ejemplo forma cónica o mayor dispersión en los extremos, significa que la variabilidad es dependiente del valor ajustado y por lo tanto no hay homocedasticidad.

ggplot(data = datos, aes(modelo$fitted.values, modelo$residuals)) +
geom_point() +
geom_smooth(color = "firebrick", se = FALSE) +
geom_hline(yintercept = 0) +
theme_bw()

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Breusch-Pagan test:

bptest(modelo)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo
## BP = 63.654, df = 4, p-value = 4.943e-13

No hay evidencias de homocedasticidad.

Matriz de correlación entre predictores.

corrplot(cor(dplyr::select(datos, farmaYabarr, parques, estTransito, hogares)),
         method = "number", tl.col = "black")

Análisis de Inflación de Varianza (VIF):

vif(modelo)

## datos$farmaYabarr     datos$parques datos$estTransito     datos$hogares 
##          4.749697          3.246777          8.358464          2.724810

Los predictores muestran una correlación lineal e inflación de varianza muy alta.

Autocorrelación:

dwt(modelo, alternative = "two.sided")

##  lag Autocorrelation D-W Statistic p-value
##    1       0.8059061     0.3865268       0
##  Alternative hypothesis: rho != 0

Identificación de posibles valores atípicos o influyentes

Existe una Autocorelación alta.

datos$studentized_residual <- rstudent(modelo)
ggplot(data = datos, aes(x = predict(modelo), y = abs(studentized_residual))) +
geom_hline(yintercept = 3, color = "grey", linetype = "dashed") +
# se identifican en rojo observaciones con residuos estandarizados absolutos > 3
geom_point(aes(color = ifelse(abs(studentized_residual) > 3, 'red', 'black'))) +
scale_color_identity() +
labs(title = "Distribución de los residuos studentized",
     x = "predicción modelo") + 
theme_bw() + theme(plot.title = element_text(hjust = 0.5))

which(abs(datos$studentized_residual) > 3)

## [1] 143 150

Se identificaron dos valores atípicos.

summary(influence.measures(modelo))

## Potentially influential observations of
##   lm(formula = confirmados ~ datos$farmaYabarr + datos$parques +      datos$estTransito + datos$hogares, data = datos) :
## 
##     dfb.1_ dfb.dt$Y dfb.dts$p dfb.dt$T dfb.dts$h dffit   cov.r   cook.d hat    
## 1   -0.03  -0.02    -0.03      0.03     0.02     -0.05    1.05_*  0.00   0.03  
## 6   -0.02   0.01    -0.02      0.00     0.00     -0.02    1.06_*  0.00   0.04  
## 10  -0.06   0.02    -0.04      0.00     0.01     -0.06    1.05_*  0.00   0.03  
## 11  -0.01   0.00     0.00      0.00     0.00     -0.01    1.06_*  0.00   0.04  
## 27  -0.16   0.28     0.24     -0.37     0.04     -0.42_*  1.03    0.04   0.06_*
## 30   0.01   0.01     0.02     -0.02    -0.01      0.02    1.11_*  0.00   0.08_*
## 31  -0.17  -0.23    -0.24     -0.05    -0.44     -0.65_*  0.95_*  0.08   0.05_*
## 32  -0.12   0.01     0.12     -0.23    -0.14     -0.39_*  0.95    0.03   0.03  
## 33  -0.02   0.48     0.48     -0.70    -0.16     -0.73_*  0.97    0.11   0.07_*
## 65   0.01   0.07     0.06     -0.03     0.06     -0.12    1.05_*  0.00   0.04  
## 72   0.01   0.04     0.04     -0.01     0.04     -0.08    1.06_*  0.00   0.04  
## 86   0.00   0.00     0.00      0.00     0.00      0.00    1.06_*  0.00   0.04  
## 143 -0.02   0.05     0.08     -0.09     0.06      0.25    0.87_*  0.01   0.01  
## 144  0.02  -0.05     0.07      0.01     0.10      0.22    0.90_*  0.01   0.01  
## 145 -0.02  -0.06    -0.01      0.07     0.10      0.19    0.94_*  0.01   0.01  
## 150 -0.01   0.01     0.08     -0.07     0.07      0.27    0.84_*  0.01   0.01  
## 151  0.02  -0.01     0.08     -0.01     0.10      0.19    0.94_*  0.01   0.01  
## 309  0.00   0.01     0.08     -0.07     0.00     -0.11    1.06_*  0.00   0.04  
## 312  0.00   0.00    -0.02      0.03     0.02      0.05    1.08_*  0.00   0.06_*
## 313  0.05  -0.16     0.19     -0.17    -0.21     -0.45_*  1.11_*  0.04   0.11_*
## 314  0.16  -0.60     0.08      0.21    -0.18     -0.70_*  1.09_*  0.10   0.12_*
## 315  0.01  -0.04     0.00      0.02    -0.01      0.05    1.05_*  0.00   0.03  
## 317  0.00   0.01     0.01     -0.01     0.00     -0.02    1.06_*  0.00   0.05  
## 320 -0.01   0.03    -0.02      0.01     0.03      0.05    1.08_*  0.00   0.06_*
## 321  0.08  -0.26     0.02      0.12    -0.06     -0.28    1.16_*  0.02   0.13_*

En la tabla generada se recogen las observaciones que son significativamente influyentes en al menos uno de los predictores (una columna para cada predictor). Las tres últimas columnas son 3 medidas distintas para cuantificar la influencia. A modo de guía se pueden considerar excesivamente influyentes aquellas observaciones para las que:

Leverages (hat): Se consideran observaciones influyentes aquellas cuyos valores hat superen \(2.5((p+1)/n)\), siendo p el número de predictores y n el número de observaciones. Distancia Cook (cook.d): Se consideran influyentes valores superiores a 1. La visualización gráfica de las influencias se obtiene del siguiente modo:

influencePlot(modelo)

##        StudRes         Hat      CookD
## 33  -2.6938534 0.069038367 0.10554109
## 143  3.2049500 0.005922210 0.01188986
## 150  3.5075487 0.005698145 0.01361413
## 314 -1.8779979 0.120464554 0.09584449
## 321 -0.7243749 0.131108312 0.01585897

Conclusión

El modelo lineal múltiple \(Casos diarios confirmados de COVID-19= 4.096farmaYabarr −3.275parques+4.533estTransito−16.843hogares\)

es capaz de explicar el 28.55% de la variabilidad observada en los casos diarios confirmados de COVID-19 en Sonora. El test F muestra que es significativo (p-value: 2.2e-16).

U2A4

María José Encinas, Rafael Gutiérrez, Carlos Alvarez, Paul Becerra

6/5/2021

Regresión Lineal multiple

Librerías

Descarga de este código.

Descarga los datos

Confirmados, sospechosos, decesos y negativos en los municipios de Sonora

Importación de datos

Graficación de datos

Gráfica

Datos de enfermedades crónicas en Sonora.

Importar datos de enfermedades crónicas en Sonora

Formateo de datos

Graficación

Tabla interactiva de datos de enfermedades crónicas

Movilidada en Sonora

Importación de datos

Gráfica movilidad

Regresión Lineal

Importación de datos

Analizar la relación entre variables

Generar el modelo

Selección de los mejores predictores

Validación de condiciones para la regresión múltiple lineal

Matriz de correlación entre predictores.

Análisis de Inflación de Varianza (VIF):

Identificación de posibles valores atípicos o influyentes

Conclusión