How to predict ln(odds) with rcs term in mixed effects logistic model?

Question

I am using "lme4" package to fit mixed-effects nonlinear logistic model to access the association of Y and X. As the response variable of my data is binary and nlmer function requires response variable to be continuous, I use glmer function and "rms" package function rcs to fit the model and visualize the nonlinear association like the R code below:

library(lme4)
library(rms)
m <- glmer(r2 ~ rcs(Anger, 5) + Gender + situ + btype + (1 | id), 
           data = VerbAgg, family = binomial("logit"),
           control=glmerControl(optimizer="bobyqa"))
p <- predict(m, newdata = VerbAgg, type = "link")
scatter.smooth(VerbAgg$Anger,p,pch='.',col="blue",lpars=list(type="l",col="red"))

I have some questions about using this code:

1. Is the code correct?
2. How to predict the ln(Odds) of r2? Is it "p <- predict(m, newdata = VerbAgg, type = "link")" ?
3. How to visualize the spline of ln(Odds) of r2 and Anger? Is it correct to use "scatter.smooth" function to plot and add a smooth curve in scatters?
4. How to get the P-nonlinear for this model?

Answer 1

1. It looks reasonable.
2. Yes. Just predict(m) should give you the predicted log-odds for each observation (since predict uses the original data by default, and type = "link" is also the default). (But being more explicit doesn't hurt.)
3. You could construct a data frame that has the baseline value for the non-focal predictors and Anger spanning an appropriate range of values, then use that as newdata in a predict call. Or you could use the emmeans or effects packages, both of which try to automate this process for you.

This worked for me (natural splines with ns seem to integrate better with lme4 than rcs)

library(lme4)
library(effects)
m <- glmer(r2 ~ splines::ns(Anger, df = 5) + Gender + situ + btype + (1 | id), 
           data = VerbAgg, family = binomial("logit"),
           control=glmerControl(optimizer="bobyqa"))
plot(Effect("Anger", m)

You can do a significance test of the spline term by using car::Anova() or by model comparison (likelihood ratio test):

m0 <- update(m, r2 ~ Gender + situ + btype + (1 | id))
anova(m, m0)