17、多级分析：R 语言中的多级建模与预测

多级分析：R 语言中的多级建模与预测

1. 多级回归概述

在数据分析中，为解决一些问题，我们可以采用能够分离出因背景因素导致的方差的分析方法，多级回归分析（也称为混合效应回归）就是这样一种方法。这里不会深入探讨这种高度复杂分析的计算细节，仅提供理解和进行基础分析所需的信息。线性回归的诊断方法同样适用于多级回归，此外还需进行额外的诊断，如检查二级残差的正态性。

1.1 回归计算基础

在一般回归中，一个观测的标准属性值的计算方式为：
- 截距（所有纳入的预测变量值都为 0 时的平均值）
- 斜率系数乘以预测变量的值（每个预测变量）
- 残差（预测值与观测值之间的差异）

回归算法的任务是找到能使整个样本残差最小化的参数。

1.2 多级建模中的计算方式

1.2.1 随机截距和固定斜率

• 仅考虑一级预测变量且所有组斜率相同时 ：一个观测的标准属性值的计算大致为以下各项之和：
  • 公共截距
  • 组特定残差（组截距与公共截距之间的差异）
  • 斜率系数乘以预测变量的值（每个预测变量）
  • 观测特定残差
  • 在这种模型中，一级属性的效应在所有组中被认为是相同的，只有截距会变化。
• 考虑一级和二级预测变量且所有组斜率相同时 ：计算大致为以下各项之和：
  • 公共截距
  • 组特定残差（组截距与公共截距之间的差异）
  • 总体斜率系数乘以一级预测变量的值（每个一级预测变量）
  • 观测特定残差

1.2.2 随机截距和随机斜率

之前我们假设一级预测变量的斜率在所有组中是相同的，但实际情况并非总是如此。下面通过一个例子来进行说明，使用从真实数据生成的模拟数据，包含倦怠（个人成就感、去个性化和情绪耗竭）和工作满意度等属性。

以下是加载协方差数据并生成数据集的代码（每个医院有 100 个观测值，共 17 家医院）：

library(MASS)
set.seed(999)
Covariances = read.table("Covariances.dat", sep = "\t", header=T)
df = data.frame(matrix(nrow=0,ncol=4))
colnames(df) = c("Hospital","Accomp","Depers","Exhaus","WorkSat")
for (i in 1:17){
  if(i == 1) {start_ln = 1}
  else start_ln = 1+((i-1)*4)
  end_ln = start_ln + 3
  covs = Covariances[start_ln:end_ln, 3:6]
  rownames(covs)=Covariances[start_ln:end_ln,2]
  dat=mvrnorm(n=100, c(rep(0,4)), covs)
  df = rbind(df,dat)
}
df$hosp = as.factor(c(rep(1,100), rep(2,100), rep(3,100),   
                      rep(4,100), rep(5,100), rep(6,100),    
                      rep(7,100),rep(8,100),rep(9,100), 
                      rep(10,100),rep(11,100),rep(12,100),
                      rep(13,100),rep(14,100),rep(15,100),
                      rep(16,100),rep(17,100)))

以下代码将绘制每家医院去个性化与工作满意度之间的关系（使用 lm() 模型）：

library(lattice)
attach(df)
xyplot(WorkSat~Depers | hosp, panel = function(x,y) { 
  panel.xyplot(x,y) 
  panel.lmline(x,y)
})

从绘制的图中可以看出，去个性化和工作满意度之间通常呈负相关，但不同组的这种模式表现程度不同，在某些情况下甚至不存在这种关系。是否在分析中考虑这些差异由分析师决定。随机斜率模型是指允许斜率在组间变化的模型。在随机斜率模型中，个体值的计算大致为以下各项之和：
- 公共截距
- 组特定残差（组截距与公共截距之间的差异）
- 对于包含随机斜率的每个预测变量，组特定系数乘以预测变量的值。组特定系数由固定部分（公共斜率）和随机部分（每个斜率对应组围绕该斜率的变化残差）组成
- 对于斜率不允许在组间变化的预测变量（如果有），斜率系数乘以预测变量的值
- 观测特定残差

多级回归的任务仍然是找到能使残差最小化的参数。需要注意的是，一些预测变量的效应可以定义为在组间变化，而其他的则定义为不变化。

1.3 多级建模在 R 语言中的应用

1.3.1 加载数据集

为了在 R 语言中构建多级模型并预测未知数据，首先加载之前使用相同方法生成的数据集（属性未进行缩放），同样是 17 家医院，每家医院有 100 个生成的观测值：

