Corresponding shiny app: https://michaeldgarber.shinyapps.io/app-wf-int/

1 Introduction

1.1 Background: Do et al. (GeoHealth, 2024)

Do and colleagues have recently estimated spatially varying effects of short-term exposure to wildfire PM2.5 on acute-care hospitalizations (Do et al. 2024) in California. They defined the exposure as a day with wildfire smoke PM2.5 ≥15 μg/m3.

The outcome was the number of daily respiratory emergency department visits and unplanned hospitalizations.

They used a case-crossover design to estimate effects at the level of the zip-code tabulation area (ZCTA) throughout 1,396 ZCTAs in California from 2006-2019. A map of these spatially varying effect estimates, expressed as a rate difference per 100,000 person-days, appears below (adapted from their figure 5).

A positive rate difference (shades of green) indicates that wildfire had an estimated harmful effect (i.e., more) hospitalizations, while a negative value (purple shades) implies that wildfire led to fewer hospitalizations.

In addition, they explored community characteristics possibly driving this spatial heterogeneity. The authors state:

We used meta‐regression to evaluate potential effect modification by community characteristics on the effect of a wildfire smoke day on acute care utilization at the ZCTA level. For each community characteristic, which was selected a priori, we ran a meta‐regression of the pooled ZCTA‐specific rate difference on the community characteristic. To preserve statistical power, we excluded 100 ZCTAs without complete data for 14 community characteristics other than A/C prevalence, and we excluded 274 ZCTAs for meta‐regression of the A/C prevalence. Our estimates are reported as rate difference per interquartile range increase of the community characteristic.

Below, adapted from their figure 6, is a plot of the increase in rate difference per interquartile range (IQR) increase for each socio-demographic community characteristic they considered:

This analysis of effect modification describes how the estimated effect of wildfire smoke on acute-care utilization varies based on variation in the community characteristics. Results show, for example, that the ZCTA-level wildfire-utilization effect was negatively associated with air conditioning (more air conditioning –> lower rate difference) and insurance status (more insurance –> lower rate difference).

Interpreting a specific value from the figure, an increase in A/C proportion from the 25th percentile to the 75th percentile was associated with a 0.239 (95% confidence interval: 0.0672, 0.411) lower rate rate difference in the effect of wildfire smoke on same-day respiratory acute-care utilization.

1.2 Distinction between effect modification vs interaction

It is important to note that just because wildfire’s effects on acute-care utilization vary across these values, it does not necessarilly mean that intervenining to change these variables would change the wildfire-healthcare-utilization effect.

Tyler VanderWeele describes the important distinction between effect modification and interaction on p. 268 of his book, Explanation in Causal Inference (bold added for emphasis here):

If we found that the effect of our primary exposure varied by strata defined by the secondary factor in this way, then we might call  this “effect heterogeneity” or “effect modification.” This might be useful, for example, in decisions about which subpopulations to target in order to maximize the effect of interventions. Provided that we have controlled for confounding of relationship between the primary exposure and the outcome, these estimates of effect modification or effect heterogeneity could be useful even if we have not controlled for confounding of the relationship between the secondary factor and the outcome. (Tyler J. VanderWeele 2015)

What we would not know, however, is whether the effect heterogeneity were due to the secondary factor itself, or something else associated with it. If we have not controlled for confounding for the secondary factor, the secondary factor itself may simply be serving as a proxy for something that is causally relevant for the outcome (Tyler J. VanderWeele and Robins 2007).

VanderWeele continues (bold again added for its relevance here):

If we are interested principally in assessing the effect of the primary exposure within subgroups defined by a secondary factor then simply controlling for confounding for the relationship between the primary exposure and the outcome is sufficient. However, if we want to intervene on the secondary factor in order to change the effect of the primary exposure then we need to control for confounding of the relationships of both factors with the outcome. When we control for confounding for both factors we might refer to this as “causal interaction” in distinction from mere “effect heterogeneity” mentioned above (T. J. VanderWeele 2009).

That is, effect modification can be considered a description of variation in an effect of a primary exposure on the outcome over values over a third variable, whereas assessment of interaction requires the third variable to be considered as a causal factor with adjustment for confounding as needed.

To inform policy, it would be useful assess interaction to know what environmental factors could be intervened upon to mitigate wildfire’s health effects.

Both tree canopy and air conditioning could plausibly attenuate wildfire’s health effects and in that way, interact with wildfire on the outcome. In the case of tree canopy, previous research has shown that green space can reduce air pollution.(Diener and Mudu 2021) And previous research in California has suggested that air conditioning use could mitigate the health effects of wildfires.(Stowell et al. 2025)

In addition, they are socio-environmental pathways that are relatively amenable to intervention. It would be a realistic policy option to encourage more tree planting, for example, or institute policies to encourage more air conditoning use during wildfire events.

However, as both of these variables may be associated with other factors on the causal pathway from wildfire to health, it would be important to control for potential confounding.

1.3 The need for locally developed health-impact assessment

Furthermore, many health-impact assessment studies make a strong transportability assumption–using dose-response functions developed elsewhere and applying it to the local context. There is a need for locally tailored health-impact assessment studies.

1.4 Analysis goals

