问题描述
在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函数更改公式。