Using B-splines within a linear mixed-effects model in R

Question

I am using linear mixed-effect model (run with the lme() function in the nlme package in R) that has one fixed effect, one random intercept term (to account for groups) and an AR(1) correlation structure to account for temporal autocorrelation. The model is a cubic polynomial model specified as so:

d2 = df$iv^2
d3 = df$iv^3

M1 = lme(dv ~ iv + d2 + d3, data=df, random= ~1|group, method="ML",
         correlation = corAR1(form =~iv|group))

Unfortunately, the fitted values from this model do not give me a nice smooth curve (instead I get what looks like a number of joined straight-line segments).

I have searched for a way to make the curve look smoother (less disjointed) and have found that the bs() function returns slightly smoother curves. From my understand, the bs function fits polynomial splines. I am firstly wondering whether the use of B-splines from bs() is technically correct within a linear mixed-effects model? If so, why do the fitted values differ between the two approaches when bs also seems to use cubic polynomials as a default?

M2 = lme(dv ~ bs(iv, df=5), data=df, random= ~1|group, method="ML",
         correlation = corAR1(form =~iv|group))

Any advice would be much appreciated!

Answer 1

Q1: Yes, using a basis expansion such as a B-spline basis is acceptable in a linear mixed model. The model is really a linear fit to the response over a set of derived variables rather than the actual data, in much the same way as you created two new variables from iv, its square and cube, and fitted a model for those instead of just iv itself.

Q2: That said, these are two very different models. The cubic polynomial is a global model, the best fitting cubic polynomial over the range of the data. The model using the B-spline basis is a piecewise polynomial of some degree $k$, say $k = 3$ for cubic splines. The individual pieces are fitted over portions of the data, the boundaries between portions are interior knots. Some constraints on the basis functions force the joins at the knots to be smooth. The key point however is that the B-spline model fits polynomials over local sections of the covariate iv whereas the polynomial model fits a global cub polynomial model over the entire range of iv.

In both cases, the fitted curve should be smooth. What you may be seeing is the effect of irregular sampling over the range of iv. The fitted values will be for the observed values of iv and then a line is drawn to join up these fitted values. If there are varying gaps in the (ordered) iv data, this can show up as low smoothness. If so, and what I tend to do by default regardless, is to predict from the model for $N$ new iv points spaced regularly over the range of iv. You can make $N$ as large as you need to achieve a sufficiently smooth representation of the fitted model.