The primary goal of this analysis is to assess how the distribution of the effect of wildfire onhealthcare utilization effect could change were there interventions changing environmental factors that influence the effect. We do so estimating the effect of tree canopy and air conditioning on the spatially varying rate difference.

A secondary goal, related to training and skill development, is to explore the utility of R Shiny tools for presenting high-dimensional results with which a user can interact. For example, a user could choose to raise the distribution of the effect modifier to the 50th percentile and see how results might change.)

2 Methods

The overall approach is as follows:

  1. Estimate the effect of both tree canopy and air conditioning (separately) on the spatially varying rate difference, controlling for potential confounders that may influence these associations such as urban vs rural status, biome, and socioeconomic characteristics. In these models, we also consider interaction variables, allowing the relationship between these main exposures (tree canopy and AC) and the rate difference to vary across strata of a third variable.
  2. Use these models to simulate the overall burden and spatial distribution of healthcare utilization attributable to wildfire that could be averted by intervening upon these variables.that

2.1 Data exploration

To inform modeling, we explored the distribution of all variables as well as bivariate associations between tree canopy, air conditioning, the rate difference, and other covariates.

July 9, 2025 note: I’m moving much of this over to a separate markdown doc to create a supplemental word doc

~/Library/CloudStorage/Dropbox/Work/CSU-UCSD post-doc/post-doc-research/climate-health-california/climate-health/docs/supplement.Rmd

2.1.1 Outcome: Rate difference per 100,000 person-days

STATE rd_100k_pt_med rd_100k_pt_mean_wt rd_100k_25th rd_100k_75th
CA -0.059 0.179 -0.809 0.809
ruca_cat rd_100k_pt_med rd_100k_pt_mean_wt
(0,3] -0.033 0.196
(3,6] -0.035 -0.164
(6,9] -0.522 0.141
(9,10] -0.732 -0.306

2.1.2 Intervention variables

2.1.2.1 Tree canopy

2.1.2.2 Air conditioning

2.1.2.3 Impervious surface

(exploring as a potential intervention variable)

Impervious surfaces vs. tree canopy

2.1.3 Covariates

2.1.3.1 Proportion above poverty and proportion insured

2.1.3.2 Biome and rural-urban classification

Definitions of rural-urban commuting codes:

  1. Metropolitan area core: primary flow within an urbanized area (UA)
  2. Metropolitan area high commuting: primary flow 30% or more to a UA
  3. Metropolitan area low commuting: primary flow 10% to 30% to a UA
  4. Micropolitan area core: primary flow within an Urban Cluster of 10,000 to 49,999 (large UC)
  5. Micropolitan high commuting: primary flow 30% or more to a large UC
  6. Micropolitan low commuting: primary flow 10% to 30% to a large UC
  7. Small town core: primary flow within an Urban Cluster of 2,500 to 9,999 (small UC)
  8. Small town high commuting: primary flow 30% or more to a small UC
  9. Small town low commuting: primary flow 10% to 30% to a small UC
  10. Rural areas: primary flow to a tract outside a UA or UC

2.1.4 RD x intervention vars and by other covariates

2.1.4.1 RD x intervention variables

2.1.4.1.1 RD x Tree canopy

Here are scatterplots plotting the rate difference against the (square root of) proportion tree canopy:

Observations 1270
Dependent variable rd_100k_quo_pt
Type Linear regression
Est. S.E. t val. p
(Intercept) 0.396 0.095 4.159 0.000
tree_canopy_prop_sqrt -1.051 0.297 -3.545 0.000
Standard errors: MLE

There is a slight negative association between tree canopy and the RD, but it is rather weak.

2.1.4.1.2 RD x Tree canopy x rural-urban

We also stratified this association by urban-rural category:

Observations 1270
Dependent variable rd_100k_quo_pt
Type Linear regression
Est. S.E. t val. p
(Intercept) 0.599 0.113 5.293 0.000
tree_canopy_prop_sqrt -1.831 0.384 -4.767 0.000
ruca_cat(3,6] -0.895 0.367 -2.439 0.015
ruca_cat(6,9] -1.381 0.522 -2.646 0.008
ruca_cat(9,10] 0.663 0.545 1.217 0.224
tree_canopy_prop_sqrt:ruca_cat(3,6] 3.318 0.854 3.884 0.000
tree_canopy_prop_sqrt:ruca_cat(6,9] 2.763 1.402 1.970 0.049
tree_canopy_prop_sqrt:ruca_cat(9,10] -0.861 1.078 -0.799 0.425
Standard errors: MLE

The stratified scatterplot suggests that in the most urban areas and in the least urban areas, higher tree canopy is associated with a lower rate difference, whereas in the other two RUCA categories, the association is in the other direction: higher tree canopy is associated with a higher rate difference (i.e., suggesting a more harmful effect of wildfire on the outcome).

2.1.4.1.3 RD x Tree canopy x biome

We also stratified by biome, which illustrates, as above, that the tree-canopy distribution differs starkly by biome. The RD x tree canopy association appears to be negative in Mediterran Forests, Woodlands & Scrub and in Deserts & Xeric shrublands, but slightly positive in the other two.

