not all stanreg point estimates are posterior medians

Summary:

Not all point estimates in stanreg objects are actually posterior median values (notably fitted.values, residuals, and, for mixed effects regressions, estimates returned by coef).

Description:

In the rstanarm:::stanreg function, the median values of the regression coefficients are taken and then used to calculate the fitted.values element of the stanreg object. However, since the median (unlike the mean) is not a linear function, the resulting fitted values are not the medians of the posterior distribution for the fitted values. This then effects the residuals element of the object, and a similar issue exists for the calculation of the estimates returned by coef for mixed effects models, since the median fixed effects values and median random effects values are being added together in rstanarm:::coef.stanreg.

Reproducible Steps:

I will demonstrate with the fitted values for linear mixed effects regressions, because this is a simple case; but the issue generalizes at least to the other cases mentioned above, and for other regression families as well. The following two functions should return TRUE and a message indicating the Mean relative difference, respectively, for any model fit with stan_lmer:

rstanarm_fitted_values <- function(object) {
  if (!inherits(object, "stanreg") || !inherits(object, "lmerMod") ||
  !isTRUE(all.equal(object$family, gaussian()))) {
    stop("Function only tested on linear mixed effects regressions.")
  }
  
  modmat <- posterior_linpred(object, XZ = TRUE)
  coefmat <- as.matrix(object)[, 1:ncol(modmat), drop = FALSE]
  coefs <- t(apply(coefmat, 2, median))
  point_estimates <- as.vector(coefs %*% t(modmat))
  
  return(all.equal(unname(point_estimates), unname(object$fitted.values)))
}

median_fitted_values <- function(object) {
  if (!inherits(object, "stanreg") || !inherits(object, "lmerMod") ||
  !isTRUE(all.equal(object$family, gaussian()))) {
    stop("Function only tested on linear mixed effects regressions.")
  }
  
  modmat <- posterior_linpred(object, XZ = TRUE)
  coefmat <- as.matrix(object)[, 1:ncol(modmat), drop = FALSE]
  mu <- coefmat %*% t(modmat)
  posterior_medians <- apply(mu, 2, median)
  
  return(all.equal(unname(posterior_medians), unname(object$fitted.values)))
}

RStanARM Version:

2.15.3

R Version:

3.3.3

Operating System:

Windows 10

Start with rstanarm:::stanreg and coef.stanreg, then run the supplied stan_lmer reproducible example to compare fitted.values with posterior medians. Trace the related residuals and coef calculations as well. Done means the affected point estimates represent the intended posterior medians and the regression-model behavior is covered by tests.