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!

  • 1
    $\begingroup$ FYI, if your data aren't sampled densely and evenly over the range of the covariate iv, then the fitted values (at the observed data) won't ever be smooth, even if the model described by the coefficients itself is smooth. $\endgroup$ – Gavin Simpson Jan 30 '15 at 19:34

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.

  • $\begingroup$ @ Gavin. Thank you so much for that excellent answer. Extremely helpful and much appreciated. Can I ask how to predict from a cubic polynomial model when there is a group effect? I tried the following but can't get it to work: MyData = data.frame(iv = seq(from = 75, to = 334, by = 1),d2 = seq(from = 75, to = 334, by = 1),d3 = seq(from = 75, to = 334, by = 1)), Pred = predict(M1, newdata = MyData, level=0, type="response"),plot(…),lines(MyData$iv, Pred). Any advice would be much appreciated. $\endgroup$ – jjulip Jan 31 '15 at 0:12
  • $\begingroup$ @ Gavin. I have now written a separate question on plotting predicted values from polynomial models here: stackoverflow.com/questions/28254403/… as it seems like it might be something that can't be answered within a comment. Thanks for your help with my question above! $\endgroup$ – jjulip Jan 31 '15 at 17:38
  • $\begingroup$ @jjulip not online with computer having R on it, but will take a look on Monday if not before. recall you can get population or individual predictions from lme() though... $\endgroup$ – Gavin Simpson Jan 31 '15 at 17:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.