Observations 1270
Dependent variable rd_100k_quo_pt
Type Linear regression
Est. S.E. t val. p
(Intercept) 0.774 0.124 6.262 0.000
tree_canopy_prop_sqrt -1.936 0.417 -4.643 0.000
biome_name_freqTemperate Conifer Forests -1.386 0.684 -2.027 0.043
biome_name_freqTemperate Grasslands, Savannas & Shrublands -1.833 0.349 -5.257 0.000
biome_name_freqDeserts & Xeric Shrublands 0.319 0.527 0.605 0.545
tree_canopy_prop_sqrt:biome_name_freqTemperate Conifer Forests 2.712 1.170 2.319 0.021
tree_canopy_prop_sqrt:biome_name_freqTemperate Grasslands, Savannas & Shrublands 4.020 1.219 3.298 0.001
tree_canopy_prop_sqrt:biome_name_freqDeserts & Xeric Shrublands -5.005 3.311 -1.512 0.131
Standard errors: MLE

2.1.4.1.4 RD x A/C

The association with air conditioning is weakly negative.

Observations 1104 (166 missing obs. deleted)
Dependent variable rd_100k_quo_pt
Type Linear regression
Est. S.E. t val. p
(Intercept) 0.326 0.087 3.740 0.000
ac_prop -0.365 0.135 -2.703 0.007
Standard errors: MLE

2.1.4.1.5 RD x A/C x rural-urban

Like with tree canopy, the RD x A/C association varies by rural-urban category. It is more negative in more populous areas, and positive–suggesting air conditioning is associated with a higher difference effect of wildfire on acute-care utilization in less populous areas.

Observations 1104 (166 missing obs. deleted)
Dependent variable rd_100k_quo_pt
Type Linear regression
Est. S.E. t val. p
(Intercept) 0.392 0.092 4.247 0.000
ac_prop -0.451 0.143 -3.157 0.002
ruca_cat(3,6] 0.327 0.344 0.951 0.342
ruca_cat(6,9] -0.963 0.622 -1.550 0.122
ruca_cat(9,10] -1.992 0.470 -4.236 0.000
ac_prop:ruca_cat(3,6] -0.449 0.563 -0.797 0.425
ac_prop:ruca_cat(6,9] 0.550 0.864 0.637 0.524
ac_prop:ruca_cat(9,10] 3.296 0.731 4.507 0.000
Standard errors: MLE

2.1.4.1.6 RD x A/C x biome
Observations 1104 (166 missing obs. deleted)
Dependent variable rd_100k_quo_pt
Type Linear regression
Est. S.E. t val. p
(Intercept) 0.376 0.094 4.005 0.000
ac_prop -0.263 0.154 -1.705 0.088
biome_name_freqTemperate Conifer Forests -0.550 0.277 -1.988 0.047
biome_name_freqTemperate Grasslands, Savannas & Shrublands -1.423 0.553 -2.575 0.010
biome_name_freqDeserts & Xeric Shrublands 0.983 0.980 1.003 0.316
ac_prop:biome_name_freqTemperate Conifer Forests -0.003 0.637 -0.005 0.996
ac_prop:biome_name_freqTemperate Grasslands, Savannas & Shrublands 0.815 0.651 1.253 0.211
ac_prop:biome_name_freqDeserts & Xeric Shrublands -0.633 1.184 -0.534 0.593
Standard errors: MLE

2.1.4.1.7 RD x Impervious surfaces

Sep 23, 2025: Considering impervious surfaces as an alternative effect modifier

Observations 1270
Dependent variable rd_100k_quo_pt
Type Linear regression
Est. S.E. t val. p
(Intercept) -0.289 0.103 -2.800 0.005
imperv_prop 0.837 0.209 4.014 0.000
Standard errors: MLE

There is a slight negative association between tree canopy and the RD, but it is rather weak.

2.1.4.1.8 RD x Impervious surfaces x rural-urban

We also stratified this association by urban-rural category:

Observations 1270
Dependent variable rd_100k_quo_pt
Type Linear regression
Est. S.E. t val. p
(Intercept) -0.496 0.132 -3.753 0.000
imperv_prop 1.208 0.252 4.791 0.000
ruca_cat(3,6] 1.196 0.330 3.626 0.000
ruca_cat(6,9] 0.718 0.625 1.149 0.251
ruca_cat(9,10] -0.241 0.393 -0.612 0.541
imperv_prop:ruca_cat(3,6] -2.894 1.130 -2.562 0.011
imperv_prop:ruca_cat(6,9] -4.128 2.412 -1.711 0.087
imperv_prop:ruca_cat(9,10] 6.408 3.309 1.936 0.053
Standard errors: MLE

2.1.4.1.9 RD x Impervious surfaces x biome
Observations 1270
Dependent variable rd_100k_quo_pt
Type Linear regression
Est. S.E. t val. p
(Intercept) -0.211 0.133 -1.584 0.113
imperv_prop 0.916 0.252 3.638 0.000
biome_name_freqTemperate Conifer Forests -0.338 0.282 -1.197 0.232
biome_name_freqTemperate Grasslands, Savannas & Shrublands -0.297 0.352 -0.843 0.399
biome_name_freqDeserts & Xeric Shrublands 1.814 0.437 4.146 0.000
imperv_prop:biome_name_freqTemperate Conifer Forests 1.686 1.272 1.325 0.185
imperv_prop:biome_name_freqTemperate Grasslands, Savannas & Shrublands -0.884 0.768 -1.150 0.250
imperv_prop:biome_name_freqDeserts & Xeric Shrublands -6.058 1.262 -4.802 0.000
Standard errors: MLE

