使用 reformulate() 后更新 nlme 模型

Question

在 lme()

的公式参数中使用 reformulate 后更新 nlme 模型时出现小问题

这是一些数据

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 时对公式参数进行编码的方式。

非常感谢任何解决方案。

Answer 1

问题是，当您不在 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 函数更改公式。

使用 reformulate() 后更新 nlme 模型

Updating nlme models after using reformulate()

r

nlme