How to simulate quantities of interest using arm or rstanarm packages?

[Danilo Freire]

Hello everyone,

Very newbie question here. I would like to know how to simulate quantities of interest out of a regression model estimated using either the arm or the rstanarm packages in R. I am just starting to use Bayesian methods and R and have been using the Zelig package for some time. I asked a similar question before, but I would like to know if it is possible to simulate those quantities using the posterior distribution estimated by those packages.

In Zelig you can set the values you want for the independent values and it calculates the results for the outcome variable (expected value, probability, etc). An example:

# Creating a dataset:
x <- rnorm(100,20,10)
z <- rnorm(100,10,5)
e <- rnorm(100,0,1)
y <- 2*x+3*z+e
df <- data.frame(x,z,e,y)

# Loading Zelig
require(Zelig)

# Model
m1.zelig <- zelig(y ~ x + z, model="ls", data=df)
summary(m1.zelig)

# Simulating z = 10
s1 <- setx(m1.zelig, z = 10)
simulation <- sim(m1.zelig, x = s1)
summary(simulation)

So Zelig keeps x at its mean (20.56), and simulates the quantity of interest with z = 10. In this case, y is approximately 71.

The same model using arm:

# Model
require(arm)
m1.arm <- bayesglm(y ~ x + z, data=df)
summary(m1.arm)

And using rstanarm:

# Model
require(rstanarm)
m1.stan <- stanlm(y ~ x + z, data=df)
print(m1.stan)

Is there any way to simulate z = 10 and x equals to its mean with the posterior distribution estimated by those two packages and get the expected value of y? 

Thank you very much for your attention,

Danilo

[Bob Carpenter]

rstanarm is an experimental interface that I don't think
is finished.

You can do the simulation directly in Stan using
the generated quantities block and Stan's _rng functions
for distributions.

But you'll need to provide the model error and coeffs
as data.

- Bob

[Ben Goodrich]

Like Bob said, you can do things like this in the generated quantities block of a .stan file (the ones in rstanarm are hardcoded without generated quantities in addition to being unfinished).

To do it in R, you need to do setx() yourself, but the posterior simulations of the parameters can be obtained in various ways. Probably the route that is most similar to Zelig would be

stanfit <- stan(...)
posterior <- as.matrix(stanfit)

You could then matrix multiply the x vector by (some subset of the columns of) posterior and then go on to construct your quantities of interest.

For arm, there is also a sim() function which draws parameters from the estimated sampling distribution.

Ben

<pre class="lang-r prettyprint prettyprinted" style="margin-bottom: 10px; padding: 5px; overflow: auto; width: auto; max-height: 600px; word-wrap: normal; line-height: 17.804800033569336px; background-image: initial; background-attachment: initial; background-size: initial; background-origin: initial; b

...

[Danilo Freire]

Dear Bob and Ben,

Thank you very much for your replies. Could I ask you for a working example using the code I wrote above setting x and z at different values? It would probably clarify some of my doubts. 

Best,

Danilo