2.1.4.2 RD x other covariates

2.1.4.2.1 RD x proportion above poverty

2.1.4.2.2 RD x proportion insured

The association between the RD and proportion with insurance is stronger than that between proportion above poverty.

2.1.4.2.3 RD x biome

The RD distribution is in the positive direction in Mediterranean Forests, Woodlands & Scrub and more negative in all the others.

biome_name_freq n_zcta rd_100k_pt_med rd_100k_pt_mean_wt
Mediterranean Forests, Woodlands & Scrub 908 0.154 0.341
Temperate Conifer Forests 101 -0.282 0.043
Temperate Grasslands, Savannas & Shrublands 225 -0.498 -0.444
Deserts & Xeric Shrublands 62 -0.756 -0.534

2.1.4.2.4 RD x urban-rural category

The distribution is more negative in more rural areas, suggesting in those areas, wildfires had a preventive effect on healthcare utilization

ruca_cat n_zcta rd_100k_pt_med rd_100k_pt_mean_wt
(0,3] 1087 -0.033 0.196
(3,6] 99 -0.035 -0.164
(6,9] 51 -0.522 0.141
(9,10] 59 -0.732 -0.306

2.2 Modeling

Considering these empirical trends in the data along with substantive considerations of the plausible causal pathway, we decided upon the following models:

2.2.1 Tree canopy

Control variables:

  • proportion insured

  • rural-urban category

  • biome

Interaction (product terms):

  • rural-urban category

  • biome

glm_rd_100k_tree_conf_ins_ruca_biome_int_ruca_biome =glm(
  rd_100k_quo_pt~tree_canopy_prop_sqrt +insured_prop + ruca_cat +
    biome_name_freq+
    tree_canopy_prop_sqrt*ruca_cat+
    tree_canopy_prop_sqrt*biome_name_freq,
  family = gaussian,
  na.action = "na.exclude", 
  weights = rd_quo_weight,
  data = wf_eff_emm_wide_no_outliers
)

jtools::summ(
  glm_rd_100k_tree_conf_ins_ruca_biome_int_ruca_biome,
  model.fit=T,
  digits=3)
Observations 1270
Dependent variable rd_100k_quo_pt
Type Linear regression
𝛘²(14) 144.392
p 0.000
Pseudo-R² (Cragg-Uhler) 0.081
Pseudo-R² (McFadden) 0.023
AIC 4523.742
BIC 4606.090
Est. S.E. t val. p
(Intercept) 1.390 0.629 2.212 0.027
tree_canopy_prop_sqrt -2.332 0.479 -4.865 0.000
insured_prop -0.593 0.735 -0.807 0.420
ruca_cat(3,6] -0.305 0.376 -0.811 0.418
ruca_cat(6,9] -0.597 0.531 -1.125 0.261
ruca_cat(9,10] 1.352 0.567 2.387 0.017
biome_name_freqTemperate Conifer Forests -1.731 0.717 -2.415 0.016
biome_name_freqTemperate Grasslands, Savannas & Shrublands -1.976 0.359 -5.499 0.000
biome_name_freqDeserts & Xeric Shrublands 0.207 0.545 0.380 0.704
tree_canopy_prop_sqrt:ruca_cat(3,6] 2.292 0.878 2.612 0.009
tree_canopy_prop_sqrt:ruca_cat(6,9] 1.495 1.420 1.053 0.293
tree_canopy_prop_sqrt:ruca_cat(9,10] -2.033 1.131 -1.798 0.072
tree_canopy_prop_sqrt:biome_name_freqTemperate Conifer Forests 2.941 1.230 2.391 0.017
tree_canopy_prop_sqrt:biome_name_freqTemperate Grasslands, Savannas & Shrublands 4.551 1.240 3.669 0.000
tree_canopy_prop_sqrt:biome_name_freqDeserts & Xeric Shrublands -5.054 3.335 -1.515 0.130
Standard errors: MLE

2.2.1.1 Interpreting selected components of the model

Model fit statistics suggest models explain relatively little of the variation. For example:

Pseudo-R² (Cragg-Uhler) = 0.080

The point estimate of -2.150 for tree_canopy_prop_sqrt suggests that, within the referent-group stratum for rural-urban category (the most urban category: RUCA 0-3) and biome (Mediterranean Forests), tree canopy has a negative association with the rate difference.

coefficients(glm_rd_100k_tree_conf_ins_ruca_biome_int_ruca_biome)[2]
## tree_canopy_prop_sqrt 
##             -2.332187

For each interaction term, we can plug in 1 to calculate the association between tree canopy and the outcome within stratum. For example, in ruca_cat(6,9] , the association is still negative, but less so. The effect can be calculated by adding that stratum’s coefficient to the referent group’s coefficient.

