martes, 12 de mayo de 2020

Programando un comando de estimación en Stata: Funciones Mata.

David M. Drukker, Director Ejecutivo de Econometría.

Mostraré cómo escribir una función en Mata, el lenguaje de programación matricial que forma parte de Stata. Esta publicación utiliza conceptos introducidos en Programación de un comando de estimación en Stata: Mata 101.

Esta es la décimo segunda publicación de la serie Programando un comando de estimación en Stata. Te recomiendo que empieces por el principio.

Funciones Mata

Los comandos funcionan en Stata. Las funciones funcionan en Mata. Los comandos operan en objetos Stata, como variables, y los usuarios especifican opciones para alterar el comportamiento. Las funciones de Mata aceptan argumentos, operan en los argumentos y pueden devolver un resultado o alterar el valor de un argumento para contener un resultado.

Considere myadd() definido a continuación.


myadd() acepta dos argumentos, X y Y, coloca la suma de X y Y dentro de A, y devuelve A. Por ejemplo,

Ejemplo 1: Definir y usar una función




Después de definir myadd(), creo las matrices C y D, y paso C y D a myadd(), que devuelve su suma. Mata asigna la suma devuelta a W, la cual muestro. Tenga en cuenta que dentro de la función myadd(), C y D se conocen respectivamente como X e Y.

Solo se puede acceder a la A creada en myadd() dentro de myadd(), y no entra en conflicto con una A definida fuera de myadd(); es decir, A es local a la función myadd(). Ahora ilustraré que A es local para myadd().

Ejemplo 2: A es local a myadd()



Habiendo ilustrado que la A definida dentro de myadd() es local a myadd(), debo señalar que las matrices C y D que definí en el ejemplo 1 están en la memoria global de Mata. Al igual que en los programas ado, no queremos utilizar nombres fijos en la memoria global de Mata en nuestros programas para no destruir los datos de los usuarios. Afortunadamente, este problema se resuelve fácilmente escribiendo funciones Mata que escriben sus resultados en Stata y no devuelven resultados. Proporcionaré discusiones detalladas de esta solución en los comandos que desarrollé en publicaciones posteriores.

Cuando una función Mata cambia el valor de un argumento dentro de la función, eso cambia el valor de ese argumento fuera de la función; en otras palabras, los argumentos se pasan por dirección. Las funciones de Mata pueden calcular más de un resultado almacenando estos resultados en argumentos. Por ejemplo, sumproduct() devuelve tanto la suma como el producto basado en elementos de dos matrices.



sumproduct() retorna la suma de los argumento X y Y en el argumento S y el producto por elementos en P.

Ejemplo 3: Devolviendo resultados en argumentos






Después de definir sumproduct(), uso I() para crear A y uso rowshape() para crear B. Luego creo W y Z; cada uno es un valor escalar perdido. Debo crear W y Z antes de pasarlos como argumentos; de lo contrario, estaría haciendo referencia a argumentos que no existen. Después de llamar a sumproduct(), muestro W y Z para ilustrar que ahora contienen la suma y el producto de X y Y por elementos.

En myadd() y sumproduct(), no especifiqué qué tipo de cosa debe ser cada argumento, ni especifiqué qué tipo de cosa devolvería cada función. En otras palabras, utilicé declaraciones implícitas. Las declaraciones implícitas son más fáciles de escribir que las declaraciones explícitas, pero hacen que los mensajes de error y el código sean menos informativos. Recomiendo declarar explícitamente devoluciones, argumentos y variables locales para que su código y sus mensajes de error sean más legibles.

myadd2() es una versión de myadd() que declara explícitamente el tipo de cosa devuelta, el tipo de cosa que debe ser cada argumento y el tipo que debe ser cada cosa local a la función.



myadd2() devuelve una matriz numérica que se construye agregando la matriz numérica X a la matriz numérica Y. El objeto local A de la función es también una matriz numérica. Una matriz numérica es una matriz real o una matriz compleja; no puede ser una matriz de cadena.

La declaración explícita de devoluciones, argumentos y variables locales hace que el código sea más informativo. Inmediatamente veo que myadd2() no funciona con matrices de cadenas, pero esta propiedad está oculta en el código de myadd().

No puedo enfatizar lo suficiente la importancia de escribir código fácil de leer. Leer el código de otras personas es una parte esencial de la programación. Es instructivo y le permite adoptar las soluciones que otros han encontrado o implementado. Si eres nuevo en la programación, es posible que aún no te des cuenta de que después de unos meses, leer tu propio código es como leer el código de otra persona. Incluso si nunca le da su código a nadie más, es esencial que escriba un código fácil de leer para poder leerlo en una fecha posterior.

Las declaraciones explícitas también hacen que algunos mensajes de error sean más fáciles de rastrear. En los ejemplos 4 y 5, paso una matriz de cadena a myadd() y myadd2(), respectivamente.

Ejemplo 4: Pasando una matriz de cadena a myadd()



Ejemplo 5: Pasando una matriz de cadena a myadd2()



