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 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.