coefficients(glm_rd_100k_tree_conf_ins_ruca_biome_int_ruca_biome)[[2]]+
  coefficients(glm_rd_100k_tree_conf_ins_ruca_biome_int_ruca_biome)[[11]]*1
## [1] -0.837296

And, compared with Mediterranean Forests, the association is more strongly negative in Deserts & Xeric shrublands

coefficients(glm_rd_100k_tree_conf_ins_ruca_biome_int_ruca_biome)[[2]]+
  coefficients(glm_rd_100k_tree_conf_ins_ruca_biome_int_ruca_biome)[[15]]*1
## [1] -7.385928

2.2.2 Air conditioning

Control variables:

  • proportion insured

  • rural-urban category

  • biome

Interaction (product terms):

  • rural-urban category

  • biome

glm_rd_100k_e_ac_prop_conf_ins_ruca_biome_int_ruca_biome =glm(
  rd_100k_quo_pt~ac_prop +insured_prop + ruca_cat + biome_name_freq+
    ac_prop*ruca_cat+
    ac_prop*biome_name_freq,
  family = gaussian,
  na.action = "na.exclude", 
  weights = rd_quo_weight,
  data = wf_eff_emm_wide_no_outliers
)

jtools::summ(
  glm_rd_100k_e_ac_prop_conf_ins_ruca_biome_int_ruca_biome ,
  model.fit=T,
  digits=4)
Observations 1104 (166 missing obs. deleted)
Dependent variable rd_100k_quo_pt
Type Linear regression
𝛘²(14) 121.5059
p 0.0000
Pseudo-R² (Cragg-Uhler) 0.0772
Pseudo-R² (McFadden) 0.0215
AIC 3948.0031
BIC 4028.1102
Est. S.E. t val. p
(Intercept) 2.2535 0.6672 3.3774 0.0008
ac_prop -0.2248 0.1557 -1.4431 0.1493
insured_prop -2.1072 0.7429 -2.8365 0.0046
ruca_cat(3,6] 0.2403 0.3660 0.6566 0.5116
ruca_cat(6,9] -0.9904 0.6608 -1.4988 0.1342
ruca_cat(9,10] -2.1581 0.5247 -4.1130 0.0000
biome_name_freqTemperate Conifer Forests 0.1403 0.3406 0.4118 0.6805
biome_name_freqTemperate Grasslands, Savannas & Shrublands -1.5273 0.5517 -2.7684 0.0057
biome_name_freqDeserts & Xeric Shrublands 1.4276 1.0023 1.4243 0.1546
ac_prop:ruca_cat(3,6] -0.0070 0.5838 -0.0120 0.9904
ac_prop:ruca_cat(6,9] 0.7941 0.9021 0.8803 0.3789
ac_prop:ruca_cat(9,10] 3.7917 0.8048 4.7114 0.0000
ac_prop:biome_name_freqTemperate Conifer Forests -1.4898 0.7246 -2.0561 0.0400
ac_prop:biome_name_freqTemperate Grasslands, Savannas & Shrublands 0.8417 0.6468 1.3014 0.1934
ac_prop:biome_name_freqDeserts & Xeric Shrublands -1.2556 1.2017 -1.0449 0.2963
Standard errors: MLE

2.2.3 Other modeling notes

  • Weight more precise estimates more highly. Considering the outcome in the model is itself a modeled result (with a point estimate and a standard error), we use meta-regression to weight observations (ZCTAs) with more precise estimates more highly. Specifically, we weight ZCTAs in each model proportional to the inverse of their standard error.

  • Choice of modeling distribution. We used simple linear regression because the outcome (the rate difference) is roughly normally distributed. Unlike a rate, for example, it is not bounded by zero, so Poisson regression is not needed.

  • Integrating uncertainty. To integrate uncertainty into results, we re-sampled residuals of the model-predicted values from a normal distribution (mean=predicted value, SD=residual standard deviation). We then took the 5th and 95th percentiles of the corresponding values as the 95% confidence intervals.

2.3 Predicting values over several scenarios

2.3.1 Scenario definitions

We then used the model coefficients to predict a new rate difference between wildfire and acute-care utilization over various scenarios.

To define the alternative scenarios, for each intervention variable (tree canopy and A/C), we found quintile values (i.e., 20th percentile, 40th percentile, …, 80th percentile) over ZCTAs within biome. We found quintiles relative to biome because tree canopy is strongly associated with biome, and it may not be realistic, for example, to raise tree canopy in the desert up to the 80th percentile of a temperate forest.

Air conditioning also varies considerably by biome. Presumably it is feasible to raise A/C regardless of biome, but to ensure the intervention is appropriate for the local context, we also defined A/C quintiles relative to biome.

In hypothetical scenarios, in each ZCTA with a value below the corresponding quintile value, we raised its value to that quintile.

And we plugged in the alternative value for the intervention variable into the model, holding all of the other values constant, to predict the alternative rate difference.

To wash out the influence of modeling error, we also used the model to predict the status-quo value for each ZCTA, using each ZCTA’s baseline values.

2.3.2 Calculating the difference between the two rate differences

