Understanding Prior in rstanarm

Question

I want to to fit y~intercept+x.1*x+x.2*x^2.

Here is the indata. The items "ab" and "ad" have the same underlying formula as "aa" and "ac".

library(data.table)
library(ggplot2)
set.seed(123)
DT.1 <- data.table(item= "aa", x=1:100)
DT.1[, y:=0+x*40-x^2*0.4+rnorm(100, mean=0, sd=100)]
DT.2 <- data.table(item= "ab", x=1:10)
DT.2[, y:=0+x*40-x^2*0.4+rnorm(10, mean=0, sd=100)]
DT.3 <- data.table(item= "ac", x=1:100)
DT.3[, y:=0+x*50-x^2*0.4+rnorm(100, mean=0, sd=100)]
DT.4 <- data.table(item= "ad", x=61:70)
DT.4[, y:=0+x*40-x^2*0.4+rnorm(10, mean=0, sd=100)]
DT <- rbind(DT.1, DT.2, DT.3, DT.4)
p <- ggplot(DT, aes(x=x, y=y, group=item))
p <- p + geom_point()
p <- p + facet_wrap(~item)
p

But when doing a normal lm, the fitting of "ab" and "ad" is not what I need.

DT[, c("intercept", "x.1", "x.2"):=as.list(coef(lm("y~x+I(x^2)", .SD))), by = list(item)]
DT[, x.index:=1:.N, by=item]
DT.coef <- DT[x.index==1, list(item, intercept, x.1, x.2)]
DT.1 <- data.table(item="aa", x=1:100)
DT.1[, y.lm:=DT.coef[item=="aa", intercept]+x*DT.coef[item=="aa", x.1]+x^2*DT.coef[item=="aa", x.2]]
DT.2 <- data.table(item="ab", x=1:100)
DT.2[, y.lm:=DT.coef[item=="ab", intercept]+x*DT.coef[item=="ab", x.1]+x^2*DT.coef[item=="ab", x.2]]
DT.3 <- data.table(item="ac", x=1:100)
DT.3[, y.lm:=DT.coef[item=="ac", intercept]+x*DT.coef[item=="ac", x.1]+x^2*DT.coef[item=="ac", x.2]]
DT.4 <- data.table(item="ad", x=1:100)
DT.4[, y.lm:=DT.coef[item=="ad", intercept]+x*DT.coef[item=="ad", x.1]+x^2*DT.coef[item=="ad", x.2]]
DT.lm <- rbind(DT.1, DT.2, DT.3, DT.4)
p <- ggplot(DT.lm, aes(x=x, y=y.lm, group=item))
p <- p + geom_point()
p <- p + facet_wrap(~item)
p

So I would like to set a prior that makes the "ab" and "ad" be much more similar to the "aa" and "ac".

Trying to replicate examples from rstanarm documentation and vignettes, I get stuck in how to define the prior. This is the closest I get, and it produces the error "Error in prior$location : $ operator is invalid for atomic vectors".

library(rstanarm)
y.posterior <- stan_lm(y~x+I(x^2), data=DT[item=="aa", ], prior=mean(0.5), chains=1, cores=1, iter=1000, seed=12345)

What (conceptually) is the prior of y~intercept+x.1*x+x.2*x^2? Should I define mean and sd for intercept, x.1, and x.2? I can "draw" how the prior curve should look (peak around x=50 with y around 1100), which translates to intercept=18, x.1=40, and x.2=-0.4.

How do I translate that to the prior rstanarm wants?

I hope my trying to use advanced statistical methods with so little knowledge of statistical science will not insult anyone here. Please do not comment on the input data, it is only there to make the code run.

Answer 1

First stan_lm only accepts one type of prior:

Must be a call to R2 with its location argument specified or NULL, which would indicate a standard uniform prior for the R^2.

If you want to specify any of the other types of priors, you'll need to switch to stan_glm. With the default family (gaussian), it will be equivalent to stan_lm.

Second, you've tried to specify a prior of mean? My first thought is to make sure you understand what a prior is. It's a probability distribution of values for (in this case) your parameters. You have three parameters: Intercept, $x_1$ and $x_2$, so you will either need to specify three priors or specify one prior that applies to all three.

I can't tell, but I think you're thinking of your function's curve as something that's driven directly by priors, but that's not the case. Priors influence parameters, parameters drive your function.

When I run your code, I get:

item intercept x.1 x.2 aa 18.00042 38.978120 -0.3873962

so perhaps you're looking for a prior for your intercept of something like normal (20, 5), and for $x_1$ of normal (40, 10) and for $x_2$ of normal (0, 4). Or something like that. You can see the implications with something like:

curve (dnorm (x, 20, 5), 0, 40)

which is the probability distribution that you'd get for intercept if you used normal (20, 5) in stan_glm. This prior is saying that the intercept is most likely to be around 20, and very unlikely to be as small as 5 or as large as 35. (Based on eyeballing the curve graph.) stan_glm will combine this prior probability with the actual data (that you generated, per your original curve) and will tell you what values of intercept are likely based on the prior and the data. That's your posterior.

Your curiosity is good, but you won't be able to master a tool as powerful as rstanarm without a more solid understanding of the principles on which it operates -- Bayesian inference. You're on the right track, since Bayesian Inference is extremely powerful and flexible, just work more on your foundations. If you want to get an excellent book, I can recommend Statistical Rethinking: A Bayesian Course with Examples in R and Stan by Richard McElreath. It's a bit pricey, and it will take serious study to get to the end of it, but it starts at the basics, is well-written, and pretty much tells you everything you need to know about Bayesian Inference.