I’m fitting a mixed effect model in rstanarm of the form

```
m1<-stan_glmer(y~gender+country+ rcs(age) + rcs(AnotherContinousVar) + (1|subjectID),family="binomial")
```

I’d like to plot the relationship between the restricted cubic splines of each continous var VS the predicted probabilities. But I’m not sure I’m doing it correctly. Here’s what I do

First, create predictions on the same data using posterior_linpred(m1,transform=TRUE) which gives me predicted probabilites.

Next, I take the mean of these probabilities for each data point and store them in the original df so as to match the original datapoints

```
preds<-posterior_linpred(m1,transform=TRUE)
pred<-colMeans(preds)
dta%>%
mutate(probability=pred)%>%
ggplot(.,aes(Age,probability))+geom_smooth(method = "loess")
```

What I’m confused about is

- Is this the probability given Age only or is this the Pr given Age AND all other factors ? If so, what level are those factors at? Are these marginal effects or conditional?
- Is loess smoothing the appropriate way of visualizing the relationship between the continous var and the predicted probabilities?
- Is there another way that I can visualize the relationship and the uncertainty of the continous var vs probabilities but also, how do I interpret them.
- Should I use something like
`r expand.grid()`

to create all possible combinations and then get predicted probabilities on those data points instead?

Another thing I’m still not clear about is, subjects are treated as random effect with a separate intercept for each (mutliple response by each subject). Do the predicted probabilities take into account the variance explained by these random effects?