Difference in rate differences. We then calculated the difference between the alternative rate difference and the baseline (status-quo) predicted rate difference, defined as alternative - status quo.

A negative value in this comparison implies that the alternative RD is lower than the status-quo predicted value, meaning wildfire is estimated to cause fewer healthcare visits under the scenario.

A positive value in this comparison implies the alternative RD is higher than the status-quo value, meaning wildfire is estimated to cause more healthcare visits under the scenario.

Difference in number of attributable cases. We also calculate the difference in the total number of attributable cases per wildfire day under each scenario.

3 Results

3.1 Tree canopy results

Tables summarizing results 1) over all ZCTAs and 2) by RUCA are below.

Target percentile refers to the percentile of the effect modifier to which all ZCTAs with values below that percentile value would be raised under the target scenario. For example, if the target scenario corresponds to the 60th percentile of air conditioning, all ZCTAs below the 60th percentile would have their A/C raised to the 60th percentile.

3.1.1 Summary table by scenario

Interpreting selected values:

By raising tree canopy of those ZCTAs below the 80th percentile for their biome to the 80th percentile…

  • N, ZCTAS intervened: 1,014 ZCTAs would be intervened upon.

  • Diff. in RD per 100k, alt- status quo: 0.19 (95% CI: 0.18, 0.20) acute-care visits attributable to wildfire could be prevented per 100,000 person-days in those ZCTAs. Because of the use of the word ‘prevented’, we flip the sign

  • Diff. in N, cases, alt- status quo: A total of 63.98 (95% CI: 59.94, 67) wildfire-attriburtable acute-care visits could be prevented per wildfire day in those ZCTAs. Again, because of the use of the word ‘prevented’, we flip the sign when writing the results.

Target percentile N, ZCTAS intervened RD per 100k, predicted status-quo RD per 100k, predicted status-quo (95% LL) RD per 100k, predicted status-quo (95% UL) RD per 100k, predicted alternative RD per 100k, predicted alternative (95% LL) RD per 100k, predicted alternative (95% UL) Diff. in RD per 100k, alt- status quo Diff. in RD per 100k, alt- status quo (95% LL) Diff. in RD per 100k, alt- status quo (95% UL) Attributable N, cases, predicted baseline Attributable N, cases, predicted baseline (95% LL) Attributable N, cases, predicted baseline (95% UL) Attributable N, cases, predicted alternative Attributable N, cases, predicted alternative (95% LL) Attributable N, cases, predicted alternative (95% UL) Diff. in N, cases, alt- status quo Diff. in N, cases, alt- status quo (95% LL) Diff. in N, cases, alt- status quo (95% UL)
0.2 255 0.31 0.09 0.52 0.28 0.07 0.50 -0.03 -0.32 0.31 28.16 8.00 47.49 25.69 6.77 46.23 -2.47 -29.43 28.75
0.4 508 0.27 0.11 0.41 0.23 0.08 0.38 -0.05 -0.26 0.16 49.30 20.40 74.38 41.04 14.70 68.00 -8.27 -46.10 29.15
0.6 762 0.23 0.11 0.35 0.15 0.04 0.27 -0.08 -0.24 0.10 61.67 29.42 93.82 40.30 9.88 71.71 -21.37 -63.15 26.03
0.8 1014 0.19 0.09 0.30 -0.01 -0.13 0.09 -0.20 -0.36 -0.06 64.91 31.11 100.69 -3.92 -43.50 32.08 -68.83 -121.45 -19.72
1.0 1266 0.15 0.04 0.24 -0.74 -0.84 -0.65 -0.89 -1.02 -0.74 57.80 14.55 93.87 -286.08 -324.47 -249.29 -343.88 -393.70 -286.56

3.1.2 Summary table by rural-urban category for 80th percentile scenario

Considering the scenario where tree canopy values are raised to the 80th percentile within biome, in the most urban ZCTAs (rural-urban category: (0,3])…

  • N, ZCTAS intervened: 880 ZCTAs would be intervened upon.

  • Diff. in RD per 100k, alt- status quo: 0.21 (95% CI: 0.20, 0.20) acute-care visits attributable to wildfire could be prevented per 100,000 person-days in those ZCTAs. Because of the use of the word ‘prevented’, we flip the sign.

  • Diff. in N, cases, alt- status quo: A total of 68.44 (95% CI: 64.63, 71.97) wildfire-attriburtable acute-care visits could be prevented per wildfire day in those ZCTAs. Again, because of the use of the word ‘prevented’, we flip the sign when writing the results.

Most of the impact is occurring in the most urban ZCTAs. In fact, in the next-most urban category ((3,6]), the intervention is causing more wildfire-attributable visits, as Diff. in N, cases, alt- status quo is positive.

