R中重复测量数据多层线性模型：时间编码差异引发的结果疑问

在学习《Discovering Statistics Using R》第899页的重复测量多层线性模型案例时，因Time变量编码方式不同，导致模型对比（anova）结果出现明显差异：我采用实际月份编码（0,6,12,18），而教材勘误采用序数编码（0,1,2,3），两者仅在最后一步加入AR(1)相关结构的模型（ARModel）对比中出现显著性差异。以下是可复现问题的代码及两次编码下的anova输出：

library(reshape2)

data_url <- "https://studysites.sagepub.com/dsur/study/DSUR%20Data%20Files/Chapter%2019/Honeymoon%20Period.dat"

satisfactionData <- read.delim(data_url, header = TRUE)

restructuredData <- melt(satisfactionData, 
                         id = c("Person", "Gender"), 
                         measured = c("Satisfaction_Base", "Satisfaction_6_Months", "Satisfaction_12_Months", "Satisfaction_18_Months"))

names(restructuredData) <- c("Person", "Gender", "Time", "Life_Satisfaction")

restructuredData$Time<-as.numeric(restructuredData$Time)-1

# ***************************************************************
# 第二次运行时，取消下方注释并重新执行整个脚本
# ***************************************************************
# restructuredData$Time = restructuredData$Time*6


# 逐步构建模型
intercept <- gls(Life_Satisfaction~1, 
                 data = restructuredData, 
                 method = "ML", 
                 na.action = na.exclude)

randomIntercept <- lme(Life_Satisfaction ~1, 
                       data = restructuredData, 
                       random = ~1|Person, 
                       method = "ML",  
                       na.action = na.exclude, 
                       control = list(opt="optim"))

timeRI <- update(randomIntercept, .~. + Time)

timeRS <- update(timeRI, random = ~Time|Person)

ARModel <- update(timeRS, correlation = corAR1(0, form = ~Time|Person))

anova(intercept, randomIntercept, timeRI, timeRS, ARModel)

序数编码（0,1,2,3）下的anova结果 #

Model df      AIC  BIC         logLik   Test   L.Ratio p-value
intercept           1  2 2064.053 2072.217 -1030.0263                         
randomIntercept     2  3 1991.396 2003.642  -992.6978 1 vs 2  74.65704  <.0001
timeRI              3  4 1871.728 1888.057  -931.8642 2 vs 3 121.66714  <.0001
timeRS              4  6 1874.626 1899.120  -931.3131 3 vs 4   1.10224  0.5763
ARModel             5  7 1872.891 1901.466  -929.4453 4 vs 5   3.73564  0.0533

实际月份编码（0,6,12,18）下的anova结果 #

Model df      AIC  BIC         logLik   Test   L.Ratio p-value
intercept           1  2 2064.053 2072.217 -1030.0263                         
randomIntercept     2  3 1991.396 2003.642  -992.6978 1 vs 2  74.65704  <.0001
timeRI              3  4 1871.728 1888.057  -931.8642 2 vs 3 121.66714  <.0001
timeRS              4  6 1874.627 1899.120  -931.3135 3 vs 4   1.10151  0.5765
ARModel             5  7 1876.627 1905.203  -931.3135 4 vs 5   0.00001  0.9978

问题解析 #

1. 为何出现结果差异 #

核心原因是AR(1)相关结构对Time变量的尺度极度敏感。AR(1)模型假设相邻时间点的残差相关系数为$\rho$，这里的“相邻”由Time变量的数值差定义：

• 序数编码下，相邻时间点的数值差为1，模型估计的是间隔1个单位（对应实际6个月）的残差相关性；
• 实际月份编码下，相邻时间点的数值差为6，此时AR(1)模型的相关结构变为$\rho^{6$（AR(1)的k步相关是$\rho}k$）。

案例中实际时间间隔固定为6个月，用月份编码时，模型尝试估计的是间隔6个单位（对应实际36个月）的相关性，完全偏离了数据的真实时间结构，最终导致ARModel的拟合效果几乎没有提升（似然值几乎不变），对比timeRS模型时似然比趋近于0。

2. 为何差异仅在最后一步显现 #

前面的模型（intercept、randomIntercept、timeRI、timeRS）未引入时间相关结构，仅涉及固定效应、随机截距/斜率：

• 固定效应中Time的系数会随尺度线性缩放（比如月份编码下的系数是序数编码的1/6），但模型的拟合优度和似然比检验结果不受线性缩放影响——线性变换仅改变系数尺度，不改变模型对数据的解释能力；
• 随机效应部分的方差也会随Time尺度缩放，但同样不影响似然比检验的结果。

只有当引入AR(1)这种依赖时间间隔尺度的相关结构时，编码方式的差异才会直接改变模型对残差相关性的假设，进而影响模型拟合效果和对比结果。

3. 实际应用中月份编码的意义与注意事项 #

意义

• 保持解释性：用真实月份编码能让固定效应系数直接反映“每过1个月”的因变量变化，结果解释更直观；
• 适配后续分析：如果后续要加入其他时间相关协变量（如实际日期变量），统一尺度可避免转换带来的麻烦。

注意事项

• 匹配时间相关结构的尺度：使用AR、MA等时间序列相关结构时，必须保证Time变量的尺度与实际时间间隔一致。比如案例中实际间隔是6个月，若用月份编码，应将Time变量转换为以“6个月”为单位的序数，或调整相关结构的参数；
• 注意系数解释的尺度：序数编码下的系数是“每增加1个时间单位（6个月）”的变化，而月份编码下是“每增加1个月”的变化，解读时需明确尺度；
• 模型对比的一致性：对比涉及时间相关结构的模型时，必须确保Time变量的编码方式一致，否则对比结果无意义。

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