Welcome, @questionmark! I’m sorry it’s taken a bit for your question to get answered.
plot_nonlinear shows you the marginal smooths. A way to interpret these is “the additional effect of {smooth predictor} on the response variable.” I think these will always include both positive and negative values as a result—they are deviations from the value predicted by the other terms in the model.
It sounds like you may be looking for a plot of the expected or predicted posterior distributions as a function of your smooth predictor instead. I’m going to use some simulated data to illustrate this idea.
library(rstanarm)
library(modelr)
library(tidybayes)
# simulate data
set.seed(946)
dat <- mgcv::gamSim(1, n = 400, scale = 2)
dat$fac <- fac <- as.factor(sample(1:20, 400, replace = TRUE))
dat$y <- dat$y + model.matrix(~ fac - 1) %*% rnorm(20) * .5
# fit model
br <- stan_gamm4(y ~ s(x0) + x1 + s(x2), data = dat, random = ~ (1 | fac),
chains = 2, cores = 2)
This model has population-level (“fixed”) effects for the intercept and x1, smooths for x0 and x2, and a group-level (“random”) intercept for fac. Looking at the population-level terms with fixef(br), we see
(Intercept) x1 s(x0).1 s(x0).2 s(x0).3
4.78493241 6.47683065 0.09113933 0.22991484 0.59671445
s(x0).4 s(x0).5 s(x0).6 s(x0).7 s(x0).8
0.24298612 -0.11503229 -2.85732828 0.05528272 -1.19817080
s(x0).9 s(x2).1 s(x2).2 s(x2).3 s(x2).4
0.03140069 -52.01086897 30.73428316 6.38281162 -12.58550181
s(x2).5 s(x2).6 s(x2).7 s(x2).8 s(x2).9
28.40463723 14.48973771 -15.46259395 -13.32189433 3.52970797
And then looking at the marginal smooths with plot_nonlinear(br), we see
Based on these, we should be able to estimate the expected value of y when all of the predictors are equal to zero. This will be (Intercept) plus the marginal smooth values at 0 for x0 and x2. Eyeballing the output table and plot, I’d guess y = 4.8 - 1.2 - 3.5 = 0. If we continue to hold all the predictors at zero, but increase the value of x2 = 0.5we should expect y = 4.8 - 1.2 - 0.7 = 2.9. Finally, with x2 = 1, we wexpect y = 4.8 -1.2 - 4.5 = -0.9.
We can do these calculations properly for rstanarm models using tools from the tidybayes package (link goes to a vignette for making plots from rstanarm fits).
postx2 <- dat |>
data_grid(x0 = 0, # generate prediction dataset
x1 = 0,
x2 = seq_range(x2, n = 100)) |>
add_epred_draws(br, re_formula = ~0) # add expected value draws from the posterior
ggplot(data = postx2, aes(x = x2, y = .epred))+
stat_lineribbon()+ # visually summarize the posterior of expected values
annotate("point", x = 0, y = 0, size = 4, color = "red")+ # eye-ball estimates
annotate("point", x = 0.5, y = 2.9, size = 4, color = "blue")+
annotate("point", x = 1, y = -0.9, size = 4, color = "yellow")+
scale_fill_brewer(palette = "Greys")
Hot dog! We got the expected values for a range of x2 with all other predictors at zero and our three validation points are pretty close for a rough eyeball estimate. You can supply any combination of predictor values you want. My goal in choosing zeros for the other predictors was to make the calculations easier to follow. Marginal plots like these more typically hold the other predictors at their means instead of at zero.
I wasn’t able to run your code without example data, but I admire your determination to dig in there. In general, I recommend using the tools that have been built by the wonderful Stan community to post-process and understand your fits. This way, you’ll avoid both the effort of writing a different custom function to handle every model you fit and avoid the risk of coincidentally passing visual inspection for the wrong reasons.
I hope this helps.