Introduction
This post draws inspiration from the scientific paper titled “Excess of COVID-19 cases and deaths due to fine particulate matter exposure during the 2020 wildfires in the United States.” The study investigates how wildfire-related air pollution may have influenced COVID-19 cases and deaths across 92 counties in the western United States. To quantify this relationship, the authors employed a Bayesian hierarchical zero-inflated negative binomial distributed lag model, linking short-term exposure to fine particulate matter (\(\text{PM}_{2.5}\)) with increased risks of COVID-19 infections and fatalities. This modeling approach allowed the researchers to tackle challenges inherent to the data while extracting meaningful insights. In the following sections, we explore the study’s strengths, limitations, and key takeaways.
The study’s methodology is particularly compelling due to its use of natural spline terms and random effects, which introduce flexibility and help capture complex relationships between the outcomes and explanatory variables. The authors also controlled for major confounders, accounted for inter-county variability, and incorporated both the incubation period and the time from symptom onset to death. In addition, a thorough sensitivity analysis was performed to evaluate the robustness of the findings. Nonetheless, some potential limitations remain. For example, cross-county differences, such as the proportion of residents aged 65 and older, were not included, despite their likely impact on COVID-19 risk. Moreover, studies of this nature should consider variations in containment measures and health system policies. The Oxford COVID-19 Government Response Tracker offers a framework for quantifying such policy stringency across counties.
Modeling Procedure
The original study emphasizes that, due to the COVID-19 incubation period, changes in cases or deaths may not appear immediately after exposure to \(\text{PM}_{2.5}\). Accounting for the time from symptom onset to death is essential for obtaining accurate estimates. With incubation periods ranging from 1 to 14 days, we consider both the average (6 days) and maximum (14 days) durations. Similarly, the time from symptom onset to death is represented by both its average (18 days) and half of its maximum (28 days). Our analysis also adjusts for mobility, temperature, humidity, longitude, and latitude to control for potential confounders. The approach is detailed below.
We adopt a streamlined method for modeling COVID-19 deaths and cases without constructing multiple lagged variables. Specifically, rolling windows of 18 and 28 days are used for deaths, while 6- and 14-day windows are used for cases. Within each window, we calculate the rolling average of \(\text{PM}_{2.5}\), with the outcome defined as the number of deaths or cases on the final day of that window. This provides a single, interpretable effect estimate for \(\text{PM}_{2.5}\) exposure, avoiding the complexity of numerous lagged variables. Once these variables are created, standard zero-inflated or mixed models can be applied. In this study, a Hurdle mixed-effects model with a truncated generalized Poisson distribution is used to appropriately model the count data.
Import the dataset
library(tidyverse)
library(zoo)
library(glmmTMB)
moddat <- read_csv(file = "https://raw.githubusercontent.com/xiaodan-zhou/covid_wildfire/master/data/moddat_Feb2021.csv")
moddat_org <- moddat
moddat %>% glimpse()## Rows: 25,484
## Columns: 31
## $ FIPS <chr> "06001", "06007", "06009", "06011", "06013", "0601~
## $ date <date> 2020-03-15, 2020-03-15, 2020-03-15, 2020-03-15, 2~
## $ cases <dbl> 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 16, 0, 4, 0, 0, 0, 0~
## $ cumu_cases <dbl> 15, 0, 2, 0, 29, 2, 0, 0, 0, 0, 69, 1, 4, 0, 0, 0,~
## $ deaths <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,~
## $ cumu_deaths <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,~
## $ population <dbl> 1671329, 219186, 45905, 21547, 1153526, 999101, 28~
## $ pm25 <dbl> 2.713462, 1.930769, 2.357692, 0.000000, 4.840385, ~
## $ Long <dbl> -122.1176, -121.7002, -120.6803, -121.9989, -122.1~
## $ Lat <dbl> 37.76750, 39.62281, 38.20185, 39.18919, 37.94935, ~
## $ tmmx <dbl> 11.509947, 10.274879, 9.114934, 13.913608, 11.9774~
## $ rmax <dbl> 83.79869, 84.90317, 98.29721, 83.91637, 82.94662, ~
## $ sph <dbl> 0.005371786, 0.005077436, 0.005813920, 0.005089228~
## $ hazardmap <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,~
## $ relative_change_feb <dbl> -0.26310, -0.17022, -0.28117, 0.19178, -0.24632, -~
## $ ratio_travelers <dbl> 0.25689, 0.27878, 0.34894, 0.18750, 0.25855, 0.226~
## $ month <dbl> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,~
## $ day <dbl> 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15~
## $ md <chr> "3-15", "3-15", "3-15", "3-15", "3-15", "3-15", "3~
## $ pm25_history_raw <dbl> 6.123077, 4.434615, 8.430769, 4.053846, 6.653846, ~
## $ County <chr> "Alameda County", "Butte County", "Calaveras Count~
## $ State <chr> "CA", "CA", "CA", "CA", "CA", "CA", "CA", "CA", "C~
## $ StateFIPS <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,~
## $ dayofweek <chr> "Sunday", "Sunday", "Sunday", "Sunday", "Sunday", ~
## $ date_num <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,~
## $ date_str <date> 2020-03-15, 2020-03-15, 2020-03-15, 2020-03-15, 2~
## $ pm25_history <dbl> 6.123077, 4.434615, 8.430769, 4.053846, 6.653846, ~
## $ pm25_raw <dbl> 2.713462, 1.930769, 2.357692, 0.000000, 4.840385, ~
## $ wildfire <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, F~
## $ pm_wildfire <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,~
## $ pm_ambient <dbl> 2.713462, 1.930769, 2.357692, 0.000000, 4.840385, ~moddat <- moddat %>%
dplyr::select(FIPS, deaths, pm25, population, relative_change_feb, Long, Lat, tmmx, rmax) %>%
na.omit()## Modeling COVID-19 deaths
## Average time from symptom onset to death is 18 days (2-8 weeks range)
## Here, I consider two values; 18 and 28 days
## For 18 days time window, set the cum_len to 18 and rerun the code
cum_len <- 28
## Prepare the data
final_data <- c()
for (cid in unique(moddat$FIPS)) {
myDf <- moddat %>% filter(FIPS == cid) %>%
mutate(cum_pm25 = c(rep(NA, (cum_len-1)), rev(rollmean(rev(pm25), k = cum_len))))
final_data <- rbind(final_data, myDf)
}
final_data <- final_data %>% na.omit() %>% mutate(FIPS = factor(FIPS))
## Fit model
fit_mod <- glmmTMB(deaths ~ relative_change_feb + Long + Lat +
tmmx + rmax + cum_pm25 + (1 | FIPS),
offset = log(population),
ziformula = ~ relative_change_feb + Long + Lat +
tmmx + rmax + cum_pm25,
family = truncated_genpois,
data = final_data)
# summary(fit_mod)
# confint(fit_mod)
#################################################
## Modeling COVID-19 cases
moddat <- moddat_org
moddat <- moddat %>%
dplyr::select(FIPS, cases, pm25, population, relative_change_feb, Long, Lat, tmmx, rmax) %>%
na.omit()
## Incubation period is 14 days (median around 6)
## Here, I consider two values; 6 and 14 days
## For 6 days time window, set cum_len to 6 and rerun the code
cum_len <- 14
## Prepare the data
final_data <- c()
for (cid in unique(moddat$FIPS)) {
myDf <- moddat %>% filter(FIPS == cid) %>%
mutate(cum_pm25 = c(rep(NA, (cum_len-1)), rev(rollmean(rev(pm25), k = cum_len))))
final_data <- rbind(final_data, myDf)
}
final_data <- final_data %>% na.omit() %>% mutate(FIPS = factor(FIPS))
## Fit model
fit_mod <- glmmTMB(cases ~ relative_change_feb + Long + Lat +
tmmx + rmax + cum_pm25 + (1 | FIPS),
offset = log(population),
ziformula = ~ relative_change_feb + Long + Lat +
tmmx + rmax + cum_pm25,
family = truncated_genpois,
data = final_data)
# summary(fit_mod)
# confint(fit_mod)library(kableExtra)
# Create a data frame for the first table
my_data <- tibble(
'Time window' = c("18 days", "28 days"),
'Estimate' = c(0.010, 0.007),
'Standard error' = c(0.002, 0.003),
'95% Wald-type confidence interval' = c("(0.005, 0.015)", "(0.002, 0.013)")
)
## Create the first table
my_first_table <- kable(my_data, format = "html",
align = "c",
caption = "Effect of $\\text{PM}_{2.5}$ on COVID-19 deaths under the Hurdle mixed-effects model") %>%
kable_styling(full_width = FALSE)
# Create a data frame for the second table
my_data <- tibble(
'Time window' = c("6 days", "14 days"),
'Estimate' = c(0.002, 0.003),
'Standard error' = c("$3 \\times 10^{-4}$", "$3 \\times 10^{-4}$"),
'95% Wald-type confidence interval' = c("(0.002, 0.003)", "(0.002, 0.003)")
)
## Create the second table
my_second_table <- kable(my_data, format = "html",
align = "c",
caption = "Effect of $\\text{PM}_{2.5}$ on COVID-19 cases under the Hurdle mixed-effects model") %>%
kable_styling(full_width = FALSE)The results are presented in the following tables. From the first and second tables, we can calculate the percentage changes in COVID-19 deaths and COVID-19 cases for the respective time windows. For an 18-day time window, a 10-unit increase in the average \(\text{PM}_{2.5}\) is associated with a 10.52% increase in COVID-19 deaths. For a 28-day time window, the same 10-unit increase corresponds to a 7.25% rise in COVID-19 deaths. Similarly, a 10-unit increase in the average \(\text{PM}_{2.5}\) leads to an increase in COVID-19 cases of 2.02% and 3.05% for time windows of 6 and 14 days, respectively.
| Time window | Estimate | Standard error | 95% Wald-type confidence interval |
|---|---|---|---|
| 18 days | 0.010 | 0.002 | (0.005, 0.015) |
| 28 days | 0.007 | 0.003 | (0.002, 0.013) |
| Time window | Estimate | Standard error | 95% Wald-type confidence interval |
|---|---|---|---|
| 6 days | 0.002 | \(3 \times 10^{-4}\) | (0.002, 0.003) |
| 14 days | 0.003 | \(3 \times 10^{-4}\) | (0.002, 0.003) |
Conclusion
In this post, we examine the relationship between short-term exposure to \(\text{PM}_{2.5}\) and the increased risk of COVID-19 deaths and cases. The results indicate a significant statistical association between short-term \(\text{PM}_{2.5}\) exposure and COVID-19 outcomes. To ensure the robustness of these findings, we recommend conducting a comprehensive sensitivity analysis to assess the effect of different time windows on the results.