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.
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.