Target percentile Rural-urban category N, ZCTAS intervened RD per 100k, predicted status-quo RD per 100k, predicted status-quo (95% LL) RD per 100k, predicted status-quo (95% UL) RD per 100k, predicted alternative RD per 100k, predicted alternative (95% LL) RD per 100k, predicted alternative (95% UL) Diff. in RD per 100k, alt- status quo Diff. in RD per 100k, alt- status quo (95% LL) Diff. in RD per 100k, alt- status quo (95% UL) Attributable N, cases, predicted baseline Attributable N, cases, predicted baseline (95% LL) Attributable N, cases, predicted baseline (95% UL) Attributable N, cases, predicted alternative Attributable N, cases, predicted alternative (95% LL) Attributable N, cases, predicted alternative (95% UL) Diff. in N, cases, alt- status quo Diff. in N, cases, alt- status quo (95% LL) Diff. in N, cases, alt- status quo (95% UL)
0.8 (0,3] 880 0.20 0.10 0.31 -0.02 -0.14 0.08 -0.22 -0.39 -0.08 65.38 32.00 100.34 -7.72 -45.29 26.82 -73.11 -126.38 -24.55
0.8 (3,6] 58 0.08 -0.40 0.50 0.49 0.02 0.95 0.41 -0.21 1.05 0.84 -4.10 5.06 4.95 0.21 9.68 4.11 -2.15 10.60
0.8 (6,9] 40 -0.48 -0.96 0.04 -0.31 -0.89 0.21 0.18 -0.66 0.83 -1.53 -3.03 0.12 -0.97 -2.82 0.66 0.56 -2.07 2.62
0.8 (9,10] 36 0.22 -0.32 0.81 -0.18 -0.80 0.38 -0.40 -1.25 0.45 0.22 -0.32 0.79 -0.18 -0.78 0.37 -0.39 -1.22 0.44

3.2 Air conditioning results

3.2.1 Summary table by scenario

By raising air conditioning levels ofthose ZCTAs below the 80th percentile for their biome to the 80th percentile…

  • N, ZCTAS intervened: 861 ZCTAs would be intervened upon.

  • Diff. in RD per 100k, alt- status quo: 0.05 (95% CI: 0.03, 0.06) acute-care visits attributable to wildfire could be prevented per 100,000 person-days in those ZCTAs. Because of the use of the word ‘prevented’, we flip the sign

  • Diff. in N, cases, alt- status quo: A total of 13.17 (95% CI: 9.41, 16.88) wildfire-attriburtable acute-care visits could be prevented per wildfire day in those ZCTAs. Again, because of the use of the word ‘prevented’, we flip the sign when writing the results.

Target percentile N, ZCTAS intervened RD per 100k, predicted status-quo RD per 100k, predicted status-quo (95% LL) RD per 100k, predicted status-quo (95% UL) RD per 100k, predicted alternative RD per 100k, predicted alternative (95% LL) RD per 100k, predicted alternative (95% UL) Diff. in RD per 100k, alt- status quo Diff. in RD per 100k, alt- status quo (95% LL) Diff. in RD per 100k, alt- status quo (95% UL) Attributable N, cases, predicted baseline Attributable N, cases, predicted baseline (95% LL) Attributable N, cases, predicted baseline (95% UL) Attributable N, cases, predicted alternative Attributable N, cases, predicted alternative (95% LL) Attributable N, cases, predicted alternative (95% UL) Diff. in N, cases, alt- status quo Diff. in N, cases, alt- status quo (95% LL) Diff. in N, cases, alt- status quo (95% UL)
0.2 205 0.28 0.02 0.55 0.27 0.02 0.51 -0.01 -0.38 0.34 16.60 1.13 32.29 15.72 1.17 29.91 -0.88 -22.37 19.74
0.4 434 0.26 0.10 0.40 0.22 0.05 0.41 -0.04 -0.29 0.22 34.89 13.04 54.61 29.84 7.11 54.73 -5.05 -39.06 29.41
0.6 654 0.21 0.07 0.34 0.15 0.01 0.28 -0.06 -0.25 0.13 44.35 13.84 71.59 30.80 2.59 58.99 -13.55 -53.08 26.55
0.8 861 0.20 0.08 0.31 0.11 0.00 0.22 -0.08 -0.25 0.09 55.42 23.02 86.69 32.05 0.18 62.68 -23.37 -70.67 24.50
1.0 971 0.18 0.07 0.29 0.07 -0.03 0.18 -0.11 -0.26 0.03 58.74 23.96 92.63 22.84 -9.69 58.56 -35.90 -84.98 9.29

3.2.2 Summary table by rural-urban category for 80th percentile scenario

  • The most urban areas have the strongest estimated effect in terms of the number of total attributable cases: Diff. in N, cases, alt- status quo: -12.26 (95% CI: -15.80, -8.53)

  • On a per-population basis, effects are strongest in the next-most urban stratum: Diff. in RD per 100k, alt- status quo: -0.18 (95% CI: -0.29, -0.06)

