erc
erc

Reputation: 10133

How to access plot data of GAM (tensor product smooth) in R?

I have such a GAM model with a tensor product smooth of location and time:

dat <- gratia::bird_move
dat$longitude <- dat$latitude

mod <- mgcv::gam(count ~ species + te(longitude, latitude, week, d = c(2,1), k = c(8,8)), data = dat)

enter image description here

I know it's possible to get the data from plot(mod) but how can I extract (and plot) data for e. g. week 48 only?

plotdata <- plot(mod) 
> str(plotdata)
List of 1
 $ :List of 12
  ..$ scale  : logi FALSE
  ..$ se     : logi FALSE
  ..$ m      : num 40
  ..$ nc     : num 3
  ..$ nr     : num 3
  ..$ lo     : num [1:3] 5 5 0
  ..$ hi     : num [1:3] 60 60 52
  ..$ vname  : chr [1:3] "longitude" "latitude" "week"
  ..$ main   : chr "te(longitude,latitude,week,37.46)"
  ..$ exclude: logi [1:14400] FALSE FALSE FALSE FALSE TRUE TRUE ...
  ..$ fit    : num [1:14400, 1] 4.98 4.9 4.83 4.74 NA ...
  ..$ plot.me: logi TRUE

Upvotes: 2

Views: 462

Answers (1)

Allan Cameron
Allan Cameron

Reputation: 173803

Typically, one one create a data frame of the input variables then run predict using these data. For example, to get a complete grid of predictions for counts of species 1 at each latitude / longitude in week 48, we could do:

pred_df <- data.frame(species = 'sp1', 
                      longitude = rep(0:60, 61), 
                      latitude = rep(0:60, each = 61), 
                      week = 48)

pred_df$count <- predict(mod, newdata = pred_df, 
                         terms = "te(longitude,latitude,week)")

To print this we could use ggplot and geomtextpath:

library(geomtextpath)

ggplot(pred_df, aes(longitude, latitude)) +
  geom_raster(aes(fill = count), interpolate = TRUE) +
  geom_textcontour(aes(z = count, 
                       label = ifelse(after_stat(level) < 0,
                                      paste0('\u2010', abs(after_stat(level))),
                                      after_stat(level))),
                   family = 'serif') +
  scale_fill_gradientn(colours = c('blue4', 'orange', 'yellow', 'white')) +
  coord_cartesian(expand = FALSE) +
  theme_minimal(base_size = 16) +
  ggtitle('Relative effect of latitude\nand longitude on counts at week 48')

enter image description here

Upvotes: 1

Related Questions