R中estimate_means计算EMM时遗漏分组的问题排查

问题背景 #

我需要为4所学校、在校时长不同的儿童群体计算得分realscore的Estimated Marginal Means（EMM，边际均值估计），但复杂模型的结果中始终缺失School C的数据。

数据样例 #

ChildTIS ChildSex   realscore   SchoolCode
8.02        2               7           B
11.04       1               6           C

复杂交互模型下的EMM结果缺失 #

使用包含三因素交互项的线性模型计算EMM，School C的结果完全空白：

# 基于学校、在校时长、性别的EMM计算
realmodel <- lm(realscore ~ SchoolCode * ChildTIS * ChildSex, data = gcsubsetdata)
means_complex_real <- estimate_means(realmodel, at = "SchoolCode")
means_complex_real

输出结果：

SchoolCode | Mean |   SE |       95% CI
---------------------------------------
A          | 0.94 | 0.20 | [0.53, 1.35]
B          | 0.68 | 0.29 | [0.10, 1.25]
C          |      |      |             
D          | 1.79 | 0.29 | [1.22, 2.37]

无效的尝试方案 #

• 数值编码学校ID：将学校转为1、2、3、4的数值型变量，能得到完整结果，但均值与分类编码时不一致，而学校属于分类变量，不能用数值编码替代：

realmodel <- lm(realscore ~ SchoolID * ChildTIS * ChildSex, data = gcsubsetdata)
means_complex_real <- estimate_means(realmodel, at = list(SchoolID = c(1,2,3,4)))
means_complex_real

输出：

SchoolID | Mean |   SE |       95% CI
-------------------------------------
1.00     | 0.83 | 0.19 | [0.45, 1.20]
2.00     | 1.10 | 0.13 | [0.84, 1.37]
3.00     | 1.38 | 0.16 | [1.06, 1.70]
4.00     | 1.66 | 0.24 | [1.17, 2.15]

• 手动指定分类学校代码：即使在estimate_means中明确列出所有学校代码，School C仍无结果：

realmodel <- lm(realscore ~ SchoolCode * ChildTIS * ChildSex, data = gcsubsetdata)
means_complex_real <- estimate_means(realmodel, at = list(SchoolCode = c("A","B","C","D")))
means_complex_real

输出：

SchoolCode | Mean |   SE |       95% CI
---------------------------------------
A          | 0.94 | 0.20 | [0.53, 1.35]
B          | 0.68 | 0.29 | [0.10, 1.25]
C          |      |      |             
D          | 1.79 | 0.29 | [1.22, 2.37]

• 子集化数据：仅保留School C和D的数据重新建模，School C依旧没有结果：

gcsubset2 <- subset(gcsubsetdata, SchoolCode=="C" | SchoolCode == "D")

realmodel <- lm(realscore ~ SchoolCode * ChildTIS * ChildSex, data = gcsubset2)
means_complex_real <- estimate_means(realmodel, at = "SchoolCode")
means_complex_real

输出：

SchoolCode | Mean |   SE |       95% CI
---------------------------------------
C          |      |      |             
D          | 1.78 | 0.34 | [1.05, 2.50]

问题根源排查 #

• 数据验证：realscore、SchoolCode、ChildTIS、ChildSex均无缺失值，未调整均值计算正常；School C有10个数据点，重命名学校代码后问题依旧。
• 模型系数检查：三因素交互模型中，School C的两个交互项系数为NA，说明模型存在完全共线性，无法估计这两个系数：

Call:
lm(formula = realscore ~ SchoolCode * ChildTIS * ChildSex, data = gcsubsetdata)

Coefficients:
(Intercept)                    SchoolCodeSchB
3.26627                          -3.19701
SchoolCodeSchC                    SchoolCodeSchD
-0.54741                           2.73317
ChildTIS                          ChildSex
-0.22716                          -0.95491
SchoolCodeSchB:ChildTIS           SchoolCodeSchC:ChildTIS
0.19012                                NA
SchoolCodeSchD:ChildTIS           SchoolCodeSchB:ChildSex
-0.16173                           1.42028
SchoolCodeSchC:ChildSex           SchoolCodeSchD:ChildSex
0.77606                          -1.66379
ChildTIS:ChildSex  SchoolCodeSchB:ChildTIS:ChildSex
0.09136                          -0.07284
SchoolCodeSchC:ChildTIS:ChildSex  SchoolCodeSchD:ChildTIS:ChildSex
NA                           0.14938

有效解决方法：简化模型 #

改用仅包含主效应的线性模型后，所有学校的EMM结果正常显示：

# 基于学校、在校时长、性别的主效应模型计算EMM
realmodel <- lm(realscore ~ SchoolCode + ChildTIS + ChildSex, data = gcsubsetdata)

means_complex_real <- estimate_means(realmodel, at = "SchoolCode")

means_complex_real

输出：

Estimated Marginal Means

SchoolCode | Mean |   SE |       95% CI
---------------------------------------
A          | 0.96 | 0.19 | [0.58, 1.34]
B          | 0.69 | 0.25 | [0.19, 1.19]
C          | 1.62 | 0.29 | [1.04, 2.20]
D          | 1.79 | 0.27 | [1.24, 2.34]

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