Reputation: 10133
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)
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
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')
Upvotes: 1