在使用reformulate之后更新nlme模型

如何解决在使用reformulate之后更新nlme模型

在nlme的公式参数中使用重新格式化后，更新lme()模型存在一个小问题

这是一些数据

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)

现在说我想将变量指定为lme函数之外的对象...

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

...然后重新格式化它们...

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 

到目前为止一切顺利。但是，如果我们想使用update()在模型的固定效果部分添加术语呢？

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

我们收到错误消息

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

很明显，我可以不用使用update()再次写出模型，但是我只是想知道是否有一种方法可以使更新工作。

我认为问题在于使用lme时reformulate对公式参数进行编码。

任何解决方案都值得赞赏。

解决方法

问题在于，当您未在对lme的调用中放入公式文字时，某些类型的函数将不起作用。特别是错误的来源是

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

nlme:::formula.lme试图在错误的环境中评估参数。构建第一个模型的另一种方法是

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

执行此操作时，会将公式注入调用中

formula(mod1)
# score ~ time

这将允许update函数更改公式。