NursesML = read.table("NursesML.dat", header = T, sep = " ")

1.3.2 零模型

我们将医院和观测值作为分析单位，来研究属性的变化情况，即比较医院层面和观测层面的变化哪个更大。可以手动进行计算，以下代码计算了每家医院我们想要预测的属性（工作满意度 WorkSat ）的均值：

means = aggregate(NursesML[,4], by=list(NursesML[,5]), 
                  FUN=mean)[2]

可以通过以下方式显示医院和观测值中工作满意度的方差：

var(unlist(means)) # 医院层面
var(NursesML[,4]) # 观测层面

输出结果显示，医院层面的方差为 0.0771365，观测层面的方差为 0.7914461。观测层面的方差远大于医院层面的方差，但两个层面的方差相互影响，因此这些结果并不可靠。使用多级建模来研究这种差异才是正确的比较方法，为此需要拟合一个只包含常数和聚类变量的多级模型，即零模型。首先安装并加载 lme4 包，该包中的 lmer() 函数可用于拟合多级模型：

install.packages("lme4"); library(lme4)
NursesML$hosp = factor(NursesML$hosp)
null = lmer(WorkSat ~ 1 + (1|hosp), data=NursesML)
summary(null)

零模型的输出结果如下：
| 效应类型 | 方差 |
| ---- | ---- |
| 随机效应 - 截距方差 | 0.06988（医院层面） |
| 残差方差 | 0.72564（观测层面） |

工作满意度的总方差为 0.06988 + 0.72564 = 0.79552。可以计算医院层面的方差比例（即组内相关系数）为 0.06988 / 0.79552 = 0.08784191，约 9% 的方差在医院层面。一般来说，如果二级（这里指医院）的方差小于 5%，可以使用传统回归进行分析，而不必过于担心数据的嵌套问题（前提是二级没有预测变量）。在固定效应中，只有截距，其值为 5.10679，这意味着观测值的平均值约为 5.10，与使用 mean() 函数计算的结果 5.106792 相同。可以通过以下代码获取每家医院的截距：

coef(null)

同时，之前已经计算了医院层面的均值（存储在 means 对象中），可以轻松显示这些值。以下是截距和均值的对比：
| 医院编号 | 截距 | 均值 |
| ---- | ---- | ---- |
| 1 | 5.313038 | 5.334455 |
| 2 | 4.945022 | 4.928223 |
| 3 | 5.810795 | 5.883899 |
| 4 | 5.113663 | 5.114376 |
| 5 | 5.184669 | 5.192756 |
| 6 | 5.055792 | 5.050496 |
| 7 | 5.359504 | 5.385746 |
| 8 | 4.975973 | 4.962389 |
| 9 | 4.759738 | 4.7237 |
| 10 | 5.16671 | 5.172932 |
| 11 | 5.330523 | 5.353755 |
| 12 | 4.829187 | 4.80036 |
| 13 | 5.042653 | 5.035993 |
| 14 | 5.061342 | 5.056623 |
| 15 | 4.999852 | 4.988747 |
| 16 | 5.056085 | 5.05082 |
| 17 | 4.810918 | 4.780195 |

这里显示的截距由固定部分（总体截距）和随机部分（每家医院一个值，对应每家医院的偏差）组成。随机部分可以使用 ranef() 函数获取：

ranef(null)

需要注意的是，由于目前分析中尚未包含预测变量， coef() 函数仅返回上述截距。在引入随机斜率后，将研究如何检验残差的正态性。

1.3.3 随机截距和固定斜率模型

现在进行第一次分析，研究个人成就感、去个性化和情绪耗竭对工作满意度的相对影响，暂时不考虑预测变量效应在医院之间的潜在变化。在运行分析之前，通常会将预测变量围绕总均值进行中心化处理。以下是加载训练和测试数据集的代码（每家医院有 50 个观测值）：

NursesMLtrain = read.table("NursesMLtrain.dat",
                           header = T, sep = " ")
NursesMLtest = read.table("NursesMLtest.dat",
                          header = T, sep = " ")

确保两个数据集中的 hosp 属性被视为因子：

NursesMLtrain$hosp = factor(NursesMLtrain$hosp)
NursesMLtest$hosp = factor(NursesMLtest$hosp)

在训练数据中拟合模型：

model = lmer(WorkSat ~ Accomp + Depers + Exhaust + (1|hosp),
             data=NursesMLtrain, REML = F)

