# Errors fitting a Markvov Random Field GAM with multiple samples per Polygon

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.

``````library(mgcv)
library(sf)
library(dplyr)
library(SpatialEpi)
library(spdep)

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.

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

#neighborhoods

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

Let's see it!

``````plot(penn_sf[["geometry"]])
plot(nb, coords =coords, col = "red", add = TRUE)
`````` 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
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.

Very useful post.

A GAM with Poisson likelihood and identity offset does indeed work if you account for the county (#180) with zero population - the obscure error alert from `mgcv` is referring to trying to take the log of zero. So try this (although other likelihoods might be more appropriate here):

```gam_mrf <- gam(cases ~ s(county, bs = 'mrf', xt = list(nb = nb)) + race + gender + age + offset(log(population+0.01)), data = penn_sf, method = "REML", family = poisson)```

Or nonproportional risk as:

```gam_mrf1 <- gam(cases ~ s(county, bs = 'mrf', xt = list(nb = nb)) + race + gender + age + s(population), data = penn_sf, method = "REML", family = poisson)``` Either works.

Also the code was somewhat incomplete - for instance the following is missing:

`coords <- st_centroid(st_geometry(penn))`

and I think you meant this:

```penn_sf <- penn_sf %>% mutate(pred_cases = fitted(gam_mrf))```

... rather than this:

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

The results for the GAM with Poisson likelihood `gam_mrf` using this code are shown below in addition to your code where appropriate (such as `penn` and `nb`):

``` require(tidyverse) require(colorspace)```

```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)```

```p.pred<- ggplot(penn_predict_agg) + geom_sf(aes(fill = pred_cases)) + scale_fill_continuous_sequential(palette="Viridis",rev=FALSE, breaks=c(100,500,900,1300)) + guides(fill=guide_colorbar(title="cancer\n cases",ticks=FALSE, nbin=30,barwidth=0.8,barheight=7)) + theme_minimal() + theme(legend.position="right")```

```p.obs<- ggplot(penn) + geom_sf(aes(fill = cases)) + scale_fill_continuous_sequential(palette="Viridis",rev=FALSE, breaks=c(100,500,900,1300)) + guides(fill=guide_colorbar(title="cancer\n cases",ticks=FALSE, nbin=30,barwidth=0.8,barheight=7)) + theme_minimal() + theme(legend.position="right")```

```cowplot::plot_grid(p.obs,p.pred,labels=c("observed","predicted"), ncol=1,align="v",hjust=c(-3,-3.1),vjust=2)``` Milani

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)

#aggregate
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")
`````` • 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. – Reinstate Monica - G. 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. – Reinstate Monica - G. Simpson May 6 at 22:35