Target percentile Rural-urban category N, ZCTAS intervened RD per 100k, predicted status-quo RD per 100k, predicted status-quo (95% LL) RD per 100k, predicted status-quo (95% UL) RD per 100k, predicted alternative RD per 100k, predicted alternative (95% LL) RD per 100k, predicted alternative (95% UL) Diff. in RD per 100k, alt- status quo Diff. in RD per 100k, alt- status quo (95% LL) Diff. in RD per 100k, alt- status quo (95% UL) Attributable N, cases, predicted baseline Attributable N, cases, predicted baseline (95% LL) Attributable N, cases, predicted baseline (95% UL) Attributable N, cases, predicted alternative Attributable N, cases, predicted alternative (95% LL) Attributable N, cases, predicted alternative (95% UL) Diff. in N, cases, alt- status quo Diff. in N, cases, alt- status quo (95% LL) Diff. in N, cases, alt- status quo (95% UL)
0.8 (0,3] 738 0.20 0.08 0.32 0.12 0.00 0.23 -0.08 -0.25 0.09 53.83 21.51 84.43 31.24 -0.31 61.04 -22.59 -67.33 24.50
0.8 (3,6] 69 0.29 -0.19 0.74 0.13 -0.30 0.53 -0.16 -0.83 0.44 2.97 -1.91 7.43 1.30 -3.07 5.40 -1.67 -8.42 4.44
0.8 (6,9] 29 -0.57 -1.17 0.07 -0.52 -1.09 0.06 0.05 -0.85 0.84 -1.41 -2.89 0.17 -1.28 -2.69 0.15 0.13 -2.09 2.08
0.8 (9,10] 25 0.03 -0.76 0.80 0.85 0.10 1.57 0.82 -0.30 2.00 0.03 -0.69 0.73 0.78 0.09 1.44 0.76 -0.28 1.84

3.3 Comparing interventions: tree canopy vs A/C

At the 80th percentile for each, the tree canopy scenario has about four times the effect strength as air conditioning regarding its attenuation of the wildfire-healthcare-utilization effect.

rd_100k_diff_pred_pt_tree=tree_summary_overall_for_kable %>% 
  filter(target_percentile==.8) %>% 
  dplyr::select(rd_100k_diff_pred_pt)  %>% 
  pull()

rd_100k_diff_pred_pt_ac=ac_summary_overall_for_kable %>% 
  filter(target_percentile==.8) %>% 
  dplyr::select(rd_100k_diff_pred_pt) %>% 
  pull()

rd_100k_diff_pred_pt_tree/rd_100k_diff_pred_pt_ac
## [1] 2.442241

3.4 Selected maps of results

3.4.1 Tree canopy

Map of difference in RD per population by ZCTA for the 80th percentile scenario.

Negative values denote the intervention lowers the effect of wildfire on acute-care visits.

3.4.2 Air conditioning

Same scenario

4 Discussion points

  • Advantage over traditional health-impact assessment. In most health-impact assessment studies, the exposure-response is drawn from the literature and applied to the local setting. This requires a very strong transportability assumption that the exposure-response estimate, which is often estimated in different settings, is applicable to the local context. By contrast, we use a locally-derived estimate, which relaxes this transportability assumption.

  • Control potential confounding.

    • We take a two-stage approach:

      • first, estimate rate difference in each area (Do and colleagues, GeoHealth, 2024),

      • second, model influence of environmental factors on that rate difference)

    • In this second stage, we consider potential confounding variables in n the effect of the environmental variables and the spatially varying rate difference.

  • Poor model fit. Nevertheless, in this case, the model fit is quite poor, as can be seen visually by the bivariate scatterplots and by the model-fit statistics. Results ought to be interpreted with large uncertainty in mind. As currently calculated (using predicted standard errors), confidence intervals may suggest more precision in results than is warranted. We may consider calculating intervals a different way, calculating residuals for each observation’s predicted value (predicted-truth) and then calculating confidence intervals by re-sampling residuals.

  • Why useful. Keeping this poor model fit in mind, results can be used to inform targeted environmental interventions to mitigate the negative health effects of wildfire.

  • Consideration: consider limiting analyses to those ZCTAs where the proposed intervention would have a beneficial impact (based on values of interaction coefficients)

References

Diener, Arnt, and Pierpaolo Mudu. 2021. “How Can Vegetation Protect Us from Air Pollution? A Critical Review on Green Spaces’ Mitigation Abilities for Air-Borne Particles from a Public Health Perspective - with Implications for Urban Planning.” Science of The Total Environment 796 (November): 148605. https://doi.org/10.1016/j.scitotenv.2021.148605.
Do, V., C. Chen, T. Benmarhnia, and J. A. Casey. 2024. “Spatial Heterogeneity of the Respiratory Health Impacts of Wildfire Smoke PM 2.5 in California.” GeoHealth 8 (4): e2023GH000997. https://doi.org/10.1029/2023GH000997.
Stowell, Jennifer D, Ian Sue Wing, Yasmin Romitti, Patrick L Kinney, and Gregory A Wellenius. 2025. “Emergency Department Visits in California Associated with Wildfire PM2.5 : Differing Risk Across Individuals and Communities.” Environmental Research: Health 3 (1): 015002. https://doi.org/10.1088/2752-5309/ad976d.
VanderWeele, T J. 2009. “On the Distinction Between Interaction and Effect Modification.” Epidemiology (Cambridge, Mass.) 20 (6): 863.
VanderWeele, Tyler J. 2015. “An Introduction to Interaction Analysis.” In, 249–85. New York, NY: Oxford University Press.
VanderWeele, Tyler J, and James M Robins. 2007. “Four Types of Effect Modification: A Classification Based on Directed Acyclic Graphs.” Epidemiology 18 (5): 561–68. https://doi.org/10.1097/EDE.0b013e318127181b.