Análisis espacio-temporal

En el notebook anterior hiucimos una análisis exploratorio básico de la estructura espacial pandemia en México. Ahora vamos a analizar con algo más de detalle la evolución espacio-temporal. Para hacerlo vamos a emplear técnicas que vienen del análisis regional, en general son técnicas para estudiar la convergencia en el desarrollo económico entre regiones. Sin embargo, aunque no estén nrelacionadas diréctamente con la epidemiología, las técnicas que vamos a usar resultan ilustrativas de cómo la epidemia evoluciona en el espacio.

Para correr el notebook en colab hay que instalar los paquetes:

# Important library for many geopython libraries
!apt install gdal-bin python-gdal python3-gdal 
# Install rtree - Geopandas requirment
!apt install python3-rtree 
# Install Geopandas
!pip install git+git://github.com/geopandas/geopandas.git
# Install descartes - Geopandas requirment
!pip install descartes 
# Install Folium for Geographic data visualization
!pip install folium
# Install plotlyExpress
!pip install mapclassify
!pip install palettable
!pip install libpysal
!pip install esda
!pip install giddy
!pip install splot
!pip install pickle5

También necesitamos bajar los datos:

url = "https://www.dropbox.com/s/kf9dldnqgo4eidu/agregados_semana_municipio.pkl?dl=1"
r = requests.get(url, allow_redirects=True)
open('/content/semana_municipio.pkl', 'wb').write(r.content)

url = "https://www.dropbox.com/s/46h5pnrhgnw9qca/pca_vulnerability.csv?dl=1"
r = requests.get(url, allow_redirects=True)
open('/content/pca_vulnerability.csv', 'wb').write(r.content)

url = "https://www.dropbox.com/s/md8gl5oy3mpequg/municipios.gpkg?dl=1"
r = requests.get(url, allow_redirects=True)
open('/content/municipios.gpkg', 'wb').write(r.content)

%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from libpysal.weights.contiguity import Queen
import giddy
import numpy as np
import pandas as pd
import geopandas as gpd
import splot
from splot.esda import moran_scatterplot
from esda.moran import Moran, Moran_Local, Moran_Local_BV
from splot.esda import lisa_cluster
import seaborn as sns
import contextily as ctx
import pickle5 as pickle
plt.style.use('ggplot')

Vamos a realizar el análisis para la tasa de casos confirmados, porque ya vimos que tiene mayor estructura espacial. Empezamos por leer los datos y preprocesarlos (unir a población, sacar tasas)

with open("data/semana_municipio.pkl", "rb") as fh:
  confirmados_por_semana = pickle.load(fh)

confirmados_por_semana = confirmados_por_semana.groupby(['CLAVE_MUNICIPIO_RES','FECHA_SINTOMAS'])[['Nuevos Casos']].sum()
confirmados_por_semana_wide = (confirmados_por_semana
                   .reset_index()
                   .pivot_table("Nuevos Casos", "CLAVE_MUNICIPIO_RES", "FECHA_SINTOMAS")
                  )
muns_geo = gpd.read_file("data/municipios.gpkg")
confirmados_por_semana_wide = (confirmados_por_semana_wide
                   .merge(muns_geo, left_on='CLAVE_MUNICIPIO_RES', right_on='municipio_cvegeo', how='right')
                   .fillna(0))
confirmados_por_semana_wide.head()

	2020-01-05 00:00:00	2020-01-12 00:00:00	2020-01-19 00:00:00	2020-01-26 00:00:00	2020-02-02 00:00:00	2020-02-09 00:00:00	2020-02-16 00:00:00	2020-02-23 00:00:00	2020-03-01 00:00:00	2020-03-08 00:00:00	...	2022-01-02 00:00:00	2022-01-09 00:00:00	2022-01-16 00:00:00	2022-01-23 00:00:00	_uid_	entidad_cvegeo	municipio_cvegeo	municipio_nombre	pob	geometry
