1

I am wondering if one can apply colors directly to parameters from rstanarm models with bayesplot

For example from https://mc-stan.org/bayesplot/articles/plotting-mcmc-draws.html:

library("bayesplot")
library("ggplot2")
library("rstanarm")

fit <- stan_glm(mpg ~ ., data = mtcars, seed = 1111)
posterior <- as.array(fit)

color_scheme_set("red")
mcmc_intervals(posterior, pars = c("cyl", "drat", "am", "sigma"))

returns:

enter image description here

The documentation describes custom schemes, or mixing, but they still apply to every parameter in the plot:

color_scheme_set("mix-blue-red")

plot(fit, pars = c("cyl", "drat", "am", "sigma"))

enter image description here

The plot function (with rstanarm model) no longer accepts a col argument to be able to specify each point. Is there anyway to specify a string of colors (or schemes) for each parameter in the plot?

Edit: Based on comments from @Limey, you can modify plots as ggobjects, but I don't see a straightforward way to access the original plot aesthetics.

For example, scale_colour_manual (or any other scale_colour_... for that matter) does not alter the plotted values:

plot(fit, pars = c("cyl", "drat", "am", "sigma"))+
  scale_colour_manual(values=c("blue", "purple", "orange","red"))

(see plot above, it looks the same). Another example, I can plot points, on top of the original:

plot(fit, pars = c(names(mtcars)[-1]))+
  geom_point(aes(x=fit$coefficients[-1],y=names(fit$coefficients)[-1]),color="black")

enter image description here

and make them bigger to match (note I have also made them pink here):

plot(fit, pars = c(names(mtcars)[-1]))+
  geom_point(aes(x=fit$coefficients[-1],y=names(fit$coefficients)[-1]),color="pink",size=3)

enter image description here

I could do this for each of the errorbar objects as well, but at this point I am not really modifying the original plot, but rather plotting right over the top of it. In this case it makes sense to ditch the rstanarm plot function all together and simply do it manually, or create my own function, which is fine. However, if anyone knows how to do modify the plot(rstanarm_object) colors directly, I would still be interested to hear it.

Dylan_Gomes
  • 2,066
  • 14
  • 29
  • I'm not familar with `rstanarm`, but the [online doc[(https://mc-stan.org/rstanarm/reference/plot.stanreg.html) tells me that `plot()` returns "... a ggplot object that can be further customized using the ggplot2 package", so you can manipulate it with standard `ggplot2` functions once `rstanarm` has finished with it. – Limey Jun 16 '20 at 06:55
  • Yes, but I am not familiar with any way to access aesthetics from the `rstanarm` object. Adding geom_point(), for example, tries to add new points, and needs new aesthetics, so it seems like it doesn't allow you to change the rstanarm object. – Dylan_Gomes Jun 16 '20 at 14:00
  • `geom_point()` would indeed add new points. But `scale_colour_discrete()`, for example, allows you to customise the appearance of colour aesthetics and the legends associated with them. There are similar functions for every aesthetic, and for both discrete and continuous scales. `scale_XXXX_manual()` is the most complex but also the most flexible. There are other functions that allow you to customise axes and you can repeat `geom_XXXX` calls to modify choices `rstanarm` has made, just as you have done above.. Please take a look at the vignettes and other info on the tidyverse website. – Limey Jun 16 '20 at 14:56
  • Thanks for the comment, and I have used all of these functions before. But there is something different about the plotted objects that come from `rstanarm`, as none of the `scale_` functions change the plotted values. Although I can still scale axes. I have added that in my question above. – Dylan_Gomes Jun 16 '20 at 16:18

1 Answers1

1

This doesn't solve the issue directly, but it may still be a workaround that could be improved upon, if no other solutions arise.

Stemming from the edits to the OP via comments from @Limey:

You can just plot everything yourself from the model object, which gets close to the plotted values. I realize that this isn't perfect, and perhaps someone has a suggestion as to how these values can perfectly overlap. rstanarm plots default to 50% and 90% credible intervals. I use t tables here and 50% and 90% confidence intervals, which aren't the same, but it comes out close.

color_scheme_set("gray")

