Obtaining effect size from “rstanarm” package's linear regression

I have read this thread and have an additional question.

I have a multilevel factor in a multivariable stan_lm model.

A journal editor requests that the total contribution to outcome variance be estimated for this factor.

I have estimated the POV with the code shown above for each level.

PPD_v <- posterior_predict(fit, seed = same value for all POV)
nd <- model.frame(fit)
nd$x1 <- ‘factor level’
PPD_c <- posterior_predict(fit, newdata = nd, seed = same value for all POV)
POV <- 1 - apply(PPD_c, 1, var) / apply(PPD_v, 1, var)
quantile(POV, c(.025, .975))

I can sum the individual POVs and then apply median and quantile functions.

The values appear reasonable, but have I violated relevant assumptions?


I would say that is the right code, but it is a dubious request by the journal editor and the endpoints of a 95% interval are not estimated particularly accurately with the default number of draws.

My posterior sample size is quite large.

My PPD objects have draws set at 1000.

The values for the POV are approximately median 0.27 (0.26, 0.28).

I would have thought that that number of draws would be sufficient for accurate estimation of the endpoints.

Also, the Bayesian R2 for the model is 0.42.

I have a more philosophical question. I have been convinced from the expositions by you and Gelman and others that I should move from null hypothesis testing to parameter estimation.

In the medical journals (anesthesia) where I publish, Bayesian analyses are very rare.

I will probably have to add the factor POV that has been requested. But what are the main points with which I should argue against this request?


I wouldn’t bother arguing against it, but I have a hard time thinking of a reason why it would matter if the proportion were 0.27 vs 0.1 or whatever.