如何在R的lcmm包hlme中移除潜类别增长模型的非显著高阶项

问题背景 #

我使用R语言lcmm包的hlme函数拟合潜类别增长模型，初始代码如下：

gridsearch(rep = 100, maxiter = 10, minit = lcga1_c,
 m=hlme(cHairCortisolConclog~Visit+I(Visit^2)+I(Visit^3), subject = "ID", ng = 3, 
data = data_long, mixture = ~ Visit+I(Visit^2)+I(Visit^3)))

模型纳入了Visit的三次项，但输出显示仅类别2和3的二次、三次项显著，类别1的对应项无统计学意义：

Fixed effects in the longitudinal model:

coef      Se   Wald p-value

intercept class1 -2.87572 2.32800 -1.235 0.21673
intercept class2 1.61864 0.52613 3.077 0.00209
intercept class3 2.34416 0.49703 4.716 0.00000
Visit class1 4.00429 2.80767 1.426 0.15381
Visit class2 1.30048 0.70896 1.834 0.06660
Visit class3 1.46184 0.64264 2.275 0.02292
I(Visit^2) class1 -1.34900 1.01259 -1.332 0.18279
I(Visit^2) class2 -0.76651 0.27105 -2.828 0.00469
I(Visit^2) class3 -0.75885 0.24353 -3.116 0.00183
I(Visit^3) class1 0.14569 0.11131 1.309 0.19060
I(Visit^3) class2 0.09950 0.03055 3.257 0.00112
I(Visit^3) class3 0.09487 0.02730 3.475 0.00051

希望仅移除类别1的二次项和三次项，得到如下目标输出：

Fixed effects in the longitudinal model:

coef      Se   Wald p-value

intercept class1 -2.87572 2.32800 -1.235 0.21673
intercept class2 1.61864 0.52613 3.077 0.00209
intercept class3 2.34416 0.49703 4.716 0.00000
Visit class1 4.00429 2.80767 1.426 0.15381
Visit class2 1.30048 0.70896 1.834 0.06660
Visit class3 1.46184 0.64264 2.275 0.02292
I(Visit^2) class2 -0.76651 0.27105 -2.828 0.00469
I(Visit^2) class3 -0.75885 0.24353 -3.116 0.00183
I(Visit^3) class2 0.09950 0.03055 3.257 0.00112
I(Visit^3) class3 0.09487 0.02730 3.475 0.00051

解决方案 #

lcmm包没有像SAS PROC TRAJ那样的ORDER=语句，但可以通过构造类别与轨迹项的交互项，将特定类别的二次、三次项系数约束为0，具体实现如下：

步骤1：提取或生成类别归属变量 #

从初始全模型中提取每个个体的类别归属，生成因子变量（以class1为参照）：

# 从全模型lcga1_c中提取类别，生成因子变量
data_long$class <- factor(lcga1_c$pprob$class, levels = c(1, 2, 3))

步骤2：拟合带约束的模型 #

调整hlme的公式，让二次项和三次项仅与类别变量交互（class1的对应项系数会被自动设为0），同时保留所有类别的截距和Visit项：

gridsearch(rep = 100, maxiter = 10, minit = lcga1_c,
           m=hlme(cHairCortisolConclog ~ Visit + I(Visit^2):class + I(Visit^3):class,
                  subject = "ID", ng = 3,
                  data = data_long,
                  mixture = ~ Visit + I(Visit^2):class + I(Visit^3):class))

替代写法（手动构造哑变量） #

如果偏好手动控制哑变量，也可以生成class2和class3的哑变量，再与二次、三次项交互：

# 生成class2和class3的哑变量
data_long$class2 <- as.integer(lcga1_c$pprob$class == 2)
data_long$class3 <- as.integer(lcga1_c$pprob$class == 3)

# 拟合模型
gridsearch(rep = 100, maxiter = 10, minit = lcga1_c,
           m=hlme(cHairCortisolConclog ~ Visit + I(Visit^2)*class2 + I(Visit^2)*class3 + I(Visit^3)*class2 + I(Visit^3)*class3,
                  subject = "ID", ng = 3,
                  data = data_long,
                  mixture = ~ Visit + I(Visit^2)*class2 + I(Visit^2)*class3 + I(Visit^3)*class2 + I(Visit^3)*class3))

以上两种方法都能实现仅保留class2和class3的二次、三次项，class1的对应项系数为0，与目标输出一致。

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