首先要确定刚计算的模型是否比零模型更适合数据，这需要比较两个模型的 -2loglikelihood 值。由于现在使用的数据集与计算零模型时不同，需要重新拟合零模型。另外，使用受限最大似然估计（REML， lmer() 默认使用的估计器）比较 -2loglikelihood 值不可靠，因此使用最大似然估计（ML），通过设置 REML = F 来实现：

null = lmer(WorkSat ~ 1 +  (1|hosp), data=NursesMLtrain, REML = F)
anova(null, model)

模型比较的输出结果如下：
| 模型 | AIC | BIC | -2loglikelihood | deviance | Chisq | 自由度 | p 值 |
| ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- |
| null | | | | | | | |
| model | | | | | | | *** |

AIC 和 BIC 列分别指赤池信息准则和贝叶斯信息准则，与 -2loglikelihood 值一样，它们都是衡量数据与模型拟合程度的指标，但会考虑模型的复杂度（包含的参数数量）。偏差列也是衡量模型拟合度的指标。AIC、BIC 和偏差值越小越好，而 -2loglikelihood 值应随着拟合度的提高而增加。Chisq 列指模型之间 -2loglikelihood 的差异（近似服从卡方分布）。卡方检验的自由度计算为两个模型自由度的差异。最后一列是检验的 p 值，三个星号表示模型在接近 0 的值处显著，由此可以得出包含三个预测变量的模型比零模型更好的结论。

可以使用 MuMIn 包中的 r.squaredLR() 函数计算模型解释的额外方差部分，首先安装并加载该包：

install.packages("MuMIn"); library(MuMIn)
r.squaredLR(model,null)

输出结果显示，（伪）R 平方值为 0.189，这意味着我们的模型预测了零模型中约 19% 的方差。可以查看模型的摘要信息：

summary(model)

不幸的是， lmer() 函数不会提供系数的 p 值，而在假设检验中获取 p 值通常被认为是必不可少的（尽管有时更倾向于检查置信区间）。为了获取 p 值，需要自己进行计算。需要注意的是，多级建模中 p 值计算的可靠性存在一些争议，这也是输出中默认不包含 p 值的原因。首先提取 t 值，然后将其与正态分布进行比较并输出：

tvals = coef(summary(model))[,3]
tvals.p <- 2 * (1 - pnorm(abs(tvals)))
round(tvals.p,3)

输出结果显示，截距、个人成就感、去个性化和情绪耗竭在 p < .001 时显著不同于 0：
| 变量 | p 值 |
| ---- | ---- |
| (Intercept) | 0.000 |
| Accomp | 0.000 |
| Depers | 0.019 |
| Exhaust | 0.000 |

查看固定效应，可以看到在每个预测变量的均值水平（使用中心化属性）下，当所有预测变量都处于平均水平时，平均工作满意度为 5.11854。个人成就感每增加一个单位，工作满意度增加 0.17611；去个性化和情绪耗竭每增加一个单位，工作满意度分别减少 0.07335 和 0.29215。当然，这些只是估计值。可以使用 confint.merMod() 函数获取置信区间：

confint.merMod(model)

输出结果显示了 95% 置信度下的真实值，检查置信区间是否包含 0 是确定预测变量是否显著的另一种方法。可以注意到所有预测变量的置信区间都不包含 0（意味着它们是显著的）。 .sig01 和 sigma 值分别指二级（ .sig01 ）和一级（ sigma ）的标准差，两者都不为 0（因为置信区间不包含 0）：
| 变量 | 2.50% | 97.50% |
| ---- | ---- | ---- |
| .sig01 | 0.1740378 | 0.39715380 |
| sigma | 0.7222870 | 0.79515650 |
| (Intercept) | 4.9773312 | 5.25975309 |
| Accomp | 0.1013788 | 0.25092248 |
| Depers | -0.1350151 | -0.01166813 |
| Exhaust | -0.3498438 | -0.23442925 |

2. 随机截距和随机斜率模型

2.1 模型拟合与比较

在之前的模型中，我们假设所有医院的斜率是相同的。但现在，我们希望对一般医院（总体）得出结论，而不仅仅是对收集数据的医院。因此，需要允许斜率在不同医院之间变化。从数据的直观观察来看，各医院斜率存在简单变化，这也表明模型中应包含随机斜率。

我们拟合一个新的带有随机斜率的模型：

modelRS = lmer(WorkSat ~ Accomp + Depers + Exhaust + 
               (1+Accomp+Depers+Exhaust|hosp), data=NursesMLtrain, REML = F)

将这个模型与零模型进行比较：

anova(null, modelRS)

虽然输出未给出，但结果表明新模型比零模型更拟合数据。

