Importación de datos

Se leen los archivos con los que se estará trabajando en el área de COVID-19. Estos datos representan a los confirmados, muertes, sospechosos y negativos del Estado de Sonora.

# 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")
confimados_mun <- read.csv("confirmadosMunicipios.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

Tabla de confirmados en los municipios de Sonora

datatable(confimados_mun)

Se crea un vector fecha para visualizar la serie de tiempo que se abarca en la gráfica con las variables que se muestran.

# 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

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: COECYT)") +
  scale_y_continuous(labels = comma)
ggplotly(gcov)

Se aprecian las variaciones respecto a cada variable observada. Los datos respecto a los casos confirmados de COVID-19 presentan una alza después de los periodos vacacionales de verano e invierno. Podemos inferir que esto se debe a que las personas salieron de vacaciones, se reunieron familias, y esto provocó un aumento en contagios.

En el caso de los decesos, se puede decir que los aumentos concuerdan con el periodo de aumento de los casos confirmados.

Con los valores en los casos negativos se puede inferir que, conforme pasaba el tiempo más personas se realizaban pruebas y estas salían negativas.

Al inicio de la pandemia se puede apreciar más casos sospechosos. Esto puede deberse a la preocupación de la población ante el contagio de un nuevo virus.

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

Se leen los datos sobre las enfermedades.

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

## [1] "data.frame"

head(datos2)

##     Estado Jurisdiccion             Uneme       CLUES Cve.Persona Genero
## 1 Guerrero     ACAPULCO UNEME EC ACAPULCO GRSSA012354  GuACUN2751  Mujer
## 2 Guerrero     ACAPULCO UNEME EC ACAPULCO GRSSA012354  GuACUN5442  Mujer
## 3 Guerrero     ACAPULCO UNEME EC ACAPULCO GRSSA012354  GuACUN5995  Mujer
## 4 Guerrero     ACAPULCO UNEME EC ACAPULCO GRSSA012354  GuACUN6110 Hombre
## 5 Guerrero     ACAPULCO UNEME EC ACAPULCO GRSSA012354  GuACUN6113 Hombre
## 6 Guerrero     ACAPULCO UNEME EC ACAPULCO GRSSA012354  GuACUN5995  Mujer
##   Cve.Diagnóstico   Diagnostico Fecha.Diagnótico
## 1           E11.9      Diabetes       08/01/2018
## 2           E11.9      Diabetes       08/01/2018
## 3           E11.9      Diabetes       08/01/2018
## 4           E11.9      Diabetes       08/01/2018
## 5           E11.9      Diabetes       08/01/2018
## 6           E78.6 Dislipidemias       08/01/2018

Formateo de datos

Se obtienen los datos del Estado de Sonora, que es el lugar que se desea analizar. Se transversan los datos obtenerlos de una forma más entendible.

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

Graficación

A continuación, se presenta una visualización de los datos sobre las enfermedades crónicas en Sonora en las jurisdicciones de Caborca, Ciudad Obregón, Hermosillo, Santa Ana

# Gráfico agrupado
ggplot(SonoraS, aes(fill = Jurisdiccion, y = Diagnostico)) +
  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

Importación de datos

Se leen los datos sobre movilidad, se toma cada grupo y se aplica una serie de tiempo para su visualización.

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

Se presenta de manera gráfica los cambios de movilidad que ha habido en las siguientes áreas: comercio y recreación, farmacias y abarrotes. parques, estaciones de tránsito, espacios de trabajo y hogares.

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)

Comercio y recreación

En esta gráfica se puede apreciar que el cambio sobre la línea base (la cual es cero). Podemos observar que a partir del 13 de marzo, empezando la cuarentena, la movilidad empezó a bajar, pero a partir del 10 de abril empieza a fluctuar y va subiendo poco a poco. Se puede apreciar que en el mes de diciembre empieza a subir y esto podría ser debido a las compras navideñas y el año nuevo, dado que el 25 de diciembre y el 1 de enero tiene una movilidad muy baja.

Farmacias y abarrotes

En esta gráfica los valores que más se destacan son el 24 de diciembre y el 31 de enero por ser muy altos. Se puede deducir que estos valores se deben a las preparaciones para las cenas de navidad y año nuevo respectivamente, dado que la movilidad más baja son el 25 de diciembre y el 1 de enero.

Parques

La movilidad en los parques de Sonora se ha mantenido por debajo de la línea base durante la mayor parte de nuestra gráfica.

Estaciones de tránsito

Podemos apreciar que hubo una gran disminución de la movilidad a partir del 17 de marzo, esto ya que con la cuarentena se frenaron muchas actividades como escuelas y trabajos presenciales. Los valores comenzarón con una lenta subida, alcanzando sus valores máximos justamente el 23 de diciembre, justamente 33 puntos sobre la línea base. Justamente un año desde que inicio la cuarentena la movilidad se pisicionó en -31 puntos sobre la línea base.

Espacios de trabajo

Podemos aperciar que la movilidad ha fluctuado demasiado desde el inicio de nuestra gráfica el 7 de marzo de 2020 hasta un año después, pero generalmente se mantiene por debajo de la línea base y esta bajó más en los días festivos.

Hogares

En esta gráfica las personas que permanecen en case ha subido más de la línea base inicial, dado a que muchas actividades se han suspendido o reducido por la pandemia. Podemos apreciar mayores valores de en los primeros meses, pero se ha mantenido fluctuando.

Importación de datos

Se leen los datos sobre el cambio en la línea base en las diferentes áreas de movilidad, además de los contagios ocurridos en el mismo lapso de tiempo (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

La siguiente sentencia hace una selección de las variables que más influyen en los contagios por COVID-19. En este caso se van a emplear la estrategia de stepwise mixto. El valor matemático empleado para determinar la calidad del modelo va a ser Akaike(AIC).

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

Cada una de las pendientes de un modelo de regresión lineal múltiple (coeficientes parciales de regresión de los predictores) se define del siguiente modo: Si el resto de variables se mantienen constantes, por cada unidad que aumenta el predictor en cuestión, la variable (Y) varía en promedio tantas unidades como indica la pendiente.

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

• Relación lineal entre los predictores numéricos y la variable de respuesta (nuestra movilidad y nuestro número de casos confirmados):

Esta condición se puede validar bien mediante diagramas de dispersión entre la variable dependiente y cada uno de los predictores o con diagramas de dispersión entre cada uno de los predictores y los residuos del modelo. Si la relación es lineal, los residuos deben de distribuirse aleatoriamente en torno a 0 con una variabilidad constante a lo largo del eje X. Esta última opción suele ser más indicada ya que permite identificar posibles datos atípicos.

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 el 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

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).

importación de datos

indicadores <- read.csv("indicadoresSectores.csv")

Graficación

ff <- data.frame(indicadores)

f <- ff %>% 
  ggplot(aes(x=poblacion, y=personal.ocupado.total, size=produccion.bruta.total, colour=municipio)) +
  geom_point(alpha=0.8)+
geom_text(aes(label = municipio), hjust = 0, vjust = 0)+
  scale_size_continuous(range = c(1, 9))+
  labs(title="Carrera a la 'Reactivación económica'", x = "Población" , y = "Personal ocupado", colour= "Municipio", size="Producción bruta")+
  scale_y_continuous(labels = comma)+
  scale_x_continuous(labels = comma)
ggplotly(f)

Podemos apreciar que Hermosillo es mucho más grande en población, personal ocupado y en producción bruta, lo cuál tiene sentido dado que es la cápital del Estado. Está gráfica nos dá un indicativo del lugar en el que se encuntran los municipios.