R与SAS混合模型AIC计算结果不一致，求教SAS的AIC计算公式

我最近在用R的nlme包复现SAS中的重复测量双因素方差分析，模型采用复合对称(CS)协方差矩阵和REML估计方法。目前除了AIC和BIC值，其他结果（包括对数似然值、ANOVA表）都和SAS完全一致，所以想搞清楚SAS的AIC到底是怎么计算的。

我的R操作代码 #

library(nlme)
dataset_melt <- structure(list(Groupe = c("A", "A", "A", "A", "A", "B", "B", "B", "B", "B", "C", "C", "C", "C", "C", "A", "A", "A", "A", "A", "B", "B", "B", "B", "B", "C", "C", "C", "C", "C", "A", "A", "A", "A", "A", "B", "B", "B", "B", "B", "C", "C", "C", "C", "C", "A", "A", "A", "A", "A", "B", "B", "B", "B", "B", "C", "C", "C", "C", "C", "A", "A", "A", "A", "A", "B", "B", "B", "B", "B", "C", "C", "C", "C", "C"), ID = c("01/001", "01/002", "01/003", "01/004", "01/005", "02/001", "02/002", "02/003", "02/004", "02/005", "03/001", "03/002", "03/003", "03/004", "03/005", "01/001", "01/002", "01/003", "01/004", "01/005", "02/001", "02/002", "02/003", "02/004", "02/005", "03/001", "03/002", "03/003", "03/004", "03/005", "01/001", "01/002", "01/003", "01/004", "01/005", "02/001", "02/002", "02/003", "02/004", "02/005", "03/001", "03/002", "03/003", "03/004", "03/005", "01/001", "01/002", "01/003", "01/004", "01/005", "02/001", "02/002", "02/003", "02/004", "02/005", "03/001", "03/002", "03/003", "03/004", "03/005", "01/001", "01/002", "01/003", "01/004", "01/005", "02/001", "02/002", "02/003", "02/004", "02/005", "03/001", "03/002", "03/003", "03/004", "03/005"), temps = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L), .Label = c("T0", "T1", "T2", "T3", "T4"), class = "factor"), value = c(29.4, 21, 23.4, 26.2, 28.5, 27.8, 27.2, 20.6, 20.2, 25.3, 26.2, 29.2, 27.1, 23.1, 20.6, 22.9, 29.6, 20.9, 25.2, 25, 26, 26.7, 25.1, 21, 28.2, 23.4, 27.1, 29.8, 22.2, 26.6, 29.9, 29.1, 23.4, 22.6, 25.7, 24.5, 29.6, 21.5, 28.9, 20.1, 26.5, 23.4, 24.9, 25.3, 25, 27.4, 29.5, 24.6, 27.4, 24.6, 21.3, 23.6, 22.8, 23.6, 20.6, 26.5, 29.2, 20.6, 25.7, 29.1, 23.7, 24.3, 28.7, 21.9, 23.7, 29.8, 27.1, 28.7, 28.3, 20.4, 28.7, 20.3, 22.8, 23.4, 21.5)), row.names = c(NA, -75L), .Names = c("Groupe", "ID", "temps", "value"), class = "data.frame")
options(contrasts=c("contr.SAS","contr.poly"))
mon_lme <- lme(value ~ Groupe *temps, random = ~ +1 | ID, correlation=corCompSymm(form=~temps|ID), data = dataset_melt,method='REML')
anova(mon_lme) # 与SAS结果基本一致
# numDF denDF F-value p-value
# (Intercept) 1 48 6040.352 <.0001
# Groupe 2 12 0.495 0.6215
# temps 4 48 0.057 0.9938
# Groupe:temps 8 48 1.175 0.3334
summary(mon_lme)$AIC # 363.938
summary(mon_lme)$BIC # 399.5419
k <- attr(logLik(mon_lme), "df")
aic <- 2 * k -2 * logLik(mon_lme)
aic -2 * logLik(mon_lme) # 与SAS结果一致
# 'log Lik.' 329.6698 (df=18)

关键差异：AIC的参数计数方式 #

AIC的核心公式是 $AIC = -2 \times logLik + 2k$，其中k是模型中估计的参数数量。R和SAS的分歧就出在k的计数范围上：

• R的nlme包：计算AIC时，k包含所有估计参数——包括固定效应参数、随机效应方差参数、协方差结构参数。比如我的模型中，logLik(mon_lme)的自由度是18，这个数字包含了Groupe*temps的15个固定效应参数，加上随机截距方差、CS相关系数、残差方差这3个协方差参数，总共18个。所以R的AIC是 $218 - 2logLik$。

• SAS的REML模型：计算AIC时，k只统计协方差参数（随机效应和残差相关/方差参数）的数量，不包含固定效应参数。原因是REML估计中，固定效应是被“吸收”的，SAS默认在模型选择指标中只考虑协方差结构的参数数量。

验证我的情况 #

在我的代码中，手动调整k为协方差参数数量（3个）后，计算出的AIC就和SAS完全一致了：

# 假设协方差参数数量为3
sas_aic <- -2 * logLik(mon_lme) + 2*3

这个结果会和SAS输出的AIC匹配，而R默认的AIC因为多算了15个固定效应参数，所以数值会比SAS的大很多。

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