使用nlme::lme将代码转换为lme 4：：lmer

Question

我正在遵循R码(这里)在应用纵向分析由Fitzmaurice，Laird和Ware进行惩罚样条回归分析使用nlme::lme。我的实际应用程序将有更多的随机系数，所以我计划使用lme4::lmer。我很难将链接代码中的随机样条术语bf1-bf22转换为在lmer中兼容的。如果我理解lme实现，就会为分段样条(前缀为bf的变量)指定随机截取。我(天真地)尝试将s(time,k=22,bs='re')添加到lmer公式的RHS中(下面的代码显示了使用lme4函数的gamm4 )，但这似乎行不通。

数据可以从这里下载，代码如下所示。

library(foreign)

ds <- read.dta("progesterone.dta") #Wherever the file is stored)
attach(ds)

bf1 <- (time+7)*I(time > -7)
bf2 <- (time+6)*I(time > -6)
bf3 <- (time+5)*I(time > -5)
bf4 <- (time+4)*I(time > -4)
bf5 <- (time+3)*I(time > -3)
bf6 <- (time+2)*I(time > -2)
bf7 <- (time+1)*I(time > -1)
bf8 <- (time)*I(time > 0)
bf9 <- (time-1)*I(time > 1)
bf10 <- (time-2)*I(time > 2)
bf11 <- (time-3)*I(time > 3)
bf12 <- (time-4)*I(time > 4)
bf13 <- (time-5)*I(time > 5)
bf14 <- (time-6)*I(time > 6)
bf15 <- (time-7)*I(time > 7)
bf16 <- (time-8)*I(time > 8)
bf17 <- (time-9)*I(time > 9)
bf18 <- (time-10)*I(time > 10)
bf19 <- (time-11)*I(time > 11)
bf20 <- (time-12)*I(time > 12)
bf21 <- (time-13)*I(time > 13)
bf22 <- (time-14)*I(time > 14)


Const <- factor(rep(1,length(logp)))

group.time <- group*time
group.bf15 <- group*bf15
require(nlme)
model_lme <- lme(logp ~ time + group + group.time + group.bf15, 
              random=list(Const=pdIdent(~-1 + bf1 + bf2 + bf3 + bf4 + bf5 + bf6 + 
                                        bf7 + bf8 + bf9 + bf10 + bf11 + bf12 + bf13 + bf14 + bf15 + bf16 + 
                                        bf17 + bf18 + bf19 + bf20 + bf21 + bf22), 
                          id=pdSymm(~time)))


require(gamm4)
require(mgcv)
model_lmer <- gamm4(logp ~ time + group + group.time + group.bf15 +
                   s(time,k=22,bs="re") + s(time,id,k=22,bs="re"))

Answer 1

装载包

library(foreign)
library(nlme)
library(gamm4)
library(mgcv)
library(ggplot2); theme_set(theme_bw())
library(cowplot)

过程数据

我简化了创建单个线性样条基函数的过程. (attach()几乎不是一个好主意.)

ds <- read.dta("progesterone.dta") #Wherever the file is stored)
dss <- within(ds, {
    for (i in (-7):14) {
        assign(paste0("bf",i+8), (time+i)*as.numeric(time > i))
    }
})
matplot(dss[dss$id==1,5:26], type ="l")
dss$Const <- factor(rep(1,nrow(dss)))

线性样条拟合lme模型

同样，精简(使用reformulate而不是写出所有bf1 + bf2 + ...术语；使用time*group而不是手工构建所有虚拟/交互)。

model_lme <- lme(logp ~ time*group + I(group*bf15), 
                 random=list(Const=pdIdent(
                                 reformulate(c("-1",
                                               paste0("bf",1:22)))),
                             id=pdSymm(~time)),
                 data = dss)

gamm4版本

ds$bf15 <- dss$bf15
model_gamm4 <- gamm4(logp ~time*group + I(group*bf15) +
                         s(time, bs = "cs"),
                     ## need parentheses!
                     random = ~(time|id),
                     data = ds)

这里需要注意的几点是：(1)样条项的集合用s(time, bs = "cs")表示(其中"cs“表示”三次样条“：参见?smooth.construct.cs.smooth.spec代替'cs‘来查看对其他平滑选择的帮助)；(2)原始模型中的另一个RE项是关于时间的随机斜率项(id = pdSymm(~time))，在这里表示为一个单独的random = (time|id)；这个组件是在lmer-fitting侧处理的，而不是为它构造mgcv基。

比较预测

dss2 <- data.frame(ds, pred1 = predict(model_lme), pred2 = predict(model_gamm4$mer))
gg0 <- ggplot(dss, aes(time, logp, colour = factor(group))) + geom_point()
plot_grid(gg0 + geom_line(data=dss2, aes(group=id, y = pred1)) + 
                 labs(title = "lme (linear)"),
   gg0 + geom_line(data=dss2, aes(group=id, y = pred2)) + 
                 labs(title = "gamm4 (cubic)"))

​

​

我也在mgcv上试过这个

model_gam <- gam(logp ~time*group + I(group*bf15) +
                     s(time, bs = "cs") + s(id, time, bs = "re") +
                     s(id, bs = "re"),
                 data = dss, method = "REML")

唯一明显的区别是，斜率和截距区域是独立指定的，而不是相互关联的，但预测值有很大不同；我不知道这是关于模型是否合适，还是在预测过程中如何处理随机效应。