Updating nlme models after using reformulate()

Question

Having a small issue with updating nlme models after using reformulate in the formula argument of lme()

Here is some data

set.seed(345)
A0 <- rnorm(4,2,.5)
B0 <- rnorm(4,2+3,.5)
A1 <- rnorm(4,6,.5)
B1 <- rnorm(4,6+2,.5)
A2 <- rnorm(4,10,.5)
B2 <- rnorm(4,10+1,.5)
A3 <- rnorm(4,14,.5)
B3 <- rnorm(4,14+0,.5)
score <- c(A0,B0,A1,B1,A2,B2,A3,B3)
id <- rep(1:8,times = 4, length = 32)
time <- factor(rep(0:3, each = 8, length = 32))
group <- factor(rep(c("A","B"), times =2, each = 4, length = 32))
df <- data.frame(id = id, group = group, time = time,  score = score)

Now say I want to specify the variables as objects outside the lme function...

t <- "time"
g <- "group"
dv <- "score"

...and then reformulate them...

mod1 <- lme(fixed = reformulate(t, response = "score"),
            random = ~1|id, 
            data = df)

summary(mod1)

Linear mixed-effects model fit by REML
 Data: df 
       AIC      BIC    logLik
  101.1173 109.1105 -44.55864

Random effects:
 Formula: ~1 | id
        (Intercept)  Residual
StdDev:   0.5574872 0.9138857

Fixed effects: reformulate(t, response = "score") 
                Value Std.Error DF   t-value p-value
(Intercept)  3.410345 0.3784804 21  9.010626       0
time1        3.771009 0.4569429 21  8.252693       0
time2        6.990972 0.4569429 21 15.299445       0
time3       10.469034 0.4569429 21 22.911036       0
 Correlation: 
      (Intr) time1  time2 
time1 -0.604              
time2 -0.604  0.500       
time3 -0.604  0.500  0.500

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.6284111 -0.5463271  0.1020036  0.5387158  2.1784156 

Number of Observations: 32
Number of Groups: 8

So far so good. But what if we want to add terms to the fixed effects portion of the model using update()?

mod2 <- update(mod1, reformulate(paste(g,"*",t), response = "score"))

We get the error message

Error in reformulate(t, response = "score") : 
  'termlabels' must be a character vector of length at least one

Obviously I can write the model out again without using update() but I was just wondering if there is a way to make update work.

I gather the problem lies in the way that lme encodes the formula argument when using reformulate.

Any solution much appreciated.

score 2 · Accepted Answer · answered Sep 03 '20 at 23:31

The problem is that when you don't put in formula literal in the call to lme, certain types of functions don't work. In particular, the place where the error is coming from is

formula(mod1)
#  Error in reformulate(t, response = "score") : 
#  'termlabels' must be a character vector of length at least one

The nlme:::formula.lme tries to evaluate the parameter in the wrong environment. A different way to construct the first model would be

mod1 <- do.call("lme", list(
  fixed = reformulate(t, response = "score"),
  random = ~1|id, 
  data = quote(df)))

When you do this, this injects the formula into the call

formula(mod1)
# score ~ time

which will allow the update function to change the formula.

Great solve @MrFlick. And thank you for the explanation. `do.call` is one of those functions which I have seen advanced R people use but have never quite understood. Why does `do.call` result in the rendering of the formula in a way that can be understood downstream but the original does not? What does it do? — llewmills, Sep 04 '20 at 04:32

Updating nlme models after using reformulate()

1 Answers1