如何在R中进行带有多个断点的分段线性混合回归？

Question

我正在 R 中拟合分段线性混合回归。我知道我可以使用

lme

包中的

nlme

，然后使用

segmented

来执行分段线性混合回归。然而，在阅读了

segmented

包的文档后，我注意到

segmented.lme

只能处理1个断点，其中我有两个（在第30天和第90天）。

作为背景，我想对第 0、30、90 和 180 天的汽车里程 (

macars

)（变量

days

）进行建模，并以

age

作为混杂因素。请注意，该模型只是说明性的，而不是真实数据。

这是我的原型代码，它使用

lme

包，在我读到

segmented.lme

只能处理 1 个断点后感到困惑之前：

fit <- lme(macars ~ days + age, random = ~ days | id, data = df)
summary(fit)
pw.fit <- segmented(fit, seg.Z = ~ days, psi = list(days = c(30, 90)),  random = list(id = pdDiag(~1 + days))
summary(pw.fit)

编辑： 根据@user2554330提供的见解，我设法拟合模型如下：

> fit <- lme(macars ~ bs(days, knots = c(30, 90), degree = 1) + age, random = ~ days | id, data = df)
> summary(fit)
Linear mixed-effects model fit by REML

Random effects:
 Formula: ~days | id
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev      Corr  
(Intercept) 0.165393834 (Intr)
days        0.001133132 -0.222
Residual    0.292970477       

Fixed effects:  macars ~ bs(days, knots = c(30, 90), degree = 1) + age
                                                 Value  Std.Error  DF    t-value p-value
(Intercept)                                   3.370401 0.13087013 125  25.753787  0.0000
bs(days, knots = c(30, 90), degree = 1)1     -0.883785 0.07340094 125 -12.040518  0.0000
bs(days, knots = c(30, 90), degree = 1)2     -0.870973 0.11990249 125  -7.264013  0.0000
bs(days, knots = c(30, 90), degree = 1)3     -0.722164 0.10003216 125  -7.219320  0.0000
age                                           0.008423 0.00331230  60   2.542882  0.0136

Correlation: 
                                            (Intr)   b(days,k=c(30,90),d=1)1 b(days,k=c(30,90),d=1)2 b(days,k=c(30,90),d=1)3
bs(days, knots = c(30, 90), degree = 1)1     -0.305                                                                               
bs(days, knots = c(30, 90), degree = 1)2     -0.237  0.327                                                                        
bs(days, knots = c(30, 90), degree = 1)3     -0.221  0.465                   0.264                                                
age                                          -0.898 -0.005                   0.090                  -0.019                        

Number of Observations: 100
Number of Groups: 30

现在的问题是如何解释这些值？根据下图，我预计

bs(days, knots = c(30, 90), degree = 1)1

高度为负值，

bs(days, knots = c(30, 90), degree = 1)2

为轻微负值，

bs(days, knots = c(30, 90), degree = 1)3

为轻微正值，但这里的情况并非如此。有什么遗漏吗？

提前致谢

Answer 1

如果您知道断点在哪里，按照 @user255430 的建议，您可以通过 bs(days, knots = c(30, 90), degree = 1) 构建一个

线性 B 样条基础

（或者更具体地要求模型函数为您构建一个）。

但是，这并不像您想象的那样参数化。

library(splines)
days <- 1:200
X <- bs(days, knots = c(30, 90), degree = 1)
par(las = 1, bty = "l") ## cosmetic
matplot(X, type = "l")

如您所见，最后一部分 (

days > 90

) 是通过第二个（红色虚线）和第三个（绿色点线）分量的总和来预测的，而不仅仅是第三个分量。

您可以使用截断幂基础样条线来代替：

library(cplm)q
## set k=3 to suppress warning
p <- tp(days, knots = c(30, 90), k = 3, degree = 1)
X <- cbind(p$X, p$Z)
matplot(X, type = "l")

但是，这有点不方便；对于只有两个结，您可以通过

将组件包含在模型中

~ ... days + I(days*(days>30)) + I(days*(days>90))

一般来说，您至少应该考虑将所有这些术语也包含在随机效应组件中。