It appears that every invocation of `posterior_linpred`

and `posterior_epred`

gives a different permutation of the draws. This prevents users from doing posterior calculations (e.g. covariances) using the draws.

Relatedly, I was unable to find an example of using the XZ argument to manually compute the same thing using the raw draws, and my naive attempts do not seem to match the output of `posterior_epred`

.

```
# Simulate some logistic regression data
set.seed(42)
n_obs <- 1000
num_cat <- 3
rel_min <- 0.1
cat_probs <- runif(num_cat)
cat_probs[cat_probs < rel_min] <- rel_min
cat_probs <- cat_probs / sum(cat_probs)
x <- sample(num_cat, size=n_obs, replace=TRUE, prob=cat_probs)
y_prob <- runif(num_cat)
df <-
data.frame(x=x, p=y_prob[x]) %>%
mutate(y=as.numeric(runif(n()) < p))
num_draws <- 5000
# Run rstanarm
fit <- stan_glmer(y ~ (1 | x),
family = binomial(link = "logit"),
data = df,
prior = normal(0, 1, autoscale = TRUE),
prior_covariance = decov(scale = 0.50),
adapt_delta = 0.99,
refresh = 0,
seed = 1010,
iter = num_draws)
# Then every one of these repeated commands gives a different answer each time:
posterior_linpred(fit, newdata=df, draws=num_draws, seed=42, permuted=FALSE)[1,1]
posterior_linpred(fit, newdata=df, draws=num_draws, seed=42, permuted=FALSE)[1,1]
posterior_linpred(fit, newdata=df, draws=num_draws, seed=42)[1,1]
posterior_linpred(fit, newdata=df, draws=num_draws, seed=42)[1,1]
posterior_linpred(fit, newdata=df, draws=num_draws)[1,1]
posterior_linpred(fit, newdata=df, draws=num_draws)[1,1]
posterior_linpred(fit, newdata=df, draws=num_draws - 10)[1,1]
posterior_linpred(fit, newdata=df, draws=num_draws - 10)[1,1]
```

If possible, add also code to simulate data or attach a (subset of) the dataset you work with.

- Operating System: Ubuntu 22.04.4 LTS
- rstanarm Version: 2.21.4

Don’t forget to add relevant tags to your topic (top right of this form) especially for application area. Delete this text before posting your topic. Thx!

Looking forward to your topic!