I would like to extract the fitted values from my bam() model in order to be able to plot the figures in either prism or with ggplot. I use GAM to model brain responses in EEG data, the data are thus many 400ms time series of four different categories. The model contains a by-factor for these categories(uV ~ s(time, by = CatInt)). I would like to get the fitted values for these 4 categories separately.
I found the function fitted() or fitted.values() but this just outputs a string with the same length as the number of data points I put into the model, while I would like one average for the data over time per category. Ideally it would include either SD, SE or CI.
Is this possible?
If you are interested in customization of the plot.gam() output, you may keep the plot as an object.
First of all, let's generate some testing data with the built-in gamSim( ) function and construct a model:
library(mgcv)
set.seed(2) # R version 3.6.0 used
dat <- gamSim(5, n = 200, scale = 2)
b <- bam(y ~ x0 + s(x1, by = x0), data = dat)
Then you may just keep the plot object which is a list with all the data needed to generate a customized plot:
pl_data <- plot(b, pages = 1)
> str(pl_data)
List of 4
$ :List of 11
..$ x : num [1:100] 0.00711 0.01703 0.02694 0.03686 0.04678 ...
..$ scale : logi TRUE
..$ se : num [1:100] 3.25 3.04 2.83 2.64 2.47 ...
..$ raw : num [1:200] 0.185 0.702 0.573 0.168 0.944 ...
..$ xlab : chr "x1"
..$ ylab : chr "s(x1,3.33):x01"
..$ main : NULL
..$ se.mult: num 2
..$ xlim : num [1:2] 0.00711 0.9889
..$ fit : num [1:100, 1] -4.76 -4.64 -4.52 -4.4 -4.28 ...
..$ plot.me: logi TRUE
....
Finely, you may extract the pl_data components which are of interest with usual subsetting approaches, for example:
first_item_pl <- pl_data[[1]][c("x", "fit", "se")]
to extract different items of the first list or to extract the same item from all four nested lists
fit_list <- lapply(pl_data, "[[", "fit")