restricted cubic splines with emmeans

[nbarrowman]

I am trying to use emmeans with an lmer model with a restricted cubic spline term using rcs but I get the following error message:

Error in rcspline.eval(x, nk = nknots, inclx = TRUE, pc = pc, fractied = fractied) : 
  knots not specified, and < 6 non-missing observations

In the "Extending emmeans" vignette there is a brief mention of splines:

The params argument is needed when the model formula refers to variables besides predictors. For example, a model may include a spline term, and the knots are saved in the user's environment as a vector and referred to in the call to fit the model. In trying to recover the data, we try to construct a data frame containing all the variables present on the right-hand side of the model, but if some of those are scalars or of different lengths than the number of observations, an error occurs. So you need to exclude any names in params when reconstructing the data.

I feel like I'm close, but can't quite figure out how to proceed.

Thanks, I love the package!

[rvlenth]

It would help a lot to see the R code you used to fit the model. But the main point is that the names of all local variables in that model-fitting call need to be included in params.

[rvlenth]

With some effort I found out where rcs is (in rms package), but the example that accompanies it is disabled and does not work if I try to run it.

Here's an example using the ns() function:

> my.knots = c(2.5, 3, 3.5)
> mod = lm(Sepal.Length ~ Species*ns(Sepal.Width, knots = my.knots), data = iris)

We get an error if we don't specify params:

> ref_grid(mod)
Error in ref_grid(mod) : We are unable to reconstruct the data.
The variables needed are:
	Species Sepal.Width my.knots
Are any of these actually constants? (specify via 'params = ')
The dataset name is:
	iris
Does the data still exist? Or you can specify a dataset via 'data = '

The clue here is that it's looking for a variable named my.knots to include in the reference grid.
If we specify params = "my.knots", it knows that my.knots is not a predictor, but rather some other variable needed by the model object:

> ref_grid(mod, params = "my.knots")
'emmGrid' object with variables:
    Species = setosa, versicolor, virginica
    Sepal.Width = 3.0573

Now here are some results at different x values:

> emmeans(mod, ~ Species | Sepal.Width, params = "my.knots", 
+     at = list(Sepal.Width = 2:4))
NOTE: Results may be misleading due to involvement in interactions
Sepal.Width = 2:
 Species    emmean     SE  df lower.CL upper.CL
 setosa       5.09 2.1391 135    0.861     9.32
 versicolor   5.48 0.3392 135    4.812     6.15
 virginica    6.18 0.7918 135    4.613     7.74

Sepal.Width = 3:
 Species    emmean     SE  df lower.CL upper.CL
 setosa       4.65 0.1366 135    4.381     4.92
 versicolor   6.20 0.0951 135    6.011     6.39
 virginica    6.73 0.0934 135    6.547     6.92

Sepal.Width = 4:
 Species    emmean     SE  df lower.CL upper.CL
 setosa       5.41 0.1269 135    5.159     5.66
 versicolor  -5.77 9.3033 135  -24.172    12.63
 virginica    9.20 0.6478 135    7.916    10.48

Confidence level used: 0.95 

I hope this clarifies the purpose and use of params

[nbarrowman]

That's super helpful! Thank you.

I had prepared a reproducible example, but I'll read this carefully before proceeding.

[nbarrowman]

I'm happy to report that I got the rcs function (for restricted cubic splines) to work with emmeans. (Yes, rcs is from the rms package, which I should have stated.) Here's a reproducible example using a made-up data set in case someone else wants to do this.

Note that rcs can automatically choose knots. In the code below I first ask for and extract the locations of 3 knots using my.knots <- attributes(rcs(d$months,3))$parms. Note that in rcs it's parms not params! Then I call lmer with a restricted cubic spline term specified as rcs(months,parms=my.knots). Finally, I call ref_grid and I specify params="my.knots". (Just to underline the point, this time it's params not parms).

d <- structure(list(id = c(1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 4, 4, 
4, 5, 5, 5, 5), months = c(6, 12, 18, 24, 5, 15, 20, 26, 7, 12.5, 
19, 5.5, 11.5, 22.5, 6.2, 12.2, 18.2, 24.2), outcome = c(2, 4, 
5, 6, 3, 4, 4, 7, 2, 8, 5, 10, 11, 12, 1.8, 3.8, 5.2, 6.1)), row.names = c(NA, 
-18L), class = c("tbl_df", "tbl", "data.frame"))

library(lme4)
library(rms)
library(emmeans)

my.knots <- attributes(rcs(d$months,3))$parms
fit <- lmer(outcome~rcs(months,parms=my.knots)+(1|id),data=d)

rg <- ref_grid(fit,at=list(months=seq(0,32,by=0.1)),params="my.knots")

p1 <- emmip(rg,~months,CIs=TRUE,plotit=FALSE)
ggplot(data=p1,aes(x=months,y=yvar))+
  geom_line()+
  geom_ribbon(data=p1,aes(ymin=LCL,ymax=UCL),alpha=0.3) +
  xlab("Time (months)")+ylab("Outcome")

Many thanks for the help and for the emmeans package!

[rvlenth]

Terrific! I'm glad you got it to work. Plus, this is a good example; and the care with which you emphasize params vs. parms is useful.

I am realizing that there isn't an example in the vignettes of using params. Do you have an objection to linking this issue somewhere in the vignettes or help pages?

[nbarrowman]

No objection, cheers!

[rvlenth]

I think we've completed this issue, so closing. But comment further if something comes up.

[befriendabacterium]

Hi both - thanks for this thread, it helped me troubleshoot my own issue using emmeans with a model with an rcs term. I just wanted to add to this to share knowledge based on my troubleshooting of a similar issue, in case it helps someone who comes across this thread.

I've found is that, in rms::rcs, if you specify the 'parms' argument as a number and allow the rms::rcs() function to define the knots, it makes your model incompatible with emmeans functions such as emmeans, qdrg and ref_grid, because emmeans can't find the knots parameter (even though they are stored in the 'parms' object of the model. So you have to first make an object with the knots (using the same procedure as the rms::rcs() function would use within it) and then use this as the 'parms' argument in the rms::rcs bit of your model. Below are examples of the same model, but in emmeans incompatible and compatible formats, respectively (based on @nbarrowman's examples above):

#load data
d <- structure(list(id = c(1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 4, 4, 
4, 5, 5, 5, 5), months = c(6, 12, 18, 24, 5, 15, 20, 26, 7, 12.5, 
19, 5.5, 11.5, 22.5, 6.2, 12.2, 18.2, 24.2), outcome = c(2, 4, 
5, 6, 3, 4, 4, 7, 2, 8, 5, 10, 11, 12, 1.8, 3.8, 5.2, 6.1)), row.names = c(NA, 
-18L), class = c("tbl_df", "tbl", "data.frame"))

library(lme4)
library(rms)
library(emmeans)

#emmeans incompatible example - throws up error 'Error in rcspline.eval(x, nk = nknots, inclx = TRUE, pc = pc, fractied = fractied) : knots not specified, and < 6 non-missing observation'
my.knots<-3
fit <- lmer(outcome~rcs(months,parms=my.knots)+(1|id),data=d)
rg <- ref_grid(fit,at=list(months=seq(0,32,by=0.1)),params="my.knots")

#emmeans compatible example - works fine
my.knots <- attributes(rcs(d$months,3))$parms
fit <- lmer(outcome~rcs(months,parms=my.knots)+(1|id),data=d)
rg <- ref_grid(fit,at=list(months=seq(0,32,by=0.1)),params="my.knots")

[rvlenth]

Thanks. Glad this was helpful