0	1.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	...	9.0	2.0	2.0	0.0	1	16	16046	Juárez	15290.0	MULTIPOLYGON (((-100.45897 19.32952, -100.4588...
1	0.0	1.0	0.0	0.0	0.0	0.0	0.0	1.0	0.0	1.0	...	3.0	1.0	3.0	0.0	2	16	16047	Jungapeo	22358.0	MULTIPOLYGON (((-100.55704 19.51833, -100.5561...
2	0.0	0.0	0.0	0.0	0.0	0.0	0.0	1.0	0.0	0.0	...	0.0	1.0	2.0	0.0	3	16	16048	Lagunillas	5862.0	MULTIPOLYGON (((-101.38110 19.60250, -101.3811...
3	0.0	1.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	...	1.0	0.0	1.0	0.0	4	16	16049	Madero	18769.0	MULTIPOLYGON (((-101.11654 19.53720, -101.1168...
4	3.0	17.0	34.0	59.0	52.0	68.0	62.0	40.0	42.0	28.0	...	1431.0	2768.0	1759.0	13.0	5	05	05035	Torreón	744247.0	MULTIPOLYGON (((-103.42797 25.30109, -103.3861...

5 rows × 114 columns

Markov espacial

El primer análisis que vamos a hacer es ver la evolución de la epidemia como una cadena de Markov en la que el estado de cada municipio en un tiempo determinado es función de su estado en el tiempo previo (en este caso los momentos son semanas epidemiológicas) y el estado de los municipios vecinos.

Las cadenas de Markov funcionan sobre estados discretos, mientras que nosotros tenemos datos continuos, entonces antes tenemos que cuantizarlos de alguna forma. Un problema que se presenta aquí es que por la forma en la que la epidemia va creciendo, al principio hay poca variabilidad geográfica en el número (o la tasa) de casos, la mayoría son ceros. Entonces, para estudiar la epidemia con esta técnica, tenemos que empezar a partir de algún momento que nos permita discretizar la tasa en algún número razonable de intervalos, digamos, por lo menos 4. Esto implica que no podemos realizar el análisis para toda la duración de la epidemia, solo de cierto momento en adelante.

Una vez que seleccionemos el número de clases (y por lo tanto el momento en el que vamos a iniciar el análisis), Markov espacial nos permite calcular la dinámica de las transiciones de clases condicionadas por el contexto local. La dinámica de las transiciones es dividida a lo largo de las diferentes clases a las que pertenezcan los vecinos. De esta forma, tendremos un régimen temporal para municipios cuyos vecinos estén en la clase más baja de casos confirmados, otra para los que tengan vecinos en la segunda clase y así sucesivamente.

Si pensamos que queremos 4 clases, entonces es a partir de la semana 12 que podemos segmentar los datos en cuartiles.

Primero, calculamos la matriz de pesos espaciales

confirmados_por_semana_wide = confirmados_por_semana_wide.drop([1986], axis=0) # quitamos Cozumel
w = Queen.from_dataframe(confirmados_por_semana_wide)
w.T = 'r' # estandarizamos por fila

/home/plablo/miniconda3/envs/geoinformatica/lib/python3.7/site-packages/libpysal/weights/weights.py:172: UserWarning: The weights matrix is not fully connected: 
 There are 2 disconnected components.
  warnings.warn(message)

Obtenemos la matriz de datos (observaciones en filas, tiempo en columnas), calculamos las tasas y obtenemos el modelo de Markov espacial utilizando la librería giddy

list(confirmados_por_semana_wide.columns)

[Timestamp('2020-01-05 00:00:00'),
 Timestamp('2020-01-12 00:00:00'),
 Timestamp('2020-01-19 00:00:00'),
 Timestamp('2020-01-26 00:00:00'),
 Timestamp('2020-02-02 00:00:00'),
 Timestamp('2020-02-09 00:00:00'),
 Timestamp('2020-02-16 00:00:00'),
 Timestamp('2020-02-23 00:00:00'),
 Timestamp('2020-03-01 00:00:00'),
 Timestamp('2020-03-08 00:00:00'),
 Timestamp('2020-03-15 00:00:00'),
 Timestamp('2020-03-22 00:00:00'),
 Timestamp('2020-03-29 00:00:00'),
 Timestamp('2020-04-05 00:00:00'),
 Timestamp('2020-04-12 00:00:00'),
 Timestamp('2020-04-19 00:00:00'),
 Timestamp('2020-04-26 00:00:00'),
 Timestamp('2020-05-03 00:00:00'),
 Timestamp('2020-05-10 00:00:00'),
 Timestamp('2020-05-17 00:00:00'),
 Timestamp('2020-05-24 00:00:00'),
 Timestamp('2020-05-31 00:00:00'),
 Timestamp('2020-06-07 00:00:00'),
 Timestamp('2020-06-14 00:00:00'),
 Timestamp('2020-06-21 00:00:00'),
 Timestamp('2020-06-28 00:00:00'),
 Timestamp('2020-07-05 00:00:00'),
 Timestamp('2020-07-12 00:00:00'),
 Timestamp('2020-07-19 00:00:00'),
 Timestamp('2020-07-26 00:00:00'),
 Timestamp('2020-08-02 00:00:00'),
 Timestamp('2020-08-09 00:00:00'),
 Timestamp('2020-08-16 00:00:00'),
 Timestamp('2020-08-23 00:00:00'),
 Timestamp('2020-08-30 00:00:00'),
 Timestamp('2020-09-06 00:00:00'),
 Timestamp('2020-09-13 00:00:00'),
 Timestamp('2020-09-20 00:00:00'),
 Timestamp('2020-09-27 00:00:00'),
 Timestamp('2020-10-04 00:00:00'),
 Timestamp('2020-10-11 00:00:00'),
 Timestamp('2020-10-18 00:00:00'),
 Timestamp('2020-10-25 00:00:00'),
 Timestamp('2020-11-01 00:00:00'),
 Timestamp('2020-11-08 00:00:00'),
 Timestamp('2020-11-15 00:00:00'),
 Timestamp('2020-11-22 00:00:00'),
 Timestamp('2020-11-29 00:00:00'),
 Timestamp('2020-12-06 00:00:00'),
 Timestamp('2020-12-13 00:00:00'),
 Timestamp('2020-12-20 00:00:00'),
 Timestamp('2020-12-27 00:00:00'),
 Timestamp('2021-01-03 00:00:00'),
 Timestamp('2021-01-10 00:00:00'),
 Timestamp('2021-01-17 00:00:00'),
 Timestamp('2021-01-24 00:00:00'),
 Timestamp('2021-01-31 00:00:00'),
 Timestamp('2021-02-07 00:00:00'),
 Timestamp('2021-02-14 00:00:00'),
 Timestamp('2021-02-21 00:00:00'),
 Timestamp('2021-02-28 00:00:00'),
 Timestamp('2021-03-07 00:00:00'),
 Timestamp('2021-03-14 00:00:00'),
 Timestamp('2021-03-21 00:00:00'),
 Timestamp('2021-03-28 00:00:00'),
 Timestamp('2021-04-04 00:00:00'),
 Timestamp('2021-04-11 00:00:00'),
 Timestamp('2021-04-18 00:00:00'),
 Timestamp('2021-04-25 00:00:00'),
 Timestamp('2021-05-02 00:00:00'),
 Timestamp('2021-05-09 00:00:00'),
 Timestamp('2021-05-16 00:00:00'),
 Timestamp('2021-05-23 00:00:00'),
 Timestamp('2021-05-30 00:00:00'),
 Timestamp('2021-06-06 00:00:00'),
 Timestamp('2021-06-13 00:00:00'),
 Timestamp('2021-06-20 00:00:00'),
 Timestamp('2021-06-27 00:00:00'),
 Timestamp('2021-07-04 00:00:00'),
 Timestamp('2021-07-11 00:00:00'),
 Timestamp('2021-07-18 00:00:00'),
 Timestamp('2021-07-25 00:00:00'),
 Timestamp('2021-08-01 00:00:00'),
 Timestamp('2021-08-08 00:00:00'),
 Timestamp('2021-08-15 00:00:00'),
 Timestamp('2021-08-22 00:00:00'),
 Timestamp('2021-08-29 00:00:00'),
 Timestamp('2021-09-05 00:00:00'),
 Timestamp('2021-09-12 00:00:00'),
 Timestamp('2021-09-19 00:00:00'),
 Timestamp('2021-09-26 00:00:00'),
 Timestamp('2021-10-03 00:00:00'),
 Timestamp('2021-10-10 00:00:00'),
 Timestamp('2021-10-17 00:00:00'),
 Timestamp('2021-10-24 00:00:00'),
 Timestamp('2021-10-31 00:00:00'),
 Timestamp('2021-11-07 00:00:00'),
 Timestamp('2021-11-14 00:00:00'),
 Timestamp('2021-11-21 00:00:00'),
 Timestamp('2021-11-28 00:00:00'),
 Timestamp('2021-12-05 00:00:00'),
 Timestamp('2021-12-12 00:00:00'),
 Timestamp('2021-12-19 00:00:00'),
 Timestamp('2021-12-26 00:00:00'),
 Timestamp('2022-01-02 00:00:00'),
 Timestamp('2022-01-09 00:00:00'),
 Timestamp('2022-01-16 00:00:00'),
 Timestamp('2022-01-23 00:00:00'),
 '_uid_',
 'entidad_cvegeo',
 'municipio_cvegeo',
 'municipio_nombre',
 'pob',
 'geometry']

y = confirmados_por_semana_wide.iloc[:,20:107].to_numpy(copy=True) # sacamos el array de los valores: una fila por observación, las variables en columnas
y = (y / confirmados_por_semana_wide['pob'].to_numpy().reshape(2456,1)) * 1000 # tasa por 1000 habitantes
sm = giddy.markov.Spatial_Markov(y, w, fixed=True, k=5, m=5) # instanciamos el modelo
sm.summary()

--------------------------------------------------------------
                     Spatial Markov Test                      
--------------------------------------------------------------
Number of classes: 5
Number of transitions: 211216
Number of regimes: 5
Regime names: LAG0, LAG1, LAG2, LAG3, LAG4
--------------------------------------------------------------
   Test                   LR                Chi-2
  Stat.            16728.631            20038.130
    DOF                   80                   80
p-value                0.000                0.000
--------------------------------------------------------------
P(H0)           C0         C1         C2         C3         C4
     C0      0.725      0.048      0.112      0.074      0.041
     C1      0.220      0.391      0.327      0.055      0.007
     C2      0.185      0.112      0.461      0.206      0.036
     C3      0.118      0.019      0.202      0.491      0.170
     C4      0.067      0.002      0.035      0.163      0.733
--------------------------------------------------------------
P(LAG0)         C0         C1         C2         C3         C4
     C0      0.819      0.040      0.074      0.041      0.026
     C1      0.321      0.420      0.231      0.025      0.004
     C2      0.378      0.119      0.387      0.103      0.013
     C3      0.412      0.021      0.177      0.314      0.076
     C4      0.460      0.004      0.045      0.112      0.379
--------------------------------------------------------------
P(LAG1)         C0         C1         C2         C3         C4
     C0      0.700      0.067      0.131      0.072      0.031
     C1      0.234      0.427      0.299      0.036      0.004
     C2      0.224      0.144      0.463      0.152      0.018
     C3      0.203      0.029      0.246      0.425      0.098
     C4      0.236      0.004      0.068      0.194      0.497
--------------------------------------------------------------
P(LAG2)         C0         C1         C2         C3         C4
     C0      0.634      0.057      0.159      0.102      0.048
     C1      0.175      0.363      0.383      0.072      0.007
     C2      0.147      0.114      0.494      0.216      0.028
     C3      0.116      0.021      0.233      0.498      0.132
     C4      0.114      0.002      0.064      0.248      0.572
--------------------------------------------------------------
P(LAG3)         C0         C1         C2         C3         C4
     C0      0.617      0.036      0.142      0.129      0.076
     C1      0.137      0.325      0.430      0.097      0.012
     C2      0.115      0.091      0.473      0.271      0.050
     C3      0.064      0.018      0.199      0.535      0.185
     C4      0.053      0.003      0.044      0.216      0.685
--------------------------------------------------------------
P(LAG4)         C0         C1         C2         C3         C4
     C0      0.581      0.020      0.117      0.149      0.133
     C1      0.103      0.308      0.402      0.147      0.039
     C2      0.102      0.068      0.432      0.303      0.095
     C3      0.058      0.011      0.153      0.515      0.263
     C4      0.024      0.001      0.021      0.122      0.831
--------------------------------------------------------------

Lo que vemos en el resultado son las probabiliades de transición condicionadas al contexto local de los municipios. Cada LAG corresponde a una clasificación de los vecinos de los municipios en uno de los cuartiles: LAG0 el más bajo y LAG3 el más alto, más adelante vamos a visualizar estas probabilidades de una mejor forma. Por lo pronto lo importante es que la hipótesis nula (que el contexto de los municipios no es importante) se rechaza debido a los valores de $LR$, $CHI-2$ y sus respectivos $p-values$.

Para ver mejor el resultado del modelo, hagamos heatmaps con las matrices de transición

fig, axes = plt.subplots(2,3,figsize = (15,8))
for i in range(2):
    for j in range(3):
        ax = axes[i,j]
        if i==0 and j==0:
            im = sns.heatmap(sm.p, annot=True, linewidths=.5, ax=ax, cbar=True, vmin=0, vmax=1,
                          square=True, cmap="YlGn",fmt='.3f')
            ax.set_title("Global",fontsize=13) 
        else:

            p_temp = sm.P[i*3+j-1]
            im = sns.heatmap(p_temp, annot=True, linewidths=.5, ax=ax, cbar=True, vmin=0, vmax=1,
                          square=True, cmap="YlGn",fmt='.3f')
            ax.set_title("Retraso espacial %d"%(i*3+j),fontsize=13)

El primer heatmap representa la matriz global de transiciones, sin considerar el contexto. Las siguientes matrices son las transiciones condicionadas al estado de los municipios vecinos.

Lo primero que podemos ver es que los municipios en los cuantiles más altos cambián su comportamiento de forma siginificativa en función de sus vecinos, cuando sus vecinos están también en el cuantil más alto (Retraso espacial 5), tienden a permanecer así, pero cunado sus vecinos están en los cuantiles bajos, entonces pueden moverse al cuantil más bajo también.

Otra cosa interesante es que para los municipios en el cuartil más bajo, el estado de los vacinos produce transiciones hacia los cuartiles 2 y 3. No parece haber evidencia de que los municipios en el cuantil 1 vean afectadsas sus transiciones al segundo cuartil, pero sí para los que están el el cuantil 2 que aumentan la probabilidad de subir al 3 conforme sus vecinos tienen mayores tasas.

De alguna forma estos resultados, son un poco contraintuitivos, uno esperaría ver algunas transiciones hacia abajo y una mayor probabilidad de subir para los municipios más bajos. Quizá lo que estamos viendo es en parte conseuencia de la forma en la que cuantizamos los datos. Usar Fixed=True calcula los cuartiles tomando los datos de las $n*t$ observaciones, alternativamente, podemos cuantizar utilizando ínicamente las observaciones en cada tiempo. Noten que para esto vamos a tomar una semana menos

y = confirmados_por_semana_wide.iloc[:,20:107].to_numpy(copy=True) # sacamos el array de los valores: una fila por observación, las variables en columnas
y = (y / confirmados_por_semana_wide['pob'].to_numpy().reshape(2456,1)) * 1000 # tasa por 1000 habitantes
sm_r = giddy.markov.Spatial_Markov(y, w, fixed=False, k=5, m=5) # spatial_markov instance o
sm_r.summary()

--------------------------------------------------------------
                     Spatial Markov Test                      
--------------------------------------------------------------
Number of classes: 5
Number of transitions: 211216
Number of regimes: 5
Regime names: LAG0, LAG1, LAG2, LAG3, LAG4
--------------------------------------------------------------
   Test                   LR                Chi-2
  Stat.            16304.593            19187.827
    DOF                   80                   80
p-value                0.000                0.000
--------------------------------------------------------------
P(H0)           C0         C1         C2         C3         C4
     C0      0.724      0.047      0.112      0.073      0.043
     C1      0.200      0.412      0.319      0.060      0.009
     C2      0.184      0.120      0.449      0.212      0.036
     C3      0.117      0.021      0.208      0.477      0.177
     C4      0.073      0.003      0.037      0.174      0.714
--------------------------------------------------------------
P(LAG0)         C0         C1         C2         C3         C4
     C0      0.817      0.039      0.075      0.042      0.027
     C1      0.304      0.454      0.213      0.026      0.004
     C2      0.372      0.130      0.367      0.113      0.018
     C3      0.408      0.032      0.195      0.280      0.084
     C4      0.447      0.003      0.046      0.124      0.379
--------------------------------------------------------------
P(LAG1)         C0         C1         C2         C3         C4
     C0      0.705      0.061      0.129      0.071      0.035
     C1      0.210      0.447      0.301      0.037      0.006
     C2      0.227      0.153      0.449      0.153      0.019
     C3      0.209      0.030      0.246      0.412      0.103
     C4      0.217      0.009      0.076      0.202      0.497
--------------------------------------------------------------
P(LAG2)         C0         C1         C2         C3         C4
     C0      0.643      0.056      0.157      0.098      0.047
     C1      0.154      0.393      0.369      0.076      0.007
     C2      0.156      0.119      0.482      0.213      0.030
     C3      0.115      0.023      0.237      0.483      0.142
     C4      0.111      0.006      0.061      0.253      0.570
--------------------------------------------------------------
P(LAG3)         C0         C1         C2         C3         C4
     C0      0.607      0.043      0.150      0.126      0.074
     C1      0.138      0.335      0.399      0.109      0.019
     C2      0.115      0.100      0.469      0.273      0.044
     C3      0.069      0.020      0.206      0.519      0.185
     C4      0.055      0.003      0.046      0.215      0.681
--------------------------------------------------------------
P(LAG4)         C0         C1         C2         C3         C4
     C0      0.591      0.026      0.112      0.143      0.127
     C1      0.112      0.329      0.408      0.121      0.030
     C2      0.092      0.082      0.424      0.317      0.085
     C3      0.050      0.013      0.163      0.507      0.267
     C4      0.026      0.001      0.021      0.135      0.817
--------------------------------------------------------------

El modelo sigue siendo estadísticamente significativo

fig, axes = plt.subplots(2,3,figsize = (15,8)) 

for i in range(2):
    for j in range(3):
        ax = axes[i,j]
        if i==0 and j==0:
            im = sns.heatmap(sm_r.p, annot=True, linewidths=.5, ax=ax, cbar=True, vmin=0, vmax=1,
                          square=True, cmap="YlGn",fmt='.3f')
            ax.set_title("Global",fontsize=13) 
        else:
            p_temp = sm_r.P[i*3+j-1]
            im = sns.heatmap(p_temp, annot=True, linewidths=.5, ax=ax, cbar=True, vmin=0, vmax=1,
                          square=True, cmap="YlGn",fmt='.3f')
            ax.set_title("Retraso espacial %d"%(i*3+j),fontsize=13)
                

Estos resultados tienen mucho más sentido. Por ejemplo, globalmente, la probabilidad de pasar del cuartil más bajo al siguiente es 0.047, sin embargo, condicionada a que sus vecinos esten en el tercer cuartil es de 0.056. En el sentido contrario, la probabilidad de permanecer en el cuartil más alto es más baja cuando los vecinos se encuentran en el cuartil más bajo.

LISA Markov

Finalmente, una alternativa es utilizar como discretización los cuadrantes de LISA para las observaciones a nivel municipio. Tomemos, por ejemplo, el diagrama de dispersión de Moran para la última semana

ultima_confirmados = confirmados_por_semana_wide[[pd.to_datetime('2021-01-03'), 'municipio_cvegeo', 'geometry', 'pob']].copy()
ultima_confirmados['tasa'] = (ultima_confirmados[pd.to_datetime('2021-01-03')] / ultima_confirmados['pob'])*1000
local_moran = Moran_Local(ultima_confirmados['tasa'].to_numpy(), w, permutations=1000)
fig, ax = plt.subplots(figsize=(12, 12))
fig, ax = moran_scatterplot(local_moran, p=0.05, ax=ax)
ax.set_xlabel('Tasa de Confirmados')
ax.set_ylabel('Retraso Espacial')
plt.show()

El diagrama de dispersión de Moran nos muestra la dispersión de los valores de la variable de interés para cada observación contra el promedio de sus vecinos. Los cuatro cuadrantes (líneas punteadas) definen el tipo de cluster:

• HH: arriba a la izquierda
• LL: abajo a la izquierda
• LH: abajo a la derecha
• HL: arriba a la izquierda

Entonces, para cada momento en el tiempo, cada observación queda clasificada en alguno de esos cuadrantes, lo que nos da una cuantización alternativa de los datos que incluye ya la información del contexto geográfico.

Primero ajustamos el modelo y vemos si es significativo.

Nota: vamos a usar menos semanas para quedarnos con semanas que tengan cantidades importantes de conglomerados significativos

y = confirmados_por_semana_wide.iloc[:,20:107].to_numpy(copy=True) # sacamos el array de los valores: una fila por observación, las variables en columnas
y = (y / confirmados_por_semana_wide['pob'].to_numpy().reshape(2456,1)) * 1000 # tasa por 1000 habitantes
lm = giddy.markov.LISA_Markov(y, w, permutations=1000)
"Chi2: %8.3f, p: %5.2f, dof: %d" % lm.chi_2

/home/plablo/miniconda3/envs/geoinformatica/lib/python3.7/site-packages/giddy/markov.py:1399: RuntimeWarning: invalid value encountered in true_divide
  r = y / ybar
/home/plablo/miniconda3/envs/geoinformatica/lib/python3.7/site-packages/giddy/markov.py:1401: RuntimeWarning: divide by zero encountered in true_divide
  rlag = ylag / ybar
/home/plablo/miniconda3/envs/geoinformatica/lib/python3.7/site-packages/giddy/markov.py:1401: RuntimeWarning: invalid value encountered in true_divide
  rlag = ylag / ybar

'Chi2: 144786.662, p:  0.00, dof: 9'

Estamos ajustando el modelo a partir de la semana 21 en adelante y de acuerdo a los valores de $Chi^2$ y de $p$ es calro que se descarta la hipótesis de que las trancisiones no dependen del retraso espacial.

Veamos entonces la matriz de transiciones

fig, ax = plt.subplots(figsize = (7,6))
im = sns.heatmap(lm.p, annot=True, linewidths=.5, ax=ax, cbar=True, vmin=0, vmax=1,
                 square=True,  cmap="YlGn",fmt='.3f', xticklabels=['HH', 'LH', 'LL', 'HL'],
                 yticklabels=['HH', 'LH', 'LL', 'HL'])

La diagonal domina las probabilidades, aunque no por mucho, lo cual quiere decir que no es un régimen estable. Sin embargo es interesante ver que tanto los LH como los HL son menos estables. Para los LH la probabilidad de pasar a HH o a LL es casi igual, sin embargo para los HL sí es bastante más probable pasar a LL. Esto nos dice que los municipios que están en mayor riesgo son en definitiva los LH. Para darnos una idea de cuáles son, veamos una vez más el mapa para la última semana

lisa_confirmados = Moran_Local.by_col(ultima_confirmados.drop(pd.to_datetime('2021-01-03'), axis=1).fillna(0),'tasa',  w, permutations=10000, outvals=['q'])
# lisa_ultima = lisa_ultima.loc[:,['municipio_cvegeo', 'geometry', 'tasa,', 'tasa_p_sim', 'tasa_q']]
lisa_confirmados.rename({'tasa_q': 'cluster'}, inplace=True, axis=1)
lisa_confirmados.loc[lisa_confirmados['tasa_p_sim'] >= 0.5, 'cluster'] = -4
lisa_confirmados = gpd.geodataframe.GeoDataFrame(lisa_confirmados)
lisa_confirmados.crs = 'EPSG:4326'
lisa_confirmados['cluster'] = lisa_confirmados['cluster'].replace({1: 'HH', 2: 'LH', 3:'LL', 4:'HL', -4:'NS'})
lisa_confirmados.head()
color_dict = {'HH':'red', 'LL':'blue', 'HL':'orange', 'LH':'cornflowerblue', 'NS':'lightgray'}
fig, ax = plt.subplots(figsize=(14, 16))
for ctype, data in lisa_confirmados.groupby('cluster'):
    color = color_dict[ctype]
    # Plot each group using the color defined above
    data.to_crs(epsg=3857).plot(color=color,
              ax=ax,
              label=ctype,
              alpha=0.8 )

red_patch = mpatches.Patch(color='red', label='HH')
blue_patch = mpatches.Patch(color='blue', label='LL')
orange_patch = mpatches.Patch(color='orange', label='HL')
lblue_patch = mpatches.Patch(color='cornflowerblue', label='LH')
gray_patch = mpatches.Patch(color='lightgray', label='NS')
plt.legend(handles=[red_patch, blue_patch, orange_patch, lblue_patch, gray_patch], loc=(0.75,0.5))
ax.set(title='Conglomerados espaciales para la tasa de confirmados')
ax.set_axis_off()
ctx.add_basemap(ax, source=ctx.providers.OpenStreetMap.Mapnik)
plt.show()