将nlme::lme中的惩罚样条随机项转换为lme4::lmer兼容形式

我参考《应用纵向分析》中的R代码，用nlme::lme实现了惩罚样条回归分析。由于实际场景需要加入更多随机系数，计划改用lme4::lmer，但无法将代码中bf1至bf22的随机样条项转换为lmer兼容形式。我理解lme实现中是给前缀为bf的分段样条变量指定了随机截距，尝试在gamm4公式中添加s(time,k=22,bs='re')但未成功。

数据使用《应用纵向分析》配套的progesterone数据集，原始nlme代码如下：

library(foreign)
library(nlme)

ds <- read.dta("progesterone.dta")
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

# 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)))

解决方案 #

原模型的核心结构分为两部分：全局分段样条随机效应（对应Const组的22个bf变量）和个体水平的随机截距+时间斜率（对应id组的pdSymm(~time)）。用gamm4结合mgcv和lme4的语法可以完美复刻：

library(gamm4)
library(mgcv)
library(foreign)

# 读取并预处理数据
ds <- read.dta("progesterone.dta")

# 生成分段样条基变量（与原代码一致）
ds$bf1 <- (ds$time+7)*I(ds$time > -7)
ds$bf2 <- (ds$time+6)*I(ds$time > -6)
ds$bf3 <- (ds$time+5)*I(ds$time > -5)
ds$bf4 <- (ds$time+4)*I(ds$time > -4)
ds$bf5 <- (ds$time+3)*I(ds$time > -3)
ds$bf6 <- (ds$time+2)*I(ds$time > -2)
ds$bf7 <- (ds$time+1)*I(ds$time > -1)
ds$bf8 <- (ds$time)*I(ds$time > 0)
ds$bf9 <- (ds$time-1)*I(ds$time > 1)
ds$bf10 <- (ds$time-2)*I(ds$time > 2)
ds$bf11 <- (ds$time-3)*I(ds$time > 3)
ds$bf12 <- (ds$time-4)*I(ds$time > 4)
ds$bf13 <- (ds$time-5)*I(ds$time > 5)
ds$bf14 <- (ds$time-6)*I(ds$time > 6)
ds$bf15 <- (ds$time-7)*I(ds$time > 7)
ds$bf16 <- (ds$time-8)*I(ds$time > 8)
ds$bf17 <- (ds$time-9)*I(ds$time > 9)
ds$bf18 <- (ds$time-10)*I(ds$time > 10)
ds$bf19 <- (ds$time-11)*I(ds$time > 11)
ds$bf20 <- (ds$time-12)*I(ds$time > 12)
ds$bf21 <- (ds$time-13)*I(ds$time > 13)
ds$bf22 <- (ds$time-14)*I(ds$time > 14)

# gamm4模型：结合固定效应、全局随机样条和个体随机效应
model_gamm4 <- gamm4(logp ~ time + group + group:time + group:bf15 +
                       s(time, bs = "re", k = 22),  # 对应原lme的全局随机样条（pdIdent结构）
                     random = ~ (1 + time | id),   # 对应原lme的个体随机截距+时间斜率（pdSymm结构）
                     data = ds)

# 提取lmer格式的模型对象
model_lmer <- model_gamm4$mer

关键说明 #

1. 全局随机样条：s(time, bs = "re", k = 22)对应原模型中Const组的22个bf变量，bs="re"指定随机效应样条，k=22匹配原分段样条的基数量，且默认采用对角协方差结构（与pdIdent一致）。
2. 个体随机效应：random = ~ (1 + time | id)直接复刻原pdSymm(~time)的结构，允许个体截距与时间斜率之间存在协方差。
3. 原代码中错误的s(time,id,k=22,bs="re")写法混淆了分组变量的指定方式，gamm4中需用by参数指定分组，但此处个体水平的随机效应并非样条，而是标准的线性随机斜率，因此用lmer原生语法更准确。

内容的提问来源于stack exchange，提问作者stats134711