vars<-c(2,5,9,11) # select values from fit$coefficients 
                  # (this is assumed to match with names(mtcars))
                  # Note I have switched `sigma` to `carb` for ease of model retrieval

plot(fit, pars = c(names(mtcars)[vars]))+ # the original plot, now in grey

# the first error bar. With a t table, 1.645 is a 90% confidence interval,
# which isn't the same here as the 90% credible interval from rstanarm, but it 
# gets us close
  geom_errorbarh(aes(xmin=(fit$coefficients[vars]-fit$ses[vars]*1.645),
                     xmax=(fit$coefficients[vars]+fit$ses[vars]*1.645),
                     y=names(fit$coefficients)[vars]),
                 height=0,alpha=0.5,color="yellow" # transparency is so you can see grey through the yellow.. I know this isn't pretty
                     )+

  geom_errorbarh(aes(xmin=(fit$coefficients[vars]-fit$ses[vars]*.675), # this is .675, which corresponds to a 50% confidence interval (again not the same as a credible interval, but it is getting us close
                     xmax=(fit$coefficients[vars]+fit$ses[vars]*.675),
                     y=names(fit$coefficients)[vars]),
                 height=0,lwd=2,alpha=0.5,color="yellow"
                     )+
# now plot the points over the lines:
  geom_point(aes(x=fit$coefficients[vars],
y=names(fit$coefficients)[vars]),
color="yellow",size=3,alpha=0.5)

enter image description here

This isn't pretty, and it is not exact (you can see a bit of grey to the edges of the yellow), but it is close. If you can modify each value set color, then you can alternate colors:

colors=c("red","blue","red","blue") # selecting colors one by one
colors=rep(c("red","blue"),2) # or just have alternating colors repeat (2 times here)

ggplot()+ # note that I have removed `rstanarm` object, which is no longer necessary
  geom_errorbarh(aes(xmin=(fit$coefficients[vars]-fit$ses[vars]*1.645),
                     xmax=(fit$coefficients[vars]+fit$ses[vars]*1.645),
                     y=names(fit$coefficients)[vars]),
                 height=0,alpha=0.5,color=colors
                     )+
  geom_errorbarh(aes(xmin=(fit$coefficients[vars]-fit$ses[vars]*.675),
                     xmax=(fit$coefficients[vars]+fit$ses[vars]*.675),
                     y=names(fit$coefficients)[vars]),
                 height=0,lwd=2,alpha=0.5,color=colors
                     )+
  geom_point(aes(x=fit$coefficients[vars],
y=names(fit$coefficients)[vars]),
color=colors,size=3,alpha=0.5)

enter image description here

or something entirely wacky:

colors=c("red","blue","green","yellow")
colors2=rev(colors)

ggplot()+
  geom_errorbarh(aes(xmin=(fit$coefficients[vars]-fit$ses[vars]*1.645),
                     xmax=(fit$coefficients[vars]+fit$ses[vars]*1.645),
                     y=names(fit$coefficients)[vars]),
                 height=0,alpha=0.5,color=colors
                     )+
  geom_errorbarh(aes(xmin=(fit$coefficients[vars]-fit$ses[vars]*.675),
                     xmax=(fit$coefficients[vars]+fit$ses[vars]*.675),
                     y=names(fit$coefficients)[vars]),
                 height=0,lwd=2,alpha=0.5,color=colors
                     )+
  geom_point(aes(x=fit$coefficients[vars],
                 y=names(fit$coefficients)[vars]),
             color=colors2,size=3,alpha=0.5)

enter image description here

Dylan_Gomes
  • 2,066
  • 14
  • 29
  • While there may be a way to do this within the bayesplot package, I wonder if using the the plot object returned by `p = plot(fit)` is slightly easier i.e. `ggplot(p$data, aes(x=m, y=parameter, xmin=ll, xmax=hh, col=parameter)) + geom_point(size=1.5) + geom_errorbarh(height=0, size=0.5) + geom_errorbarh(aes(xmin=l, xmax=h), height=0, size=1)`. Then, similar to your answer, all the usual ggplot tools are available – user20650 Oct 07 '20 at 15:51