El mensaje de error en el ejemplo 4 indica que en algún lugar de myadd(), un operador o una función no pudieron realizar algo en dos objetos porque sus tipos no eran compatibles. No se deje engañar por la simplicidad de myadd(). Rastrear una falta de coincidencia de tipos en el código real puede ser difícil.

Por el contrario, el mensaje de error en el ejemplo 5 dice que la matriz C que pasamos a myadd2() no es una matriz real ni compleja como requiere el argumento de myadd2(). Al mirar el código y el mensaje de error, inmediatamente me informa que el problema es que pasé una matriz de cadena a una función que requiere una matriz numérica.

Las declaraciones explícitas son tan recomendables que Mata tiene una configuración que lo requiere, como se ilustra a continuación.

Ejemplo 6: Encendiendo matastrict



Establecer matastrict en on hace que el compilador Mata requiera que las variables de retorno y locales se declaren explícitamente. Por defecto, matastrict está desactivado, en cuyo caso las variables locales y de retorno pueden declararse implícitamente.

Cuando matastrict está activado, no es necesario que los argumentos se declaren explícitamente porque algunos argumentos contienen resultados cuyos tipos de entrada y salida pueden diferir. Considere makereal() definido y usado en el ejemplo 7.

Ejemplo 7: Cambiar un tipo de argumentos




La declaración de makereal() especifica que makereal() no devuelve nada porque void viene antes del nombre de la función. A pesar de que matastrict está activado, no declaré qué tipo de cosa debe ser A. Las dos líneas ejecutables de makereal() aclaran que A debe ser una cadena como elemento de entrada y que A será real como elemento de salida, lo que luego ilustraré.

Utilizo la función de que los argumentos pueden declararse implícitamente para que mi código sea más fácil de leer. Muchas de las funciones de Mata que escribo reemplazan los argumentos con resultados. Declaro explícitamente argumentos que son entradas, y declaro implícitamente argumentos que contienen salidas. Considere sumproduct2().


 


sumproduct2() no devuelve nada porque void viene antes del nombre de la función. El primer argumento X es una matriz real, el segundo argumento Y es una matriz real, el tercer argumento S se declara implícitamente y el cuarto argumento P se declara implícitamente. Mi convención de codificación de que las entradas se declaran explícitamente y que las salidas se declaran implícitamente de inmediato me informa que X y Y son entradas pero que S y P son salidas. Que X y Y son entradas y que S y P son salidas se ilustra en el ejemplo 8.

Ejemplo 8: Declarando explícitamente entradas, pero no salidas



Hecho y sin hacer

Mostré cómo escribir una función en Mata y discutí las declaraciones con cierto detalle. Escriba help m2_declarations para muchos más detalles.

En mi próxima publicación usaré las funciones de Mata para realizar los cálculos para un comando de estimación simple.

¡Gracias por leernos! 

Drukker, David. (22 de diciembre de 2015). Programming an estimation command in Stata: Mata functions. (Trad. Ángel Cruz). The Stata Blog. Not Elsewhere Classified. Recuperado de https://blog.stata.com/2015/12/22/programming-an-estimation-command-in-stata-mata-functions/

lunes, 13 de abril de 2020

Cómo crear mapas animados de coropletas usando los datos COVID-19 de la Universidad Johns Hopkins.


Chuck Huber, Director Asociado de Alcance Estadístico.

En mis publicaciones anteriores, mostré cómo descargar los datos COVID-19 del repositorio GitHub Johns Hopkins, graficar los datos a lo largo del tiempo y crear mapas coropléticos. Ahora, le mostraré cómo crear mapas animados de coropletas para explorar la distribución de COVID-19 a lo largo del tiempo y el lugar.

El siguiente video muestra la cantidad acumulada de casos de COVID-19 por cada 100,000 habitantes para cada condado en los Estados Unidos desde el 22 de enero de 2020 hasta el 5 de abril de 2020. El mapa no cambia mucho hasta mediados de marzo, cuando el virus comienza a extenderse más rápido. Entonces, podemos ver cuándo y dónde las personas están siendo infectadas. Puede hacer clic en el icono "Reproducir" en el video para reproducirlo y hacer clic en el icono en la parte inferior derecha para ver el video en modo de pantalla completa.



En mi último post, aprendimos cómo crear un mapa coroplético. Tendremos que aprender dos habilidades adicionales para crear un mapa animado: cómo crear un mapa para cada fecha en el conjunto de datos y cómo combinar la colección de mapas en un archivo de video.

Cómo crear un mapa para cada fecha
Comencemos importando y describiendo los datos sin procesar del repositorio GitHub de Johns Hopkins. Tenga en cuenta que importé estos datos el 5 de abril.



Las variables de v12 a v86 contienen el número acumulado de casos confirmados de COVID-19 en cada condado para cada día a partir del 22 de enero de 2020 y hasta el 5 de abril de 2020. Enlistemos algunas observaciones para ver los datos sin procesar.