2.2 残差正态性检验

2.2.1 二级残差检验

使用 sjPlot 包中的 sjp.lmer() 函数来检验二级残差的正态性：

install.packages("sjPlot"); library(sjPlot)
sjp.lmer(modelRS, type = "re.qq")

从绘制的图中可以看出，截距和每个预测变量的残差大致呈正态分布，但情绪耗竭存在一些偏差。

2.2.2 一级残差检验

使用 qqnorm() 函数对一级残差进行同样的操作：

qqnorm(resid(modelRS))

结果显示一级残差也大致呈正态分布。

2.3 模型方差解释与系数检验

使用以下代码观察模型解释的方差约为 19%：

r.squaredLR(model,null)

查看模型的详细信息：

summary(modelRS)

输出结果中，在随机效应部分，我们注意到出现了三个预测变量在二级的斜率方差、标准差以及它们之间的相关性；固定效应中的系数也与之前的模型不同。

再次检验预测变量对工作满意度的影响：

tvals = coef(summary(modelRS))[,3]
tvals.p <- 2 * (1 - pnorm(abs(tvals)))
round(tvals.p,3)

输出结果显示，所有三个预测变量在 p < .05 时显著：
| 变量 | p 值 |
| ---- | ---- |
| (Intercept) | 0 |
| Accomp | 0.002 |
| Depers | 0.010 |
| Exhaust | 0 |

2.4 斜率绘制

使用 languageR 包中的 plotLMER.fnc() 函数绘制斜率：

install.packages("languageR"); library(languageR)
par(mfrow=c(1,3))
plotLMER.fnc(modelRS)

由于预测变量进行了中心化处理，x 轴上的 0 值代表预测变量的均值。

2.5 抽样对估计的影响

如果想了解抽样对估计的影响，可以运行以下代码，每次运行代码时模型可能会有小的差异。注意使用不同的模型名称，以免影响后续结果：

# 加载初始数据集
NursesML = read.table("NursesML.dat", header = T,
                      sep = " ")
NursesML$hosp = factor(NursesML$hosp)
# 创建训练和测试集（各占 50%）
library(caret)
trainObs = createDataPartition(NursesML[,5], p = .5,
                               list=F)
NursesMLtrain = NursesML[trainObs,]
NursesMLtest = NursesML[-trainObs,]
# 预测变量的总均值中心化
for (i in 1:3){
  NursesMLtrain[i] = NursesMLtrain[i]-
    colMeans(NursesMLtrain[i])
  NursesMLtest[i] = NursesMLtest[i]-
    colMeans(NursesMLtest[i])
}

3. 使用多级模型进行预测

3.1 使用 predict() 函数进行预测

现在模型已经准备好，可以在测试数据集中预测工作满意度。一种方法是使用 predict() 函数， allow.new.levels 参数指定是否允许分析中出现新的医院。由于训练集和测试集中的医院相同，将其值设为 F （默认值）：

NursesMLtest$predicted = predict(modelRS, NursesMLtest,
                                 allow.new.levels = F)

3.2 多级建模分析流程总结

以下是多级建模分析的主要流程 mermaid 流程图：

graph LR
    A[加载数据集] --> B[零模型分析]
    B --> C{模型是否合适}
    C -- 否 --> B
    C -- 是 --> D[随机截距和固定斜率模型]
    D --> E[模型比较]
    E --> F{新模型是否更好}
    F -- 否 --> D
    F -- 是 --> G[随机截距和随机斜率模型]
    G --> H[残差正态性检验]
    H --> I{残差是否正态}
    I -- 否 --> G
    I -- 是 --> J[预测未知数据]

3.3 多级建模关键步骤总结

• 数据准备 ：加载数据集，确保聚类变量为因子类型，必要时对预测变量进行中心化处理。
• 零模型分析 ：拟合只包含常数和聚类变量的零模型，比较不同层面的方差，计算组内相关系数。
• 模型拟合 ：依次拟合随机截距和固定斜率模型、随机截距和随机斜率模型。
• 模型比较 ：使用 anova() 函数比较不同模型的 -2loglikelihood 值，判断模型的优劣。
• 残差检验 ：检验二级和一级残差的正态性。
• 系数检验 ：计算系数的 p 值和置信区间，判断预测变量的显著性。
• 预测 ：使用 predict() 函数对未知数据进行预测。

通过以上步骤，我们可以在 R 语言中完成多级建模分析，并对工作满意度进行有效的预测。多级建模能够更好地处理数据的嵌套结构，考虑不同层面的变异，从而提高模型的准确性和可靠性。