I'm working on a project where I investigate the relationship between PM2.5 (a form of pollution) and ischemic stroke hospitalization (i.e. number of hospitalizations total, collected per day). My model is enclosed below, as is the plot of the results:

gam_model <- gam(Stroke ~ s(PM2.5) +
        factor(districts) +
        ns(date, df= 7) +  
                 data = regional_data,
                 family = poisson(link = "log"))

plot(gam_model, trans = exp, xlab = "PM2.5", ylab = "Relative Risk", rug = TRUE)


Because I want to interpret the model in terms of relative risk (relative to baseline PM2.5 level of zero), I set trans = exp. I don't however understand how to interpret the results, which suggest that the relative risk of being hospitalized for ischemic stroke is less at PM2.5 $\approx 20$ than at 0. Have I made a mistake?

Apologies if this is a dumb question, I'm quite new to all this.

  • 1
    $\begingroup$ Some more info would help - what is PM2.5? What is the response variable "stroke" exactly (number of strokes in some regions?) Could you share the model output and the plot you are interpreting? $\endgroup$
    – N Brouwer
    Commented Apr 7 at 1:00
  • $\begingroup$ Apologies, I've added the plot. PM2.5 is simply a form of pollutant, given here in micrograms per cubic meter of air. The response variable is the number of hospitalizations due to stroke. Both variables are recorded each day. $\endgroup$ Commented Apr 7 at 8:51

3 Answers 3


I would first refit your model with some changes:

  • convert districts to a factor before you fit the model,
  • log transform (or similar) the PM2.5 covariate; very few observations underpin the estimated smooth for this covariate over a large part of the range, and
  • replace ns(date) with a penalized spline using s(date, k = 10) (setting k to whatever you need.

Having done that, run the model diagnostics with mgcv::gam.check() (and look in the R console for the basis dimension checks) and confirm the model is fitting adequately. Run diagnostics with the DHARMa package using randomised quantile residuals for additional checks. Usee a rootogram as an additional diagnostic; my {gratia} package has a rootogram() function for GAMs produced via {mgcv}.

Then you can reproduce your plot (or do so with gratia::draw() if you prefer ggplot graphics) and interpret it (using trans = exp, or tranfun = exp if using draw()) as an incidence rate ratio (it's not a relative risk AFAICT), but that seems to be what you want to achieve, the ratio of the estimated rates.

You are interpreting the plot correctly (though why the baseline would be 0 at PM2.5 = 0 in real data is behind me. People would still be admitted to hospital with stroke even if we had 0 air pollution. Although you don't have data at PM2.5 lower than ~20, the model is suggesting that the rate ratio is even lower at PM2.5 <20.

Once you've log transformed the PM2.5 covariate, things might be a bit easier to interpret.


Gavin and Nicholas have both provided useful commentary on your question, but I would also note that you should probably also include the "district" variable as a random effect. I only assume this because intuitively I'm guessing that you don't actually care about the direct effect of districts, but are more including this variable simply to account for the clustering that may occur from this factor.

One can incorporate this into the model using a variety of specifications:

  • An easy example is one that allows random intercepts, such as s(districts, bs = "re").
  • A random slopes version could be modeled as s(district, PM2.5, bs = "re") if you believe that the association between PM2.5 and RR varies by district.
  • A "smooth" version of the second option could be specified as s(district, PM2.5, bs = "fs").

I agree with Gavin here, a penalized spline of date is preferable as it gives you smoothing regularisation and access to gam.check() for diagnostics. Residual checks are also useful, not only to check distributional assumptions but also to look for unmodelled autocorrelation. As you are working with a time series here, you'll need to specifically check for temporal autocorrelation of residuals (i.e. acf(residuals(gam_model)), or something equivalent assuming your data are arranged by date). Be aware however that gam.check() is sensitive to unmodelled temporal autocorrelation.

Another point to consider is that it is fairly widely accepted that environmental exposures act cumulatively to modulate risk in these contexts. Rather than choosing one specific lag of PM2.5, you may want to look into a distributed lag formulation that allows the conditional effect of PM2.5 to not only be nonlinear, but to change smoothly over increasing lags. {mgcv} makes these models fairly easy to set up and estimate (I have a small blog post on the topic that shows how you can set up the cumulative exposure matrices).


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.