Tendremos que incluir de v12 a v86 en nuestro conjunto de datos final. Copié el bloque de código desde el final de mi última publicación y lo pegué a continuación con dos pequeñas modificaciones. La primera modificación se muestra en rojo. El código retiene y formatea los datos del caso para cada fecha guardada en v12 a v86. Tenga en cuenta que me refiero a "v12 a v86" en el código escribiendo v*. El asterisco sirve como comodín, por lo que v* se refiere a cualquier variable que comienza con la letra "v". La segunda modificación es que no creamos una variable para nuestro recuento ajustado por la población en este momento.



// Create the geographic dataset
clear
copy https://www2.census.gov/geo/tiger/GENZ2018/shp/cb_2018_us_county_500k.zip ///
    cb_2018_us_county_500k.zip
unzipfile cb_2018_us_county_500k.zip
spshape2dta cb_2018_us_county_500k.shp, saving(usacounties) replace
use usacounties.dta, clear
generate fips = real(GEOID)
save usacounties.dta, replace
// Create the COVID-19 case dataset
clear
import delimited https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_US.csv
drop if missing(fips)
save covid19_county, replace
// Create the population dataset
clear
import delimited https://www2.census.gov/programs-surveys/popest/datasets/2010-2019/counties/totals/co-est2019-alldata.csv
generate fips = state*1000 + county
save census_popn, replace
// Merge the datasets
clear
use _ID _CX _CY GEOID fips using usacounties.dta
merge 1:1 fips using covid19_county ///
    
, keepusing(province_state combined_key v*)
keep if _merge == 3
drop _merge
merge 1:1 fips using census_popn ///
    , keepusing(census2010pop popestimate2019)
keep if _merge==3
drop _merge
drop if inlist(province_state, "Alaska", "Hawaii")
format %16.0fc popestimate2019 v*
save covid19_adj, replace

Describamos algunos de nuestros conjuntos de datos finales.


El conjunto de datos contiene todas las variables que necesitaremos para crear nuestro mapa animado: la información geográfica, la población de cada condado y el número de casos para cada día en v12 a v86. Voy a explicar todos los pasos que necesitaremos seguir usando solo tres de estas variables. Luego, los juntaremos usando las 75 variables.

Tendremos que crear un mapa separado para v12, v13, v14, y así sucesivamente hasta v86. Recordemos que aprendimos acerca de los bucles en una de mis publicaciones anteriores. Podemos usar forvalues para realizar ciclos sobre estas variables.

El siguiente ciclo forvalues almacenará los números 12, 13 y 14 en la macro local time. Puede referirse a time dentro del ciclo usando comillas simples izquierda y derecha. El siguiente ejemplo describe cada variable.




También necesitaremos calcular el número de casos ajustados por la población para cada condado para cada día. En el siguiente ejemplo, comenzamos generando una variable llamada confirmed_adj que contiene valores faltantes. Dentro del bucle, reemplazamos confirmed_adj con el recuento ajustado para el día actual y resumimos confirmed_adj. Luego, borramos confirmed_adj cuando finalice el ciclo.




También podríamos querer incluir la fecha en el título de cada mapa. Recuerde que las fechas se almacenan como el número de días desde el 1 de enero de 1960. La variable v12 contiene datos para el 22 de enero de 2020. Podemos calcular el número de días desde el 1 de enero de 1960 utilizando la función date().





Nuestro bucle forvalues comienza en time = 12, y nos gustaría que la fecha para v12 sea 21936. Por lo tanto, debemos restar 12 de la fecha antes de sumar el valor de la macro local time. El siguiente ejemplo recorre los días del 22 de enero de 2020 al 24 de enero de 2020.





Entonces podemos usar la función string() para cambiar el formato de visualización de cada fecha.




También usaremos la exportación de gráficos para exportar nuestros gráficos a archivos Portable Network Graphics (.png). Estos archivos se combinarán más adelante para crear nuestro video, y los nombres de los archivos deben numerarse secuencialmente desde "1" con ceros a la izquierda. El siguiente ejemplo muestra cómo usar la función string() para crear estos nombres de archivo.



Vamos a juntar todas las piezas y mostrar los comandos básicos que usaremos para crear cada mapa.



Cada iteración del bucle hace tres cosas. La primera línea reemplaza la variable confirmed_adj con el recuento ajustado por la población para ese día en particular. La segunda línea usa grmap para crear el mapa con la fecha en el título. La tercera línea exporta el gráfico a un archivo .png con nombres de archivo secuenciales que comienzan en 001.


El siguiente bloque de código hace lo mismo, pero con algunas opciones agregadas a grmap y graph export. La mayoría de las opciones de grmap se discutieron en mi última publicación. He agregado la opción ocolor(gs8) para representar el color del contorno del mapa con una escala de grises claros y la opción osize(vthin) para hacer que las líneas del mapa sean más delgadas. Las opciones width(3840) y height(2160) después de graph export especifican que cada mapa se guarde con una resolución de 4K. Esto creará imágenes claras y detalladas para nuestro video.

