Get census data

To look at the trends on a state level, we’ll need to adjust our dataset to match the Census/ACS format instead of the CDC/WONDER table format. This will mostly consist of adding together the data from NYC and NY state and joining the data with GEOID’s from the census.

First we’ll add together the reporting areas NEW YORK CITY and NEW YORK. Then we’ll add some identifers like abbreviated state names (e.g. TX) and the census divisions.

shigella.states <- shigella %>%
  mutate(reporting_area = recode(reporting_area,   # Make NY = NYC + NY
                                 "NEW YORK CITY"="NEW YORK")) %>%
  filter(reporting_area!="UNITED STATES") %>%      # ignore total counts
  group_by(reporting_area, date, year, mmwr_week, pathogen) %>%
  summarise(
    weekly_cases = sum(weekly_cases, na.rm = T),
    cum_thisYear = sum(cum_thisYear, na.rm = T),
    cum_lastYear = sum(cum_lastYear, na.rm = T)
  ) %>%
  ungroup() %>%
  rename(state_name=reporting_area) %>%
  mutate(state_name = str_to_title(state_name)) %>%
  
  # Add abbreviated state names & census divisions
  left_join(tibble(
    state_name = c(datasets::state.name, "District Of Columbia"),
    state_abb  = c(datasets::state.abb, "DC"),
    state_div  = c(as.character(datasets::state.division), "South Atlantic")
    ))

For this portion, I’ll use the tidycensus package to call the 5-year summary ACS estimates of population per state via table B01003. This also gets the shapefiles for each state, so we can make maps later

require(tidycensus)
require(sf)
require(leaflet)

population_sf <- get_acs(geography = "state", year=2017, cache_table=T,
        variables = c(population="B01003_001"), # table = "B01001",
        shift_geo=T, geometry = TRUE, survey="acs5")

population_sf$NAME <- str_to_title(population_sf$NAME)

# Save data
population_sf %>% write_rds("geoData/population_sf.RDS")

Now we’ll add the population data to our table of state data, so we can calculate per capita cases later on

shigella.states <- shigella.states %>%
  left_join(tibble(
    GEOID      = population_sf$GEOID,
    state_name = population_sf$NAME,
    Pop        = population_sf$estimate
  )) %>%
  select(state_abb, GEOID, date:cum_lastYear, everything(), Pop)

State trends

To begin, we’ll make a table of cases per year by state. We’ll also add two columns:

• annual_avg_cases: The average number of annual cases per year
• annual_avg_IR: The average incidence rate (IR) across the years (per 100,000)

state_avg <- shigella.states %>%
  group_by(state_abb, state_div, state_name, Pop, year) %>%
  summarise(annual_cases = sum(weekly_cases, na.rm = T)) %>%
  mutate(
    annual_avg_cases = mean(annual_cases),
    annual_avg_IR    = (annual_avg_cases/Pop)*100000
  ) %>%
  spread(year, annual_cases) %>%
  arrange(desc(annual_avg_IR)) %>% 
  ungroup() %>%
  mutate(
    rankTotal = (max(row_number(annual_avg_cases)) - row_number(annual_avg_cases) + 1),
    rankIR    = max(row_number(annual_avg_IR)) - row_number(annual_avg_IR) + 1) %>%
  select(state_abb, annual_avg_cases:`2019`, rankTotal, rankIR, everything())

Now we’ll plot the individual graphs for some states. Since showing all 50 states would be overkill, let’s look at the top 15 states by their (historical) incidence rates. The Y-axis here is still in the raw number of cumulative cases, but is scaled individually for each plot.

Figure 3: YTD cases for top states

Around week 20, there’s an uptick in cases nationally, and in Louisiana, Texas, Florida, NYC (not shown), and California (also not shown). Arkansas also looks like it might have had a bit of an uptick around week 15.

To see this better, we’ll look at the trend in weekly cases over the past few years (similar to the national graph shown earlier), but instead broken down by state and controlled for population.

Figure 4: Weekly cases per million population

I had to exclude Oklahoma from the graph above, because they had a massive spike in cases during 2015 and it was throwing the scale off. Interestingly, this was around the same time that Texas saw it’s spike, as shown below in Fig 5. They actually line up quite well when you transform the y-axis by a square root

Figure 5: TX & OK cases over time

In another report, I looked at the distribution of where Shigella cases happened in Texas. It’s not related to our project, but this spike in Oklahoma cases happened around the same time there was an apparent outbreak in the Texas panhandle (see DSHS)

Division trends

Texas, Louisiana, and Arkansas are all in the same geographical area (the West South Central division of the US), so let’s just focus on this region. The figure below shows the weekly cases for each of these four states

Figure 6: West South Central division cases by year

Note that the scales differ on each one of the graphs in the figure above. It’s clear that there appears to be an increase in cases during 2019, but it’s hard to tell exactly what’s going on given some of these states have a low base rate and differing population sizes. I’ll try to tackle these problems below.

Figure 7 requires a bit more explanation, as it’s calculating the fold increase in weekly incidence rates This is essentially the incidence rate standardized to each states’ historical incidence rates (similar to standardized incidence/mortality ratios). For clarity I’ve put an explanation of the math below, but you can just skip to figure 7 if you like.

To explain how I got the figure 7 data, take a hypothetical state with the example data below:

This state has a population of 2, and I’ve already calculated the incidence rate (IR) on that table. Taking the median across all rows we get a historical median IR to be 2. We then take the difference from the historical median (IR - median_IR) to find each row’s change from the historical median. Finally, the fold change is the Diff_from_median / median_IR, which is what I have on the y-axis. This gets us this table:

So that’s how I derived the data for the next figure. Because the denominator of the y-axis can’t be zero, I had to exclude all states who had a median number of zero cases per week since 2015. I’ve also only included states who have had less than 5 weeks of being in the “top 10” states nationally during 2019, because we’re interested in outbreaks that are sustained over time

Figure 7: Fold change in IR (2019)