all, after reviewing this excellent tutorial, I'm trying to expand out my use of GAMs for space and time. I'm trying to use the Pennsylvania smoking and cancer data, but am running into a problem when I try and work with multiple values of predictor per county. Let's consider the data.


penn_sf <- st_as_sf(pennLC$spatial.polygon) %>%
  mutate(county = unique(pennLC$data$county)) %>%
  left_join(pennLC$data) %>%
  mutate(rate = cases/population*1000) %>%
  st_transform("+proj=utm +zone=17 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")

This is great, but for a MRF, we'll want a neighborhood matrix of the counties - and this data has many rows per county

> penn_sf %>% as_tibble %>% dplyr::select(-geometry)
# A tibble: 1,072 x 8
   county cases population race  gender age        rate    ID
   <fct>  <int>      <int> <fct> <fct>  <fct>     <dbl> <int>
 1 adams      0       1492 o     f      Under.40  0         1
 2 adams      0        365 o     f      40.59     0         1
 3 adams      1         68 o     f      60.69    14.7       1
 4 adams      0         73 o     f      70+       0         1
 5 adams      0      23351 w     f      Under.40  0         1
 6 adams      5      12136 w     f      40.59     0.412     1
 7 adams      5       3609 w     f      60.69     1.39      1
 8 adams     18       5411 w     f      70+       3.33      1
 9 adams      0       1697 o     m      Under.40  0         1
10 adams      0        387 o     m      40.59     0         1
# … with 1,062 more rows

which is fine - I want to look at race, gender, and age all as predictors. So, let's create a neighborhood matrix for the state, first aggregating at the county level.

penn <- penn_sf %>%
 group_by(county) %>%
  summarize(cases = sum(cases), population = sum(population))


nb <- poly2nb(penn,  row.names = as.character(penn$county))
names(nb) <- penn$county

Let's see it!

plot(nb, coords =coords, col = "red", add = TRUE)

enter image description here

However, when I then fit the appropriate model...

gam_mrf <- gam(cases ~ s(county, bs = 'mrf', 
                         xt = list(nb = nb)) +
                 race + gender + age + offset(population), # define MRF smooth
               data = penn_sf,
               family = poisson) 

I get the following error

Error in gam.fit3(x = X, y = y, sp = L %*% lsp + lsp0, Eb = Eb, UrS = UrS,  : 
  inner loop 1; can't correct step size
In addition: Warning message:
Step size truncated due to divergence 

I've tried a few variants on this based on things I've seen around - using numerics/integers and creating an ID column which is the integer transform of the county factor in penn_sf, a neighborhood matrix created from penn_sf itself, and others - but all yield either this or another error.

1) What's the issue here, and any thoughts on how to fix?

2) Is it correct to use the aggregated polygons, rather than the original one row per data point polygons? It "feels" correct, but that's never a good way to decide things!

BTW: My goal is to then use this model to get fits, sum them up, and see if the observations match total predictions (in addition to looking at it by group) as well as trying a leave-one-out and seeing how well we can predict data for a missing county.


The issue seems to be a badly specified model; you want the offset term to be on the same scale as the link function (here the log scale). Use offset(log(population)).

What the error is saying is that the optimisation failed because it couldn't correct the step size — that is the size of the step to take down the negative log likelihood function. The steps to fix the step size adaptively during fitting appear top have failed as the log-likelihood seems to diverging no matter how the optimisation proceeds. I can certainly see how this is happening if you have input populations on the raw scale like you have here.

Also, make sure you are fitting with REML; use method = "REML" in the gam() call.

If that doesn't fix the problem, remove the factor covariates and see if you can estimate the MRF model only.


Gavin's answer above starts to get at it - the big thing here is that this is a binomial model, and I fit a poisson! With the wrong offset! A binomial GAM works great here instead -

#The Binomial GAM
gam_mrf <- gam(cases/population ~ s(county, bs = 'mrf', 
                         xt = list(nb = nb)) +
                 race + gender + age, # define MRF smooth
               data = penn_sf,
               method = "REML",
               family = binomial(),
               weights = population) 

And that's with the matrix created above. Further, you can see how well it fits in the summary, or, we can aggregate it and see how well it reproduces trends seen in the data.

Here's aggregating across everything (it works for groups as well)

penn_sf <- penn_sf %>%
  mutate(pred_cases = fitted(gam_mrf)*population,
         pred_rate = fitted(gam_mrf)*1000)

penn_predict_agg <- penn_sf %>%
  group_by(county) %>%
  summarise(population = sum(population), 
            pred_cases = sum(pred_cases),
            pred_rate = sum(pred_cases)/population*1000,
            observed_rate = sum(cases)/population*1000) 

#plot predictions
ggplot(penn_predict_agg) +
  geom_sf(aes(fill = pred_cases)) +
  scale_fill_viridis_c(option = "D", guide = guide_colorbar("Cancer Cases")) +
  theme_minimal() + theme(legend.position="bottom")


Which looks rather like the observations!

ggplot(penn) +
  geom_sf(aes(fill = cases)) +
  scale_fill_viridis_c(option = "D", guide = guide_colorbar("Cancer Cases")) +
  theme_minimal() + theme(legend.position="bottom")

enter image description here

  • I'm going to argue that the model is not necessarily a binomial one. The implication for the binomial is that you conducted n Bernoulli trials per sub-population (combination of age, race, gender) & had so many successes out the n trials. But do you think there were actually n trials? Instead you can think of the data as being the count of a number of events in each region plus some covariates that might explain differences in expected counts. It is reasonable to assume that higher counts will be associated with larger populations, so we normalize to the (sub)population size using an offset. – Gavin Simpson May 6 at 22:34
  • Both models are reasonable, but the Poisson model, to me at least, seems most reasonable given that we didn't subject the entire populations of these regions to the trial. At least that's how I see it. – Gavin Simpson May 6 at 22:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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