6
$\begingroup$

I am trying to find the correlation between county population density and total COVID-19 cases in that county. This has been running for four hours so far. What's going wrong?

data = ResourceData["Epidemic Data for Novel Coronavirus COVID-19", 
   "USCounties"];

states = EntityList[
   EntityClass["AdministrativeDivision", "AllUSStatesPlusDC"]];

usaCountyCases[s_, c_] := 
 data[Select[MatchQ[Interpreter["USState"][s], #State] &]][
     Select[MatchQ[
        Interpreter["USCounty"][c <> " county, " <> s], #County] &]][
    All, #ConfirmedCases["LastValue"] &] // Normal // First

Clear[counties]

Length[counties = 
   Flatten[EntityClass["AdministrativeDivision", 
      "ContinentalUSStates"]["Subdivisions"]]];

AppendTo[counties, 
  Entity["AdministrativeDivision", {"DistrictOfColumbia", 
    "DistrictOfColumbia", "UnitedStates"}]];

g[k_] := usaCountyCases[counties[[k]][[2]][[2]], 
   StringDelete[counties[[k]][[2]][[1]], "County"]]*1.

f[k_] := QuantityMagnitude[
   Quantity[
    AdministrativeDivisionData[counties[[k]], "PopulationDensity"], (
    "People")/("Miles")^2]]*1.

t = Table[{f[k], g[k]}, {k, Length[counties]}];

Correlation[t[[All, 1]], t[[All, 2]]]
Improve this question
$\endgroup$
2
  • $\begingroup$ How are you defining usaCountyCases? $\endgroup$
    – Michael Seifert
    Commented Jun 22, 2020 at 20:50
  • $\begingroup$ Please check edit. Thank you. $\endgroup$
    – Veritas Lux
    Commented Jun 22, 2020 at 21:20

2 Answers 2

Reset to default
8
$\begingroup$

You may want to look at Bob Sandheinrich's post on the Wolfram Community Website US Counties COVID-19 confirmed cases by population density timelines for ideas on how to deal with population density. It is not super fast, but the timings are not on the order of multiple hours.

Alternatively, you could try to build your own dataset from data tables from freely available websites and bypass some of the overhead of entity framework. If you look at csv file from the NYT GitHub repository, then you will notice that counties are identified by a "fips" code. You can use this fips code join tables available from the Census Bureau.

To request population estimates from census.gov, you will need to request your own personal API key from https://api.census.gov/data/key_signup.html

There are a few missing data that you need to decide how to handle. The following workflow builds associations between fip codes and population and county area.

apiKey = "<---Your API Key--->";
census = Import[
   "https://api.census.gov/data/2018/pep/population?get=GEONAME,POP&\
for=county:*&key=" <> apiKey, "JSON"];
(* Convert POP from string to int *)
census[[2 ;; -1, 2]] = ToExpression /@ census[[2 ;; -1, 2]];
(* convert state and county code to fips *)
(* Create fips String *)
fips = Join[{"fips"}, (#1 <> #2) & @@@ census[[2 ;; -1, {3, 4}]]];
(* Add fips Column to census *)
census = Join[census, Transpose@{fips}, 2];
(* Function to convert list to association *)
convertListToAssociation = AssociationThread[census[[1]], #] &;
(* Build dataset *)
dscensus = 
  Dataset@Map[convertListToAssociation, census[[All]][[2 ;; -1]]];
(* Create fake fip for New York City *)
dscensus = 
  AppendTo[dscensus, <|"GEONAME" -> "New York City, New York", 
    "POP" -> 8622698, "state" -> "36", "county" -> "999", 
    "fips" -> "36999"|>];
(* display fake fip *)
dscensus[Select[#fips == "36999" &], All]
(* create fips\[Rule]POP association *)
fpAssoc = 
  AssociationThread[
   Evaluate[Normal@(dscensus[All, "fips"] -> dscensus[All, "POP"])]];
(* Get area fips spreadsheet from census bureau *)
dsarea = Import[
   "https://www2.census.gov/library/publications/2011/compendia/usa-\
counties/excel/LND01.xls", {"Dataset", 1}, "HeaderLines" -> 1];
(* create fips\[Rule]Area association *)
fpAreaAssoc = 
  AssociationThread[
   Evaluate@Normal[dsarea[All, "STCOU"] -> dsarea[All, "LND010190D"]]];
(* Add fake fips area for NYC *)
fpAreaAssoc = 
  AppendTo[fpAreaAssoc, <|
    "36999" -> 
     QuantityMagnitude@
      Entity["City", {"NewYork", "NewYork", "UnitedStates"}][
       EntityProperty["City", "Area"]]|>];
(* Add missing fips areas *)
fpAreaAssoc["08014"] = 
  QuantityMagnitude@
   Entity["City", {"Broomfield", "Colorado", "UnitedStates"}][
    EntityProperty["City", "Area"]];
fpAreaAssoc["02275"] = 
  QuantityMagnitude@
   Entity["City", {"Wrangell", "Alaska", "UnitedStates"}][
    EntityProperty["City", "Area"]];
(* Some Alaska Counties that had zero areas *)
fpAreaAssoc["02068"] = 12777.076215885369`;
fpAreaAssoc["02230"] = 464.03372208688813`;
fpAreaAssoc["02282"] = 9462.681804288492`;

Update

We can extend the dataset to include fips codes and population data with the following workflow:

data = ResourceData["Epidemic Data for Novel Coronavirus COVID-19", 
   "USCounties"];
(* Find fips codes by administrative division *)
fips2 = #[EntityProperty["AdministrativeDivision", "FIPSCode"]] & /@ 
   data[All, "County"];
countyFips = 
  AssociationThread[Evaluate@Normal[data[All, "County"] -> fips2]];
countyFips[Entity["City", {"NewYork", "NewYork", "UnitedStates"}]] = 
  "36999";
(* Append fips code to dataset *)
data = data[All, <|#, "fips" -> countyFips[#County]|> &];
(* Append Population, Area, and Population Density to dataset *)
data = data[
  All, <|#, "POP" -> fpAssoc[#fips], "AREA" -> fpAreaAssoc[#fips], 
    "POPDENSITY" -> fpAssoc[#fips]/fpAreaAssoc[#fips] // N|> &]

Extended dataset

Now, you have a dataset the includes the population density. It took about a minute to create the new dataset on my machine.

Update 2: 8 Second Workflow to Calculate Correlation

I discovered that generating countyFips could vary from 60 s to 2000 s depending on the time of day so I stored that association in a pastebin.

The following workflow to calculate a correlation takes about 8 s to execute on my machine for a fresh Mathematica session.

AbsoluteTiming[
 apiKey = "<---Your API Key--->";
 census = 
  Import["https://api.census.gov/data/2018/pep/population?get=GEONAME,\
POP&for=county:*&key=" <> apiKey, "JSON"];
 (* Convert POP from string to int *)
 census[[2 ;; -1, 2]] = ToExpression /@ census[[2 ;; -1, 2]];
 (* convert state and county code to fips *)
 (* Create fips String *)
 fips = Join[{"fips"}, (#1 <> #2) & @@@ census[[2 ;; -1, {3, 4}]]];
 (* Add fips Column to census *)
 census = Join[census, Transpose@{fips}, 2];
 (* Function to convert list to association *)
 convertListToAssociation = AssociationThread[census[[1]], #] &;
 (* Build dataset *)
 dscensus = 
  Dataset@Map[convertListToAssociation, census[[All]][[2 ;; -1]]];
 (* Create fake fip for New York City *)
 dscensus = 
  AppendTo[dscensus, <|"GEONAME" -> "New York City, New York", 
    "POP" -> 8622698, "state" -> "36", "county" -> "999", 
    "fips" -> "36999"|>];
 (* display fake fip *)
 dscensus[Select[#fips == "36999" &], All];
 (* create fips->POP association *)
 fpAssoc = 
  AssociationThread[
   Evaluate[Normal@(dscensus[All, "fips"] -> dscensus[All, "POP"])]];
 (* Get area fips spreadsheet from census bureau *)
 dsarea = 
  Import["https://www2.census.gov/library/publications/2011/compendia/\
usa-counties/excel/LND01.xls", {"Dataset", 1}, "HeaderLines" -> 1];
 (* create fips->Area association *)
 fpAreaAssoc = 
  AssociationThread[
   Evaluate@Normal[dsarea[All, "STCOU"] -> dsarea[All, "LND010190D"]]];
 (* Add fake fips area for NYC *)
 fpAreaAssoc = 
  AppendTo[fpAreaAssoc, <|
    "36999" -> 
     QuantityMagnitude@
      Entity["City", {"NewYork", "NewYork", "UnitedStates"}][
       EntityProperty["City", "Area"]]|>];
 (* Add missing fips areas *)
 fpAreaAssoc["08014"] = 
  QuantityMagnitude@
   Entity["City", {"Broomfield", "Colorado", "UnitedStates"}]["Area"];
 fpAreaAssoc["02275"] = 
  QuantityMagnitude@
   Entity["City", {"Wrangell", "Alaska", "UnitedStates"}]["Area"];
 (* Some Alaska Counties that had zero areas *)
 fpAreaAssoc["02068"] = 
  QuantityMagnitude@
   Entity["AdministrativeDivision", {"DenaliBorough", "Alaska", 
      "UnitedStates"}]["Area"];
 fpAreaAssoc["02230"] = 
  QuantityMagnitude@
   Entity["AdministrativeDivision", {"SkagwayMunicipality", "Alaska", 
      "UnitedStates"}]["Area"];
 fpAreaAssoc["02282"] = 
  QuantityMagnitude@
   Entity["AdministrativeDivision", {"YakutatCityAndBorough", 
      "Alaska", "UnitedStates"}]["Area"];
 fpAreaAssoc["02158"] = 
  QuantityMagnitude@
   Entity["AdministrativeDivision", {"WadeHamptonCensusArea", 
      "Alaska", "UnitedStates"}]["Area"];
 (* A missing South Dakota county *)
 fpAreaAssoc["46102"] = 
  QuantityMagnitude@
   Entity["AdministrativeDivision", {"ShannonCounty", "SouthDakota", 
      "UnitedStates"}]["Area"];
 (* Uncomment if update is necessary *)
 (*ResourceUpdate["Epidemic Data for Novel Coronavirus COVID-19"];*)
 {tdata, data} = 
  AbsoluteTiming@
   ResourceData["Epidemic Data for Novel Coronavirus COVID-19", 
    "USCounties"];
 (* I created a pastebin of the county-> fips association for \
future reference *)
 {tfips, countyFips} = 
  AbsoluteTiming[
   Uncompress[Import["https://pastebin.com/raw/syLxP3rY", "Text"]]];
 data = data[All, <|#, "fips" -> countyFips[#County]|> &];
 data = data[
   All, <|#, "POP" -> fpAssoc[#fips], "AREA" -> fpAreaAssoc[#fips], 
     "POPDENSITY" -> fpAssoc[#fips]/fpAreaAssoc[#fips] // N|> &];
 {Correlation[
   Normal@Query[Map[Identity], #ConfirmedCases["LastValue"] &][data], 
   Normal@Query[Map[Identity], #POPDENSITY &][data]],
  Correlation[
   Normal@Query[Map[Identity], #Deaths["LastValue"] &][data], 
   Normal@Query[Map[Identity], #POPDENSITY &][data]]}
 ]
(* {7.96627, {0.69391, 0.708733}} *)

Original Answer with NYT database

Now, we are ready to import the NYT county data from GitHub. The dataset has some imperfections. New York City is missing because it is composed of 5 counties, so I created a fake fips for NYC. Also, there are some "Unknown" and Kansas City seems to be missing a fips, so I excluded those data.

(* Import New York Times Database From Github Repo *)
url = "https://raw.githubusercontent.com/nytimes/covid-19-data/master/\
us-counties.csv";
dsnyt = SemanticImport[
   url, {"String", "String", "String", "String", "Integer", 
    "Integer"}];
(* Exclude Unknown From Database *)
dsnyt = dsnyt[Select[(#county != "Unknown") &], All];
(* Exclude KC From Database Since it lacks a fips code *)
dsnyt = dsnyt[Select[(#county != "Kansas City") &]];
(* Show the remaining missing all from NYC *)
dsnyt[Select[MissingQ[#fips] &], All];
dsnyt[Select[#fips == "" &], All];
(* Replace all missing fips for NYC with fake 36999 fips code *)
dsnyt = dsnyt[All, All] /. "" -> "36999";
(* Add population, area, and population density to NYT dataset*)
dsnyt = dsnyt[
  All, <|#, "POP" -> fpAssoc[#fips], "AREA" -> fpAreaAssoc[#fips], 
    "POPDENSITY" -> fpAssoc[#fips]/fpAreaAssoc[#fips] // N|> &]

New York Times Data With Population Data Added

It took about 20 s to populate the final table with my machine.

Improve this answer
$\endgroup$
2
$\begingroup$

The simplest answer may just be that the calls to retrieve the case data are taking a long time. Here's what I see from my machine:

Timing[f[1500]]
Timing[g[1500]]

{0.136274, 7.71462}
{5.78885, 2.}

Response times vary, but they seem to be in the neighborhood of 4–6 seconds per data point retrieved. Since there are 3108 counties in your data set, a bit of simple math indicates that it would take in the ballpark of 4–5 hours to retrieve all of the data you have requested.

Whether there is any way to speed up these retrievals I do not know.

Improve this answer
$\endgroup$
1
  • $\begingroup$ It just finished in 20,300 seconds. I am getting some Missing data so am working on working with that. Thanks. $\endgroup$
    – Veritas Lux
    Commented Jun 22, 2020 at 23:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.