Recovering parameters from stanreg object

Interfaces rstanarm
1

Hi everyone! I am a novice with bayesian and am I simulated data, then fit it using the stan_lmer function:


set.seed(123)
a <- 1.5
b <- 0.4
sigma_y <- 0.3

n_perspp <- 25
n_spp <- 4
n_ids <- n_perspp * n_spp
rep <- 3
N <- n_ids * rep

# set spp names and ids 
spp <- c("alninc", "betall", "betpap", "betpop")
tree_ids <- paste0(rep(spp, each = n_perspp), 1:n_perspp)  # 100 unique tree IDs
ids <- rep(tree_ids, each = rep)  # repeat each tree ID 3 times
tree_spp <- rep(rep(spp, each = n_perspp), each = rep)  # matching species for each ID

# partial pooling
sigma_ids <- 0.8 / 2.57
sigma_spp <- 1 / 2.57
a_ids_values <- rnorm(n_ids, 0, sigma_ids)
a_spp_values <- rnorm(n_spp, 0, sigma_spp)

# match to observations
a_ids <- rep(a_ids_values, each = rep)  # repeat each a_ids 3 times
a_spp <- rep(a_spp_values, each = n_perspp * rep)  # repeat each a_spp 75 times

# gdd and devide by constant
gdd <- round(rnorm(N, 1800, 100))
gddcons <- gdd / 200

# error
error <- rnorm(N, 0, sigma_y)

# calculate ring width
ringwidth <- a + a_ids + a_spp + b * gddcons + error

# set df
simcoef <- data.frame(
  ids = ids,
  spp = tree_spp,
  gddcons = gddcons,
  b = b,
  a = a,
  a_ids = a_ids,
  a_spp = a_spp,
  sigma_y = sigma_y,
  ringwidth = ringwidth
)
# run models
  fit <- stan_lmer(
    ringwidth ~ gddcons + (1 | ids) + (1 | spp),  
    data = simcoef,
    chains = 4,
    iter = 4000,
    core=4
  )

Now I want to recover my parameter estimates from my model, but am having a lot of trouble finding out the right method to do it! I used bot coef() and as.data.frame(summary(fit)) functions, but they return slightly different things. By reading the stan manual, I realized that the coef() intercept estimates were the sum of fixed and random parameters. But how do I recover each intercept estimate if I have more than 1 hierarchical group (here ids and spp)?

I couldn’t recover my standard errors with coef(), so I used the summary function, but now my model overestimates my spp intercept.

Could anyone help me with that? I tried to me as clear as possible, but let me know if I should add any details/code or something else.

Thanks!

2

Hi, welcome to the Stan forums. There are quite a few ways to get information out of these objects. Here are some useful ones:

  • as.matrix(fit) will give you all the posterior draws for all of the parameters. Each parameter will be a column in the matrix and each row represents a draw from the posterior. Using this matrix you can compute any summaries of them you want (including everything I mention below).
  • summary(fit) gives you point estimates of all model parameters including mean and median. One clarification: the mcse column gives you the Monte Carlo Standard Error (this would theoretically go to 0 if you could run infinitely long MCMC). The sd column gives you the estimated standard deviation of the posterior distribution of the parameter.
  • ranef(fit) will give you point estimates of the so-called random effects. The output has the same structure as the coef(fit) output but it doesn’t do the summation you mentioned that coef does.
  • posterior_interval(fit) gives you posterior uncertainty intervals. The prob argument can control the amount of probability mass included in the interval (the default is 90% intervals I think).

I think the most important one is as.matrix(fit) since everything else can be computed from the posterior draws.

1 Like