Y axis scale of plot_nonlinear

Interfaces
,
1

I’ve fit a stan_gamm4 model using log transformed data (that is otherwise all positive) and gaussian family (link = ‘identity’) with one smooth and when I plot the smooth using plot_nonlinear, the scale on the y-axis is not on the real world scale. So my first questions are: What is the scale on the y-axis, and why is it not on the same scale as the input (why would you not want the model fit on the same scale you’re working on)?

Further, playing around a little with the plot_nonlinear function I’ve extracted the snippet of code pasted below, with a few minor modifications, which looks like it brings back the data on the scale of my input prior to log transformation. It’s unclear to me whether this is coincidental or whether there’s a logical explanation to it, perhaps someone is able to clarify this?

plotdata.func <- function(x){
  smooths <- 's(the.smooth)'
  XZ <- x$x
  XZ <- XZ[, !grepl("_NEW_", colnames(XZ), fixed = TRUE)]
  labels <- sapply(x$jam$smooth, "[[", "label")
  xnames <- sapply(x$jam$smooth, "[[", "vn")
  names(x$jam$smooth) <- labels
  names(xnames) <- labels
  fs <- sapply(x$jam$smooth, FUN = "inherits", what = "fs.interaction")
  B <- as.matrix(x)[, colnames(XZ), drop = FALSE]
  original <- x$jam$model
  
  df_list <- lapply(x$jam$smooth[smooths], FUN = function(s) {
    incl <- s$first.para:s$last.para
    b <- B[, incl, drop = FALSE]
    xz <- XZ[, incl, drop = FALSE]
    x <- original[, s$term]
    ord <- order(x)
    x <- x[ord]
    xz <- xz[ord, , drop = FALSE]
    f <- rstanarm:::linear_predictor.matrix(b, xz)
    data.frame(predictor = x, 
               lower95 = apply(f, 2, quantile, probs = (1 - 0.95)/2), upper95 = apply(f, 2, quantile, probs = 0.95 + (1 - 0.95)/2), 
               lower90 = apply(f, 2, quantile, probs = (1 - 0.90)/2), upper90 = apply(f, 2, quantile, probs = 0.90 + (1 - 0.90)/2), 
               lower80 = apply(f, 2, quantile, probs = (1 - 0.80)/2), upper80 = apply(f, 2, quantile, probs = 0.80 + (1 - 0.80)/2), 
               middle = apply(f, 2, median), term = s$label)
  })
  plot_data <- do.call(rbind, df_list) %>% distinct()
  plot_data[, which(colnames(plot_data) == 'lower95'):which(colnames(plot_data) == 'middle')] <- plot_data[, which(colnames(plot_data) == 'lower95'):which(colnames(plot_data) == 'middle')] + mean(tmp$Cstock_kgm2)
  plot_data <- plot_data %>% distinct()
  rownames(plot_data) <- NULL
  return(plot_data)
}

Operating System: Windows 10

2

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

image
image1443×699 75.3 KB

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")

image
image1443×699 46.2 KB

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.

1 Like
3

This did indeed help a good bit! Thanks!

Can I ask what the re_formula = ~0 does? The documentation only mentions NULL and NA as far as I can see.

I do still find it strange that it would come out on the real scale, as if it had never been log-transformed, when I run my own function. I should add that you would need to change the bit where it says ‘s(the.smooth)’ to whatever your smooth is called, and that the x that you’re sending to the function should be the model object, in your case br.

4

It’s equivalent to re_formula = NA, meaning that the predictions are not conditioned on any of the group-level “random” effects in the model. See ?rstanarm::posterior_predict.stanreg for this documentation.