generate confirmed_adj = .
forvalues time = 12/86 {
    local date = 21936 - 12 + `time'
    local date = string(`date', "%tdMonth_dd,_CCYY")
    replace confirmed_adj = 100000*(v`time'/popestimate2019)
    grmap confirmed_adj,                                                    ///
        clnumber(8)                                               ///
        clmethod(custom)                                          ///
        clbreaks(0 5 10 15 20 25 50 100 5000)                     ///
        ocolor(gs8)        osize(vthin)                                            ///
        title("Confirmed Cases of COVID-19 in the United States on `date'") ///
        subtitle("cumulative cases per 100,000 population")
    local filenum = string(`time'-11,"%03.0f")
    graph export "map_`filenum'.png", as(png) width(3840) height(2160)
}
drop confirmed_adj

El bloque de código anterior creará 75 mapas y los guardará en 75 archivos.





Veamos el archivo map_001.png, que es el mapa del 22 de enero de 2020.


Figura 1: map_001.png



A continuación, veamos el archivo map_075.png, que es el mapa del 5 de abril de 2020.




Figura 2: map_075.png




Note que la leyenda en la parte inferior izquierda de ambos mapas es la misma. Utilicé la opción clmethod(custom) con grmap para especificar mis propios puntos de corte para los valores de confirmed_adj. La leyenda cambiará de un día a otro si utiliza el clmethod() predeterminado, que selecciona los puntos de corte en función de los cuantiles de confirmed_adj. Los cuantiles cambiarán de un día a otro, y la leyenda no será coherente en todos los mapas. Seleccioné los puntos de corte en clbreaks (0 5 10 15 20 25 50 100 5000) según el mapa del día final.

Puede actualizar su video con datos futuros aumentando el límite superior del ciclo forvalues. Por ejemplo, el 12 de abril de 2020, escribiría forvalues time = 12/91.

Creamos con éxito un mapa para cada día y guardamos cada mapa en un archivo separado. Ahora, necesitamos combinar estos archivos para crear nuestro video.

Cómo crea un vídeo a partir de la colección de mapas

Escribí una publicación de blog en 2014 que describe cómo usar FFmpeg para crear un video a partir de una colección de imágenes. FFmpeg es un paquete de software gratuito con versiones disponibles para Linux, Mac y Windows. Puede ejecutar comandos FFmpeg desde Stata usando shell.

El siguiente ejemplo usa FFmpeg para combinar nuestros archivos de mapas en un video llamado covid19.mp4.

shell "ffmpeg.exe" -framerate 1/.5 -i map_%03d.png -c:v libx264 -r 30 -pix_fmt yuv420p covid19.mp4

Es posible que deba especificar la ubicación de ffmpeg.exe en su computadora como en el ejemplo a continuación.

shell "C:\Program Files\ffmpeg\bin\ffmpeg.exe" -framerate 1/.5 -i map_%03d.png -c:v libx264 -r 30 -pix_fmt yuv420p covid19.mp4

Los nombres de los archivos de imagen se especifican con la opción -i map_% 03d.png, y el nombre del video de salida, covid19.mp4, se especifica al final del comando. La velocidad de fotogramas se especifica con la opción -framerate 1/.5. La relación "1 /.5" especifica una velocidad de cuadro de 2 imágenes por segundo.

FFmpeg puede tardar varios minutos en procesar las imágenes y convertirlas en un video. Su paciencia será recompensada con el siguiente video.






Conclusiones y código recopilado

¡Lo hicimos! Creamos con éxito un mapa animado de coropletas que muestra los casos confirmados ajustados por la población de COVID-19 para cada condado en los Estados Unidos desde el 22 de enero de 2020 hasta el 5 de abril de 2020. Nuestro video es una herramienta poderosa que podemos usar para estudiar La distribución de COVID-19 a lo largo del tiempo y la ubicación.
Me gustaría recordarle que he tomado decisiones algo arbitrarias al limpiar los datos que se utilizan para crear este video. Mi objetivo era mostrarte cómo crear tus propios mapas y videos. Mis resultados deben usarse solo con fines educativos, y deberá verificar los datos cuidadosamente si planea usar los resultados para tomar decisiones.

Puede reproducir el video con el siguiente código.


// Create the geographic dataset
clear
copy https://www2.census.gov/geo/tiger/GENZ2018/shp/cb_2018_us_county_500k.zip ///
    cb_2018_us_county_500k.zip
unzipfile cb_2018_us_county_500k.zip
spshape2dta cb_2018_us_county_500k.shp, saving(usacounties) replace
use usacounties.dta, clear
generate fips = real(GEOID)
save usacounties.dta, replace
// Create the COVID-19 case dataset
clear
import delimited https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_US.csv
drop if missing(fips)
save covid19_county, replace
// Create the population dataset
clear
import delimited https://www2.census.gov/programs-surveys/popest/datasets/2010-2019/counties/totals/co-est2019-alldata.csv
generate fips = state*1000 + county
save census_popn, replace
// Merge the datasets
clear
use _ID _CX _CY GEOID fips using usacounties.dta
merge 1:1 fips using covid19_county ///
    , keepusing(province_state combined_key v*)
keep if _merge == 3
drop _merge
merge 1:1 fips using census_popn ///
    , keepusing(census2010pop popestimate2019)
keep if _merge==3
drop _merge
drop if inlist(province_state, "Alaska", "Hawaii")
format %16.0fc popestimate2019 v*
save covid19_adj, replace
// Create the maps
spset, modify shpfile(usacounties_shp)
generate confirmed_adj = .
forvalues time = 12/86 {
    local date = 21936 - 12 + `time'
    local date = string(`date', "%tdMonth_dd,_CCYY")
    replace confirmed_adj = 100000*(v`time'/popestimate2019)
    grmap confirmed_adj,                                                    ///
        clnumber(8)                                                         ///
        clmethod(custom)                                                    ///
        clbreaks(0 5 10 15 20 25 50 100 5000)                               ///
        ocolor(gs8) osize(vthin)                                            ///
        title("Confirmed Cases of COVID-19 in the United States on `date'") ///
        subtitle("cumulative cases per 100,000 population")
    local filenum = string(`time'-11,"%03.0f")
    graph export "map_`filenum'.png", as(png) width(3840) height(2160)
}
drop confirmed_adj
// Create the video with FFmpeg
shell "C:\Program Files\ffmpeg\bin\ffmpeg.exe" -framerate 1/.5 -i map_%03d.png -c:v libx264 -r 30 -pix_fmt yuv420p covid19.mp4



Huber, Chuck.
(10 de abril del 2020).
How to create animated choropleth maps using the COVID-19 data from Johns Hopkins University. The Stata Blog. Not Elsewhere Classified.



Eso es todo por hoy. ¡Gracias por leernos!

miércoles, 8 de abril de 2020

Cómo crear mapas coropléticos utilizando los datos COVID-19 de la Universidad Johns Hopkins.


Chuck Huber, Director Asociado de Alcance Estadístico.

En mi última publicación, aprendimos cómo importar los datos sin procesar de COVID-19 del repositorio GitHub de Johns Hopkins y convertir los datos sin procesar en datos de series temporales. Esta publicación demostrará cómo descargar datos en bruto y crear mapas coropléticos como la figura 1.

Figura 1: Casos confirmados de COVID-19 en los Estados Unidos, ajustados por tamaño de población.




Un mapa coroplético muestra información estadística sobre áreas geográficas utilizando diferentes colores o intensidad de sombreado. La Figura 1 muestra el número ajustado por la población de casos confirmados de COVID-19 para cada condado en los Estados Unidos a partir del 2 de abril de 2020. Cada tono de azul en el mapa representa el rango del número de casos que se muestran en la leyenda inferior izquierda del mapa. Utilicé el comando grmap contribuido por la comunidad para crear la figura 1, y necesitamos tres datos sobre cada condado para crear este mapa. Necesitamos la información geográfica de cada condado, el número de casos confirmados en cada condado y la población de cada condado.

Comencemos por el final y avancemos hacia atrás para aprender cómo construir un conjunto de datos que contenga esta información. Los datos enlistados a continuación se utilizaron para crear el mapa en la figura 1. Cada observación contiene información sobre un condado individual en los Estados Unidos.





Las primeras cuatro variables contienen información geográfica sobre cada condado. La variable _ID contiene un número de identificación único para cada condado que se utiliza para vincular con un archivo especial llamado "shapefile". Los shapefiles contienen la información que se utiliza para representar el mapa y se explicarán a continuación. La variable GEOID es un código de condado del Estándar Federal de Procesamiento de Información (FIPS). Podemos usar el código FIPS como una variable clave para fusionar datos de otros archivos. Las variables _CX y _CY contienen coordenadas geográficas. Podemos descargar estos datos de la Oficina del Censo de los Estados Unidos.


La variable confirmed contiene el número de casos confirmados de COVID-19 en cada condado. Estos datos se descargaron del repositorio GitHub Johns Hopkins. Este no es el mismo archivo que utilizamos en mis publicaciones anteriores.

La variable popestimate2019 contiene la población de cada condado. Estos datos fueron descargados de la Oficina del Censo de los Estados Unidos.
La variable confirmed_adj contiene el número de casos confirmados por 100,000 habitantes. Esta variable se calcula dividiendo el número de casos para cada condado en confirmed entre la población total para cada condado en popestimate2019. El resultado se multiplica por 100,000 para convertirlo en "casos por 100,000 habitantes".

Tendremos que descargar datos de tres fuentes diferentes y fusionar estos archivos en un solo conjunto de datos para construir el conjunto de datos para nuestro mapa. Comencemos descargando y procesando cada uno de estos conjuntos de datos.

Lista de temas

Descargue y prepare los datos geográficos
Descargue y prepare los datos del caso
Descargar y preparar los datos de la población
Cómo fusionar los archivos y calcular recuentos ajustados
Cómo crear el mapa coroplético con grmap

Descargue y prepare los datos geográficos

Comencemos con los datos geográficos. Los shapefiles contienen la información geográfica que grmap usa para crear mapas. Muchos shapefiles están disponibles gratuitamente en internet, y puede localizarlos utilizando un motor de búsqueda. Por ejemplo, busqué los términos "united states shapefile", y el primer resultado me llevó a la Oficina del Censo de los Estados Unidos. Este sitio web contiene shapefiles para los Estados Unidos que muestran los límites de los estados, distritos del Congreso, áreas estadísticas metropolitanas, micropolitanas y muchos otros. Me gustaría subdividir mi mapa de los Estados Unidos por condado, así que me desplacé hacia abajo hasta que encontré el título "County". Me gustaria usar el archivo cb_2018_us_county_500k.zip.



Podemos copiar el archivo desde el sitio web a nuestra unidad local y usar unzipfile para extraer el contenido del archivo. Puede hacer clic con el botón derecho en el archivo en la página web, seleccionar "Copiar dirección de enlace" y pegar la ruta y el nombre de archivo en su comando copy.




Los archivos cb_2018_us_county_500k.shp y cb_2018_us_county_500k.dbf contienen la información geográfica que necesitamos. Podemos usar spshape2dta para procesar la información en estos archivos y crear dos conjuntos de datos Stata llamados usacounties_shp.dta y usacounties.dta.




El archivo usacounties_shp.dta es el shapefile que contiene la información que grmap usará para representar el mapa. No necesitamos hacer nada con este archivo, pero describamos y enlistemos su contenido para ver qué contiene.



El shapefile, usacounties_shp.dta, contiene 1,047,409 coordenadas que definen los límites de cada condado en nuestro mapa. Este archivo también incluye la variable _ID, que se utiliza para vincular estos datos con los datos en usacounties.dta. Necesitaremos este archivo más tarde.


A continuación, usemos y describamos el contenido de usacounties.dta.




Las primeras tres variables contienen información geográfica sobre cada estado. La variable _ID es el identificador de unidad espacial para cada condado que se utiliza para vincular este archivo con el shapefile, usacounties_shp.dta. Las variables _CX y _CY son las coordenadas x e y del centroide del área para cada condado. La variable NAME contiene el nombre del condado para cada observación. La variable GEOID es el código FIPS almacenado como una cadena. Necesitaremos un código FIPS numérico para fusionar este conjunto de datos con otros conjuntos de datos a nivel de condado. Entonces, generemos una variable llamada fips que sea igual al valor numérico de GEOID.



Guardemos nuestros datos geográficos y pasemos a los datos COVID-19.



Descargue y prepare los datos del caso


En mis publicaciones anteriores, aprendimos cómo descargar los datos en bruto de COVID-19 del repositorio GitHub Johns Hopkins. Utilizaremos un conjunto de datos diferente de otra rama del repositorio de GitHub ubicado aquí.



El archivo time_series_covid19_confirmed_US.csv contiene datos de series de tiempo para el número de casos confirmados de COVID-19 para cada condado en los Estados Unidos. Hagamos clic en el nombre del archivo para ver su contenido.

 Cada observación en este archivo contiene los datos de un condado o territorio en los Estados Unidos. Los recuentos confirmados para cada fecha se almacenan en variables separadas. Veamos los datos delimitados por comas haciendo clic en el botón "Raw" junto a la flecha roja. 




Para importar los datos de caso sin procesar, copie la URL y el nombre de archivo de la barra de direcciones en su navegador web, y pegue la dirección web en  import delimited.


Tenga en cuenta que la URL del archivo de datos es larga y se ajusta a una segunda y tercera línea en nuestro comando import delimited. Esta URL debe ser una línea en su comando import delimited. Describamos este conjunto de datos.



Este conjunto de datos contiene 3,253 observaciones que contienen información sobre condados en los Estados Unidos. Hagamos una lista de las primeras 10 observaciones para fips, combine_key, v80, v81, v82 y v83.


La variable fips contiene el código de condado FIPS que usaremos para fusionar este conjunto de datos con la información geográfica en usacounties_shp.dta. La variable combine_key contiene el nombre de cada condado y estado en los Estados Unidos. Y las variables v80, v81, v82 y v83 contienen el número de casos confirmados de COVID-19 desde el 30 de marzo de 2020 hasta el 2 de abril de 2020. Los datos de casos más recientes se almacenan en v83, así que cambiemos su nombre a confirmed.


Encontraremos problemas más adelante cuando fusionamos conjuntos de datos si fips contiene valores faltantes. Así que borremos cualquier observación que falte para fips.


Nuestra base de datos COVID-19 está completa. Guardemos los datos y movámonos a los datos de población.


Descargar y preparar los datos de la población


Podríamos detenernos aquí, fusionar los datos geográficos con la cantidad de casos confirmados y crear un mapa coroplético de la cantidad de casos para cada condado. Pero esto sería engañoso porque las poblaciones de los condados son diferentes. Podríamos preferir informar la cantidad de casos por cada 100,000 personas, y eso requeriría conocer la cantidad de personas en cada condado. Afortunadamente, esos datos están disponibles en el sitio web de la Oficina del Censo de los Estados Unidos.


Podemos seguir los mismos pasos que utilizamos para descargar e importar los datos del caso. Primero, haga clic con el botón derecho en el nombre del archivo en el sitio web, luego seleccione “Copiar dirección de enlace” y use  import delimited para importar los datos.



Normalmente describiría el conjunto de datos en este punto, pero hay 164 variables en este conjunto de datos. Entonces, describiré solo las variables relevantes a continuación.


La variable census2010pop contiene la población de cada condado según el censo de 2010. Pero esa información tiene 10 años. La variable popestimate2019 es una estimación de la población de cada condado en 2019. Usemos esos datos pues son más recientes.


A continuación, enlistamos los datos.


Este conjunto de datos no incluye una variable con un código de condado FIPS. Pero podemos crear una variable que contenga el código FIPS usando las variables state y county. La inspección visual de los datos geográficos en usacounties.dta indica que el código FIPS de condado es el código de estado seguido del código de estado de tres dígitos. Así que creemos nuestra variable de código fips multiplicando el código state por 1000 y luego agregando el código county.


Hagamos una lista de los datos del condado para verificar nuestro trabajo.


Este conjunto de datos contiene la población estimada de cada estado junto con la variable fips que usaremos como una variable clave para fusionar este archivo de datos con los otros archivos de datos. Guardemos nuestros datos.


Cómo fusionar los archivos y calcular recuentos ajustados


Hemos creado tres archivos de datos que contienen la información que necesitamos para crear nuestro mapa coroplético. El archivo de datos usacounties.dta contiene la información geográfica que necesitamos en las variables _ID, _CX y _CY. Recuerde que estos datos están vinculados al archivo shape, usacounties_shp.dta, utilizando la variable _ID. El archivo de datos covid19_county.dta contiene la información sobre el número de casos confirmados de COVID-19 en la variable confirmed. Y el archivo de datos census_popn.dta contiene la información sobre la población de cada condado en la variable popestimate2019.

Necesitamos todas estas variables en el mismo conjunto de datos para crear nuestro mapa. Podemos fusionar estos archivos usando la variable clave fips.


Comencemos usando solo las variables que necesitamos de usacounties.dta.


A continuación, combinemos el número de casos confirmados de covid19_county.dta. La opción keepusing (province_state combined_key_confirmed) especifica que fusionaremos solo las variables province_state, combined_key y confirmed a partir del archivo de datos covid19_state.dta.


La salida nos dice que 3,142 observaciones tenían valores coincidentes de fips en los dos conjuntos de datos. merge también creó una nueva variable en nuestro conjunto de datos llamada _merge, que equivale a 3 para observaciones con valores coincidentes de fips.


La salida también nos dice que 91 observaciones tenían un código fips en los datos geográficos, pero no en los datos del caso. _merge es igual a 1 para estas observaciones. Hagamos una lista de algunas de estas observaciones.


Las primeras siete observaciones tienen información geográfica pero no hay datos para confirmed. Vamos a contar el número de observaciones donde _merge es igual a 1 y done falten los datos de confirmed.



Borremos estas observaciones de nuestro conjunto de datos pues no contienen datos para confirmed.


La salida de merge también nos dice que 109 observaciones tenían un código fips en los datos del caso, pero no en los datos geográficos. _merge es igual a 2 para estas observaciones. Hagamos una lista de algunas de estas observaciones.


Las primeras siete observaciones tienen información de casos pero no datos para _ID o GEOID.


Algunas de las observaciones son de territorios estadounidenses que no son estados. Otras observaciones tienen un valor de combined_key que sugiere que la información del condado no está clara. La inspección visual de estas observaciones sugiere que la mayoría de los casos confirmados para estas observaciones son cero. Podemos verificar esto contando el número de observaciones donde _merge es igual a 2 y confirmed es igual a cero.


El resultado indica que 78 de las 109 observaciones no contienen casos confirmados. Podríamos investigar más estas observaciones si estuviéramos utilizando nuestros resultados para tomar decisiones políticas. Pero estas observaciones son una pequeña proporción de nuestro conjunto de datos, y nuestro objetivo actual es solo aprender a hacer mapas coropléticos. Así que eliminemos estas observaciones y la variable _merge y sigamos adelante.


Describamos el conjunto de datos para verificar que la fusión se realizó correctamente.


A continuación, combinemos la variable popestimate2019 del archivo de datos census_popn.dta.



La salida nos dice que 3,142 observaciones tenían valores coincidentes de fips en los dos conjuntos de datos. merge nuevamente creó una nueva variable en nuestro conjunto de datos llamada _merge, que es igual a 3 para observaciones con un valor coincidente de fips.

El resultado también nos dice que 51 observaciones tenían un código fips en los datos geográficos, pero no en los datos de la población. _merge es igual a 2 para estas observaciones. Hagamos una lista de algunas de estas observaciones.


No tenemos información geográfica o de caso para estas observaciones, por lo tanto, elimínelas de nuestro conjunto de datos y eliminemos la variable _merge.


Ahora, generemos, etiquetemos y formateemos una nueva variable llamada confirmed_adj que contenga el número ajustado por la población de casos confirmados de COVID-19.




Describamos nuestro conjunto de datos para verificar que contiene todas las variables que necesitaremos para crear nuestro mapa.






Nuestro conjunto de datos está completo! Guardemos el conjunto de datos y aprendamos cómo crear un mapa coroplético.


Cómo crear el mapa coroplético con grmap

Utilizaremos el comando grmap, contribuido por la comunidad, para crear nuestro mapa coroplético. Debe activar grmap antes de usarlo por primera vez.


La creación de un mapa que incluya Alaska y Hawái requerirá el uso de opciones que se ajusten a su gran diferencia de tamaño y a que no sean físicamente adyacentes a los otros 48 estados. Me gustaría mantener nuestro mapa lo más simple posible por ahora, así que voy a borrar las observaciones de Alaska y Hawai.



A continuación, debemos decirle a Stata que estos son datos espaciales mediante el uso de spset. La opción modify shpfile(usastates_shp) vinculará nuestros datos con el shapefile, usastates_shp.dta

Recuerde que el shapefile contiene la información que grmap usará para representar el mapa.



Ahora, podemos usar grmap para crear un mapa de coropletas para el número de casos confirmados de COVID-19 ajustado por la población.


Figura 2: Mapa de coropletas usando sectiles.




Por defecto, grmap divide los datos en cuatro grupos basados en cuartiles de confirmed_adj. He usado la opción clnumber(7) para dividir los datos en 7 grupos o sectiles. Puede cambiar el número de grupos utilizando la opción clnumber(#), donde # es el número de categorías de colores.


También puede especificar puntos de corte personalizados utilizando las opciones clmethod(custom) y clbreaks(numlist). El siguiente mapa utiliza puntos de corte personalizados en 0, 5, 10, 15, 20, 25, 50, 100 y 5000. También he agregado un título y un subtítulo.






Figura 2: Mapa coroplético utilizando puntos de corte personalizados.





Conclusiones y código recopilado


¡Lo hicimos! ¡Creamos un mapa coroplético del número ajustado por la población de casos confirmados de COVID-19 en cada condado de los Estados Unidos! Repasemos los pasos básicos. Primero, descargamos los datos geográficos de la Oficina del Censo de los Estados Unidos y los convertimos en archivos de datos Stata usando spshape2dta. En segundo lugar, descargamos, importamos y procesamos los datos COVID-19 del repositorio GitHub Johns Hopkins y los guardamos en un archivo de datos Stata. 

En tercer lugar, descargamos, importamos y procesamos los datos de población de cada condado de la Oficina del Censo de los Estados Unidos y los guardamos en un archivo de datos de Stata. Cuarto, fusionamos los archivos de datos de Stata y calculamos el número de casos COVID-19 ajustados por la población para cada condado. Y quinto, usamos spset para decirle a Stata que estos son datos espaciales, y usamos grmap para crear nuestro mapa coroplético. Puede seguir estos pasos para crear un mapa coroplético para muchos tipos de datos, para otras subdivisiones de los Estados Unidos o para otros países.

He recopilado el siguiente código que reproducirá las figuras 2 y 3.

// Create the geographic datasets
clear
copy https://www2.census.gov/geo/tiger/GENZ2018/shp/cb_2018_us_county_500k.zip ///
     cb_2018_us_county_500k.zip
unzipfile cb_2018_us_county_500k.zip
spshape2dta cb_2018_us_county_500k.shp, saving(usacounties) replace
use usacounties.dta, clear
generate fips = real(GEOID)
save usacounties.dta, replace
// Create the COVID-19 case dataset
clear
import delimited https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_US.csv
rename v83 confirmed
drop if missing(fips)
save covid19_county, replace
// Create the population dataset
clear
import delimited https://www2.census.gov/programs-surveys/popest/datasets/2010-2019/counties/totals/co-est2019-alldata.csv
generate fips = state*1000 + county
save census_popn, replace
// Merge the datasets
clear
use _ID _CX _CY GEOID fips using usacounties.dta
merge 1:1 fips using covid19_county  ///
     , keepusing(province_state combined_key confirmed)
keep if _merge==3
drop _merge
merge 1:1 fips using census_popn  ///
     , keepusing(census2010pop popestimate2019)
keep if _merge==3
drop _merge
drop if inlist(province_state, "Alaska", "Hawaii")
generate confirmed_adj = 100000*(confirmed/popestimate2019)
label var confirmed_adj "Cases per 100,000"
format %16.0fc confirmed_adj
format %16.0fc confirmed popestimate2019
save covid19_adj, replace

// Create the maps
grmap, activate
spset, modify shpfile(usacounties_shp)
grmap confirmed_adj, clnumber(7)
grmap confirmed_adj,                                           ///
     clnumber(8)                                               ///
     clmethod(custom)                                          ///
     clbreaks(0 5 10 15 20 25 50 100 5000)                     ///
     title("Confirmed Cases of COVID-19 in the United States") ///
     subtitle("cases per 100,000 population")





Eso es todo por hoy. ¡Gracias por leernos!