R 统计编程和数据建模高级教程（四）

原文：Advanced R Statistical Programming and Data Models

协议：CC BY-NC-SA 4.0

九、缺失数据

缺失数据在几乎所有现实世界的分析中都很常见。本章正式介绍缺失数据的概念，包括描述缺失的常用方法。然后，我们讨论在分析中处理缺失数据的一些潜在方法。我们将在本章中使用的主要软件包是mice软件包，该软件包提供了处理缺失数据和最小化缺失数据对分析结果的影响的强大功能[95]。

library(checkpoint)
checkpoint("2018-09-28", R.version = "3.5.1",
  project = book_directory,
  checkpointLocation = checkpoint_directory,
  scanForPackages = FALSE,
  scan.rnw.with.knitr = TRUE, use.knitr = TRUE)

library(knitr)
library(ggplot2)
library(cowplot)
library(lattice)
library(viridis)
library(VIM)

library(mice)
library(micemd)
library(parallel)

library(data.table)
library(xtable)
library(JWileymisc) # has data

options(width = 70, digits = 2)

9.1 概念背景

当一个或多个变量的观测值缺失时，就会出现数据缺失。丢失数据的潜在原因有很多。缺失数据还会导致数据分析中的问题:

• 仅从非缺失数据获得的估计值可能有偏差。

• 丢失数据会减少信息，从而降低效率。

• 许多工具和软件不具备处理缺失数据的能力。

当数据并非完全随机缺失时，就会出现偏差，因此从观察数据点获得的估计值与所有情况都观察到的估计值存在系统性差异。发生效率损失是因为丢失了一些数据，数据中的信息更少了。此外，许多处理缺失数据的简单技术，如完全病例分析(列表式删除),会丢弃任何变量缺失的病例，即使有其他变量的数据，也会导致效率损失。最后，由于许多工具无法直接处理缺失数据，因此使用缺失数据进行分析通常更加耗时和复杂，因为在使用标准工具之前，可能需要首先处理缺失数据。

有多种方法可以解决所讨论的三个缺失数据问题。但是，在实现特定的方法之前，有一些一般的注意事项。基本上，我们假设目标是估计某个参数， θ ，并获得我们估计的不确定性的估计， V AR ( θ )。我们假设数据是从我们感兴趣的人群中抽取的，尽管这不一定是整个人群。基本的问题是在什么条件和假设下, θ 的无偏估计存在，以及我们如何获得它？

数据通常被分类为

• 当缺失机制完全独立于 θ 的估计时，完全随机缺失(MCAR)

• 当缺失机制有条件地独立于 θ 的估计时，随机缺失(MAR)

• 当缺失机制与 θ 的估计相关联时，非随机缺失(MNAR)

当且仅当数据是 MCAR 时，完整的病例分析(列表式删除)将产生对 θ 的无偏估计。在所有其他情况下，我们必须考虑我们将对缺失过程做出什么样的假设，以及利用我们假设的缺失过程下的可用数据，我们是否可以获得 θ 的无偏估计。

一种表示预期模型和解释假设的因果和缺失过程的系统方法是通过使用图形模型，例如有向无环图(Dag)。Mohan、Pearl 和 Tian [66]使用 Dag 开发了识别无偏估计的方法，并将其扩展到 Mohan 和 Pearl [64]的进一步因果查询，以及测试数据是否为 MCAR 或变量水平 MAR+ [65]。

这项研究有一些实际意义。首先，根据因果过程和缺失机制， θ 的无偏估计可能可恢复，也可能不可恢复，即使是从相同的数据中。因此，考虑模型规格至关重要。简单地调整所有可能的变量是不够的，因为其他工作表明，当向 DAG [73]添加不适当的变量或路径(“后门标准”)时，先前对 θ 的无偏估计可能会变得有偏差。第二，严格解释模型和能够或不能获得无偏估计的条件是复杂的，需要比大多数研究人员和现有科学水平更详细的关于因果过程和遗漏机制的知识或期望。随着所考虑的变量数量的增加，复杂性呈指数增长，尤其是当多个变量缺失时。

虽然在 Mohan 和 Pearl 的理论工作中没有详细说明，但在实际应用时，不仅需要指定正确的变量和路径，而且需要正确指定关系的函数形式。因此，对于P(Y | f(X)，f ( X )可能是一个线性或其他一些非线性函数，而对 θ 的无偏估计将取决于正确指定这种形式。在实践中，至少应该检查观察值的关系形式，并且在必要时允许是非线性的。

虽然提供特定变量的必要背景以确定缺失的数据机制超出了本章的范围，但它们值得考虑，至少如果没有证据来指导决策，这应被视为研究中的一个限制。下一节讨论多重插补作为处理缺失数据的一种方法，并至少提供一些处理未知函数形式的选择。

多重插补

一般

存在许多插补方法。不鼓励使用单一插补(如平均值和热卡插补)和条件平均值插补，不做进一步讨论。无论采用何种插补模型，单一插补的问题在于，它们假设缺失值的预测没有错误。这导致基于估算数据的模型的不确定性估计过低，增加了 1 型错误率。影响程度将取决于缺失数据的数量以及插补模型的准确(或不准确)程度。

除了一次性输入数据集，还可以多次输入。多重插补背后的原理是，多重插补不是为插补生成单一预测值，而是从预测分布中随机抽取。预测分布的分布捕捉了预测中的不确定性，这导致了每个估算数据集中的可变性。

生成多重估算数据的一般过程如下:

• 有缺失数据的初始数据集。

• 用初始估计值填充缺失的数据值。初始估计通常很容易生成，例如每个变量的平均值或中值，或者从变量中抽取的随机值。

• 对于每个变量，建立一个模型，从剩余的变量中预测它。使用模型来预测缺失的数据。对于参数模型，如线性回归，从预测分布中提取随机值(例如，基于预测值和标准偏差的均值正态分布)或从贝叶斯模型的后验分布中提取样本。对于非参数模型，如随机森林模型[28]，不确定性可能通过其他方式引入，如 bootstrapping，以围绕预测值建立经验分布，并从中取样。

• 重复上一步，直到预测的完整数据集与上一次迭代没有实质性的不同，表明模型已经收敛。这一步是必需的，因为随着初始估计(通常相当糟糕)变得更加准确，所建立的模型可能会随着时间而改变。

• 重复所有这些步骤，以生成所需数量的多重估算数据集(例如，5，100)。

生成多个估算数据集后，分析会与无遗漏地分析单个数据集略有不同。多重估算数据分别进行分析，结果使用鲁宾规则进行汇总[58]。简单来说:

• 模型在每个估算数据集中单独进行估计。

• 对每个模型的许多参数进行平均，包括回归系数、平均差异和预测值。一些值不是通过平均直接组合的，例如残差方差或相关的估计，其通常在平均之前被转换。

• 汇总结果的不确定性估计(例如，标准误差)包括两个可变性来源:

  • 每个模型参数的平均不确定性

  • 模型之间参数估计的可变性

• 基于平均参数和不确定性估计的统计测试，从模型和不同的数据集获取不确定性。

一般来说，任何分析都应该在单个数据集上进行，只有在最后一步，结果才应该跨数据集进行汇总。例如，如果目标不仅仅是回归模型，而是生成预测值的图表，那么应该为每个估算数据集单独生成预测值，然后将预测值组合起来。不要首先合并回归模型结果，然后使用合并的模型来生成预测。关于多重插补的一个常见误解是，只应插补自变量，而不应插补因变量。模拟表明，不输入因变量会增加偏差[67]，建议输入所有变量[38]

多重插补方法

一般来说，乘以估算数据的模型采用两种形式之一:

• JM:一种联合模型(JM ),其中指定了所有变量的联合多元分布。通过从条件分布中提取来估算值。

• FCS:一个完全条件规范(FCS ),其中为每个缺失变量建立一个单独的模型，使用所有其他变量作为预测因子。每个缺失变量的值通过从该特定变量的模型的条件密度中提取来估算。

可以证明 JM 方法在理论上是有效的【85】，这是相对于 FCS 方法的一个优势，FCS 方法目前缺乏强有力的理论依据来解释为什么它们可以产生有效的结果。

JM 方法的一个挑战是，它要求多元分布(1)是特定的，(2)是数据的合理表示，以及(3)在分析上是易处理的。通常，这些条件中的一些(但不是全部)很容易满足。例如，使用多元正态(MVN)模型是相对容易指定的，并且在分析上易于估计和借鉴。然而，只有在研究案例的子集里，所有带有缺失数据的变量才能被合理地解释为正态分布。例如，一旦引入一个二元协变量(如性别)，MVN 模型就开始变得不现实。在数据不太可能真实 MVN 的情况下，MVN 模型有时仍被用作近似值，但这凸显了 JM 的一个缺点。MVN 模型不是唯一可以指定的多元分布类型。然而，涉及正态分布和二分或计数型变量的多元分布变得非常难以估计和提取，通常使它们难以分析。

我们认为，JM 方法在实际环境中特别困难，因为许多带有缺失数据的应用研究问题(如流行病学研究、纵向研究)通常具有包含许多协变量和焦点预测因子的模型，更不用说退出或缺失的潜在决定因素了，尽管这些因素不是最终分析模型的一部分，但可能是改善缺失数据模型的关键因素。如此多的不同变量都遵循任何方便的多元分布似乎很少见，这使得 JM 方法在许多情况下不切实际。

尽管缺乏强有力的理论基础来支持它作为一种合适的技术，但 FCS 方法是 JM 方法的一种令人信服的替代方法。FCS 方法也称为通过链式方程的多重插补(MICE ),因为有一系列独立的方程。在有 k 个变量的地方，不是定义一个 k 维的多元密度，而是定义 k 维的单变量密度。这既简化了计算方面，又使得使用任何可用的分布相对容易，从而适应各种各样的结果(例如，正态连续、泊松计数、二项式二元、伽玛严格正连续、结果等)。).FCS 的另一个优点是可以为预测因子和每个结果之间的关系指定非常灵活的函数形式。

进行 FCS 的一种标准方式是从每个变量的(条件)后验分布中获取 Gibbs 样本，通常形式为P(Y_I|Y_—I，θ _i )对于 1 中的 i ，… ，k 。这些可用于“填充”每个 k 结果的缺失数据。一旦所有结果都被填入，现在可以使用每个结果的预测因子的观察和填入信息重复该过程，并且该过程反复重复直到达到收敛。一旦实现收敛，就可以从条件后验分布中抽取所需数量的样本，以生成所需数量的多重估算数据集。

Van Buuren、Brand、Groothuis-Oudshoorn 和 Rubin [94]以及 Van Buuren 和 Groothuis-Oudshoorn [95]对 FCS (MICE)进行了详细描述。FCS 的一个问题是条件分布可能不兼容。当条件分布不相容时，可变插补的顺序很重要，尽管这些影响可能相对较小。Hughes 及其同事[46]已经讨论了条件分布相容的条件，从而产生了与联合分布相同的收益，他们还提供了一个模拟来检验分布相容和不相容时 FCS 的性能。总的来说，虽然理论上不如 JM 方法合理，但在模拟和实践中，FCS 似乎表现得非常好，并且，由于其限制较少的假设，当 JM 假设不成立时，它可以胜过 JM 方法。由于增加了灵活性，这一章集中于估算的 FCS 方法。

非线性效应和非正常结果

在许多应用中，默认的假设是变量之间是线性相关的。然而，这种假设应该检查，非线性函数形式是允许的，如理论或数据所规定的。即使当数据是 MAR 时，为了使估计无偏，必须正确地指定模型，包括正确的变量和正确的函数形式。

以前，所谓的被动插补用于非线性项。例如，对于两个变量之间的相互作用， X ₁ ，X ₂ ，每个变量将被单独插补，然后经过分析中的多重插补后，基于插补的 X ₁ 和 X 将形成乘积项，X₁X₂模拟表明，这种被动插补会导致有偏估计[101]。相反，应将这些非线性转换纳入多重插补，首先创建变量，然后将其作为附加变量纳入插补模型。交互作用、平方或其他转换变量有时被称为“另一个变量”(JAV)。然而，这种方法的一个限制是，X1 和)之间的关系(或交互)丢失了。Vink 和 van Buuren [100]开发了另一种方法，称为多项式组合方法，既能产生无偏估计，又能保持X1 和)之间的关系。然而，据我们所知，它目前只适用于平方项，不适用于不同变量之间的相互作用，也不适用于高阶多项式，尽管理论上这项工作可能会推广。

一个相关的问题是偏斜或非正态分布变量。一种直观的方法是尝试一些规范化转换；然而，von Hippel [102]表明，这导致了不太理想的结果，并且在许多情况下建议简单地使用正常模型。备选方案是对每个结果使用(正确的)分布(例如，伽玛、贝塔等。).目前，大多数模拟和研究仅检查了正态、二项式结果的多重插补的特性，检查有序结果的较少；因此，对于来自替代分布的变量几乎没有指导。此外，大多数软件实现允许相当有限的一组已知发行版。

正如交互作用或其他非线性项应该在需要的地方被估算一样，如果一些其他变量 Y 不仅取决于 X ₁ ， X ₂ ，而且还取决于它们的交互作用X₁X₂，为了正确地指定模型，交互作用项必须作为预测项包含在 FCS 中。如果 Y 依赖于X1，X2，X1X₂，但是X₁X₂在插补模型中被省略，结果将但是，如果 Y 依赖于 X ₁ ， X ₂ ，而是 X ₁ ， X ₂ ，X₁X_{2Bartlett、Seaman、White 和 Carpenter [3]对此进行了更详细的探讨，他们表明实质性分析模型可以是插补模型的限制性版本，但反之则不然。具体来说，他们建议确保插补模型与实体模型兼容，或者至少是半兼容的(即实体模型是插补模型的限制版本)。}

用于插补的 GLMs

FCS 方法不需要任何特定的模型，但最常用的一类模型是广义线性模型(GLMs ),如连续、近似正态分布变量的线性回归或二分变量的逻辑回归。对于感兴趣的读者，GLMs 在第 3 和 4 章节中有更深入的介绍。GLMs 的一个好处是，它们是许多分析师熟悉的分析模型，更容易理解它们如何用于多重插补。此外，大多数 glm 在现代计算机上的估计速度非常快，这使得它成为一种计算成本相对较低的模型。尽管 GLMs 很受欢迎，速度也很快，但它也有缺点。

首先，GLMs 默认所有变量和链接函数尺度上的结果之间的线性关系。尽管对任何 GLM 模型来说都是如此，但多重插补通常涉及许多关系，这一事实加剧了这一问题。例如，考虑一个相对简单的研究问题:肥胖与血压有关吗？血压本身有两个组成部分:收缩压和舒张压。此外，许多因素可能会混淆这种关联，相关的协变量也是如此，两个可能的候选因素是年龄和社会经济地位。最后，还应包括与缺失数据相关的因素。例如，如果已知被雇佣的参与者更有可能丢失数据，那么在插补模型中包括这一点是有意义的。在这个例子中，我们刚刚开始考虑与肥胖、血压和潜在缺失相关的无数因素，我们正在处理六个变量。归根结底，我们关心的只是压力、年龄和社会经济地位(三个变量)与血压(两个变量)的关系，所以我们只需要担心六个独特的两两关联。然而，对于插补模型，我们有 15 独特的成对关联。如果其中任何一个不是线性的，要么我们的插补模型将是错误的，要么我们需要识别哪个不是线性的，并包括适当的函数形式(例如，二次函数等。).许多应用研究问题将包括更多的变量，要么是因为基础研究问题涉及更多的变量，要么是因为有更多的相关协变量或更多的因素可能与遗漏有关。让分析师手动决定每个可能的成对关联的函数形式很快变得不切实际(例如，10 个变量有 45 个唯一的对；20 个变量的 190 个唯一对，等等。) .

第二，默认情况下，GLMs 不包括变量之间的任何交互。同样，随着插补模型中变量数量的增加，分析师手动指定可能的交互次数很快变得不切实际。

第三，特别是在较小的数据集或某些应用中，要包含的变量的数量相对于观察值的数量可能很大。在这些情况下，GLMs 将倾向于过度拟合数据，特别是对于二分结果，可能导致完全分离和可估计性挑战，基本问题是 GLMs 默认不包括任何变量选择方法。

插补游戏

面对未知函数形式的挑战，一个自然的解决方案是利用试图根据经验逼近未知形式的模型。一类这样的模型是广义加性模型(gam)[40，122]，它使用平滑器和罚函数来试图逼近函数形式，而不会(非常)过度拟合数据。“基本”GAM 的扩展是位置、比例和形状参数的 GAM(GAMLSS)。GAMLSS 不仅允许对位置(如正态分布中的平均值)进行建模，还允许分布的规模(方差)作为预测值以及分布形状的函数[78，88]。尽管 gam 通常发展良好(如[40，122，123])，但很少有研究检验它们在多重插补中的特性。

de Jong、van Buuren 和 Spiess [27]进行的一项模拟研究表明，在多重插补中使用 GAMLSS 可以获得良好的结果，从而允许出现偏斜和非正态分布的结果以及非线性效应。总的来说，使用 GAMLSS 的插补比使用 GLMs 提供了更好的覆盖面；然而，值得注意的是，模拟研究[27]只检查了几个变量。不收敛和其他计算挑战可能会出现在应用研究设置与许多预测。

估算的 RFs

随机森林(RFs)是一种机器学习类型，广泛用于预测建模，因为它们往往会提供出色的结果[16]。就插补而言，RFs 比之前讨论的替代方法有几个优势。

我们在插补中注意到的一个挑战是需要正确捕捉所有变量之间的函数形式，它可能不是线性的。gam 对此提供了一个可能的解决方案，但是 RFs 也可以解决这个问题，因为树允许在变量的不同点进行分割，从而允许非线性关联。

与 gam 或 GLMs 不同，RFs 还可以识别变量之间的重要交互。即使系统地筛选所有可能性是不切实际的，RFs 也能有效地识别交互。RFs 的这一特性提供了一个很大的优势，因为有十个或更多的变量，研究人员实际上不可能指定相关变量之间的所有相互作用，即使包括所有的相互作用也会导致一个复杂的模型，可能需要非常大的样本量。

最后，RFs 是有利的，因为它们本质上具有一定程度的内置变量选择。RFs 只在有助于预测结果的变量上分裂。因此，与 GLMs 或 GAMs 不同，RFs 可以更好地处理相对于数据集中的病例数而言插补模型中存在大量变量的病例。

虽然随机森林在预测和机器学习文献中非常常用，但评估其插补性能的研究相对较少。Stekhoven 和 Bühlmann [90]在较早的一篇研究 RF 插补的论文中提出了一个生物信息学领域的 RF 模型插补以及一个R包missForest。然而，Stekhoven 和 Bühlmann [90]建议仅使用 RFs 进行单一插补，如前所述，这会因缺失数据而导致低估不确定性。最近，另外两个团队独立开发了一种使用 RFs 进行多重插补的方法，并展示了这些技术的良好性能[28，86]。在这两种情况下，模型都依赖于自举来传播输入缺失值的不确定性，以生成多个估算数据集。

其他情况

除了到目前为止讨论的多种插补情况之外，还有一些更特殊的情况需要其他方法。在事件发生时间或生存结果的情况下，多重插补是复杂的，White 和 Royston [108]讨论了这些情况下协变量的插补。对于多级数据的多重插补，van Buuren [93]提供了一份最新概述。虽然不是直接多重插补，但 Goldstein、Carpenter 和 Browne [37]描述了一个多水平模型，该模型考虑了反应和协变量(包括非正态结果/协变量和交互项)中的缺失数据。然而，这种分析的复杂性在许多情况下可能令人望而却步。

9.2 R示例

为了应用所讨论的缺失数据方法，我们将使用模拟数据来密切匹配日常研究。讨论的大多数插补模型仅适用于单水平数据，因此我们将从将每日数据折叠成整个研究的平均值开始。

## load example dataset
data("aces_daily")
draw <- as.data.table(aces_daily)[order(UserID)]
davg <- na.omit(draw[, .(
  Female = na.omit(Female)[1],
  Age = na.omit(Age)[1],
  SES_1 = na.omit(SES_1)[1],
  EDU = na.omit(EDU)[1],
  STRESS = mean(STRESS, na.rm = TRUE),
  SUPPORT = mean(SUPPORT, na.rm = TRUE),
  PosAff = mean(PosAff, na.rm = TRUE),
  NegAff = mean(NegAff, na.rm = TRUE)),
  by = UserID])

接下来，我们想在数据中添加一些缺失。对于这个例子，我们将让失踪的概率取决于感知的支持和压力水平的组合。然后，对于一些变量，我们基于这些概率产生缺失。我们随机设置压力和支持为缺失。

## missing depending on support and stress
davg[, MissingProb := ifelse(
         SUPPORT < 5,
           ifelse(STRESS > 2.5, .4, .0),
           ifelse(STRESS > 2.5, 0, .4))]

set.seed(1234)
davgmiss <- copy(davg)
davgmiss[, PosAff := ifelse(rbinom(
             .N, size = 1, prob = MissingProb) == 1,
             NA, PosAff)]
davgmiss[, NegAff := ifelse(rbinom(
             .N, size = 1, prob = MissingProb) == 1,
             NA, NegAff)]
## random missingness on stress and support
davgmiss[, STRESS := ifelse(rbinom(
             .N, size = 1, prob = .1) == 1,
             NA, STRESS)]
davgmiss[, SUPPORT := ifelse(rbinom(
             .N, size = 1, prob = .1) == 1,
             NA, SUPPORT)]
davgmiss[, Age := ifelse(rbinom(
             .N, size = 1, prob = .1) == 1,
             NA, Age)]
davgmiss[, SES_1 := ifelse(rbinom(
             .N, size = 1, prob = .1) == 1,
             NA, SES_1)]
davgmiss[, Female := factor(ifelse(rbinom(
             .N, size = 1, prob = .1) == 1,
             NA, Female), levels = 0:1,
             labels = c("Male", "Female"))]
davgmiss[, EDU := factor(ifelse(rbinom(
             .N, size = 1, prob = .1) == 1,
             NA, EDU), levels = 0:1,
             labels = c("< Uni Graduate", "Uni Graduate +"))]
## drop unneeded variables to make analysis easier
davgmiss[, MissingProb := NULL]
davgmiss[, UserID := NULL]

现在开始，直观地检查数据集中缺失的数据模式是有帮助的。这可以使用图 9-1 所示的VIM包和aggr()函数来完成。左侧显示每个变量总计中缺失数据的百分比。右侧显示了不同研究变量的不同遗漏模式，以及它们的比例。至关重要的是，我们可以看到，尽管总体数据缺失率并不高(约 15%)，但只有大约三分之一的参与者实际上拥有所有变量的完整数据。这将使一个完整的案件分析有大量的信息和权力的损失。

图 9-1

通过变量和缺失模式对缺失进行可视化总结

ggr(davgmiss, prop = TRUE,
     numbers = TRUE)

VIM包还有一个功能，使用marginplot()函数帮助识别缺失是否依赖于另一个变量。图 9-2 中显示的这些图表包含了相当多的信息。中心是成对呈现数据的基本散点图。空白显示缺失数据的分布，箱线图显示每个变量的分布，根据另一个变量是否缺失进行分层。这些图表明，当情感缺失时，压力会更大，支持会更低。

图 9-2

缺失数据的二元图。中间的点表示没有丢失的数据。空白点表示缺失的数据。箱线图根据另一个变量是否存在来总结每个变量。

par(mfrow = c(2, 2))
marginplot(davgmiss[,.(STRESS, NegAff)])
marginplot(davgmiss[,.(SUPPORT, NegAff)])
marginplot(davgmiss[,.(STRESS, PosAff)])
marginplot(davgmiss[,.(SUPPORT, PosAff)])

我们可以应用统计检验来检查一个变量的缺失是否与另一个变量相关，例如，使用如下所示的 t 检验。我们没有看到压力和情感支持的显著差异。然而，一般来说，不鼓励依赖统计显著性来确定是否将一个变量纳入插补模型，因为根据影响的大小和样本大小，结果可能不具有统计显著性，但可能会对基于多重插补数据分析的模型的结果产生影响。

## does age differ by missing on negative affect?
t.test(Age ~ is.na(NegAff), data = davgmiss)$p.value

## [1] 0.9

## does age differ by missing on positive affect?
t.test(Age ~ is.na(PosAff), data = davgmiss)$p.value

## [1] 0.14

## does stress differ by missing on negative affect?
t.test(STRESS ~ is.na(NegAff), data = davgmiss)$p.value

## [1] 0.89

## does stress differ by missing on positive affect?
t.test(STRESS ~ is.na(PosAff), data = davgmiss)$p.value

## [1] 0.17

## does social support differ by missing on negative affect?
t.test(SUPPORT ~ is.na(NegAff), data = davgmiss)$p.value

## [1] 0.49

## does social support differ by missing on positive affect?
t.test(SUPPORT ~ is.na(PosAff), data = davgmiss)$p.value

## [1] 0.17

带回归的多重插补

在本例中，我们将使用回归方法对数据进行乘法估算。由于多重插补包含随机成分(如 bootstrapping)，每次进行插补时，最终结果都会随机变化。为了使结果可重复，可以设置随机种子。为了举例和加快速度，我们只生成六个估算数据集。对于实际研究，使用更多(例如，50，100)可能有助于确保即使使用不同的种子，合并模型的最终结果也将或多或少相同。

下面的代码使用来自mice包的mice()函数来估算数据集中任何缺失的数据，davgmiss。我们要求六个插补，用m = 6。默认方法指定用于连续(?? 中的数字)变量、二进制变量、名义(多进制)变量和有序分类变量的模型类型。我们从算法的十次迭代开始，为可重复性设置种子，并在插补期间关闭消息。对system.time()的调用将返回插补所用的秒数。注意mice()的结果是一个mids类对象，它包含原始数据和所有多重插补，在一个对象中。

system.time(mi.1 <- mice(
  davgmiss,
  m = 6,   maxit = 10,
  defaultMethod = c("norm", "logreg", "polyreg", "polr"),
  seed = 1234, printFlag = FALSE)
)

##    user  system elapsed
##     2.9     0.0     3.0

使用有限次数的迭代，模型可能不会收敛。收敛可以通过图来检查，类似于贝叶斯方法。每个插补用单独的颜色绘制，如果结果变得稳定，则结果收敛，并且每个单独的插补与另一个插补没有系统的不同(这可能表明收敛到单独的局部最大值)。如果对收敛性有疑问，我们可以使用mice.mids()进行进一步的迭代，而不必重新运行最初的十次迭代。结果似乎相当稳定，没有明确的迹象表明系统差异提供了合理的证据收敛。

图 9-3

会聚的 Mice 诊断

## plot convergence diagnostics
plot(mi.1, PosAff + NegAff + SUPPORT ~ .it | .ms)

## run an additional iterations
system.time(mi.1 <- mice.mids(
  mi.1, maxit = 10,
  printFlag = FALSE)
)

##    user  system elapsed
##       3       0       3

## plot convergence diagnostics
plot(mi.1, PosAff + NegAff + SUPPORT ~ .it | .ms)

除了检查模型是否收敛之外，评估估算值是否合理也很有帮助。这可以通过使用mice包中的densityplot()绘制观察值和估算值的分布图来完成。除了单变量评估，还可以使用mice包中的xyplot()函数测试双变量关系。蓝色用于原始观察数据，红色用于估算数据。一个小问题是，使用基于回归的方法，插补是在数据范围之外产生的，特别是预测影响值低于 1，这超出了范围。虽然这可能看起来不理想，但建议保留这些值，因为排除它们或强制它们为 1 会减少预测中的可变性，并会使模型看起来比它应有的更确定。

图 9-5

观察数据和估算数据的单变量密度图，通过插补分开

图 9-4

在更多次迭代后对收敛进行 Mice 诊断

densityplot(mi.1, ˜ PosAff + NegAff + SUPPORT + STRESS)

xyplot(mi.1, NegAff + PosAff ˜ STRESS + SUPPORT)

如果在这个阶段，我们确信插补模型中没有发生任何问题，我们可以继续拟合我们的主要分析。mice包包含了with()函数的方法，当与从mice()返回的mids类对象一起使用时，该函数将给定的R表达式应用于每个估算数据集。首先，我们将运行一个线性回归模型，以积极情感作为结果，压力作为主要解释因素，其他因素作为协变量。结果是一个mira类对象，它是一个包含在单独的多重估算数据集上重复的相同分析的对象类。如果我们打印这个对象，我们会得到每个回归模型的结果。

图 9-6

带有估算数据的二元散点图，由观察数据和估算数据分别着色

lm.1 <- with(mi.1, lm(PosAff ~ STRESS + Age + EDU + Female))

lm.1

## call :
## with.mids(data = mi.1, expr = lm(PosAff ~ STRESS + Age + EDU +
##     Female))
##
## call1 :
## mice.mids(obj = mi.1, maxit = 10, printFlag = FALSE)
##
## nmis :
##  Female     Age   SES_1     EDU  STRESS SUPPORT  PosAff  NegAff
##      22      24      23      15      23      21      32      45
##
## analyses :
## [[1]]
##
## Call:
## lm(formula = PosAff ~ STRESS + Age + EDU + Female)
##
## Coefficients:
##       (Intercept)             STRESS                Age
##            4.1174            -0.2194            -0.0459
## EDUUni Graduate +       FemaleFemale
##           -0.0715             0.0450
##
##
## [[2]]
##
## Call:
## lm(formula = PosAff ~ STRESS + Age + EDU + Female)
##
## Coefficients:
##       (Intercept)             STRESS                Age
##            4.0037            -0.1697            -0.0415
## EDUUni Graduate +       FemaleFemale
##           -0.0967            -0.0464 

##
##
## [[3]]
##
## Call:
## lm(formula = PosAff ~ STRESS + Age + EDU + Female)
##
## Coefficients:
##       (Intercept)             STRESS                Age
##            3.8330            -0.2022            -0.0266
## EDUUni Graduate +       FemaleFemale
##           -0.0267            -0.2389
##
##
## [[4]]
##
## Call:
## lm(formula = PosAff ~ STRESS + Age + EDU + Female)
##
## Coefficients:
##       (Intercept)             STRESS                Age
##            4.2129            -0.2214            -0.0466
## EDUUni Graduate +       FemaleFemale
##           -0.0529            -0.0799
##
##
## [[5]]
##
## Call:
## lm(formula = PosAff ~ STRESS + Age + EDU + Female)
##
## Coefficients:
##       (Intercept)             STRESS                Age 

##            3.6608            -0.1867            -0.0221
## EDUUni Graduate +       FemaleFemale
##           -0.1658            -0.0691
##
##
## [[6]]
##
## Call:
## lm(formula = PosAff ~ STRESS + Age + EDU + Female)
##
## Coefficients:
##       (Intercept)             STRESS                Age
##            3.6097            -0.1872            -0.0232
## EDUUni Graduate +       FemaleFemale
##           -0.0190            -0.0252

我们可以检查单个模型进行模型诊断，如图 9-7 所示。

图 9-7

来自第一个估算数据集的线性回归模型诊断

par(mfcol = c(2,2 ))
plot(lm.1$analyses[[1]]) 

par(mfcol = c(1,1))

然而，一般来说，我们对单个模型不太感兴趣，相反，我们希望看到跨模型的总体结果。mice包中的pool()函数就是用来做这件事的。与summary()结合，我们可以得到通常的线性回归汇总表，但基于来自多重估算数据的汇总结果。结果是汇集的回归系数、它们的标准误差、t 值、估计的自由度、p 值和置信区间。为了很好地为这本书格式化，我们使用了xtable()函数，结果在表 9-1 中。为此，不使用LaTeX，只需运行summary(pool(lm.1), conf.int = TRUE)。

表 9-1

多个估算数据的回归结果汇总

|   |

估计

|

Std。错误

|

统计的

|

df

|

p 值

|

2.5%

|

97.5%

|
| — | — | — | — | — | — | — | — |
| (截取) | Three point nine one | Zero point six eight | Five point seven one | Eighty-nine point eight nine | Zero | Two point five five | Five point two seven |
| 强调 | −0.20 | Zero point zero four | −4.48 | Fifty point five three | Zero | −0.29 | −0.11 |
| 年龄 | −0.03 | Zero point zero three | −1.07 | Eighty-eight point eight eight | Zero point two nine | −0.10 | Zero point zero three |
| 我的优势是毕业+ | −0.07 | Zero point one five | −0.50 | Eighty-three point seven eight | Zero point six two | −0.36 | Zero point two two |
| 女性女性 | −0.07 | Zero point one five | −0.45 | Nineteen point nine eight | Zero point six five | −0.39 | Zero point two five |

xtable(summary(pool(lm.1), conf.int=TRUE),
  digits = 2,
  caption = "Regression results pooled across multiply imputed data",
  label = "tmd-pooledres1")

为了从模型中得到汇集的 R ² ，我们可以使用pool.r.squared()函数。

pool.r.squared(lm.1)

##      est lo 95 hi 95 fmi
## R² 0.14 0.054  0.26 NaN

一些额外的列通过指定参数type = "all"提供关于缺失量和缺失信息部分的信息，这是缺失数据对特定系数影响程度的指标。这些结果见表 9-2 。

表 9-2

包含附加信息的多重估算数据的回归结果汇总

|   |

估计

|

Std。错误

|

统计的

|

df

|

p 值

|

2.5%

|

97.5%

|

里弗

|

希腊字母的第 11 个

|

fonds mon é taire international 国际货币基金组织

|

取消栏

|

b

|
| — | — | — | — | — | — | — | — | — | — | — | — | — |
| (截取) | Three point nine one | Zero point six eight | Five point seven one | Eighty-nine point eight nine | Zero | Two point five five | Five point two seven | Zero point one eight | Zero point one five | Zero point one seven | Zero point four | Zero point zero six |
| 强调 | −0.20 | Zero point zero four | −4.48 | Fifty point five three | Zero | −0.29 | −0.11 | Zero point three three | Zero point two five | Zero point two eight | Zero | Zero |
| 年龄 | −0.03 | Zero point zero three | −1.07 | Eighty-eight point eight eight | Zero point two nine | −0.10 | Zero point zero three | Zero point one eight | Zero point one five | Zero point one seven | Zero | Zero |
| 我的优势是毕业+ | −0.07 | Zero point one five | −0.50 | Eighty-three point seven eight | Zero point six two | −0.36 | Zero point two two | Zero point one nine | Zero point one six | Zero point one eight | Zero point zero two | Zero |
| 女性女性 | −0.07 | Zero point one five | −0.45 | Nineteen point nine eight | Zero point six five | −0.39 | Zero point two five | Zero point eight one | Zero point four five | Zero point four nine | Zero point zero one | Zero point zero one |

xtable(summary(pool(lm.1), type = "all", conf.int=TRUE),
  digits = 2,
  caption = "Regression results pooled across multiply imputed data with additional information",
  label = "tmd-pooledres1alt")

最后，假设我们想要得到压力和积极情感之间关系的预测回归线。我们首先建立一个数据集，然后从每个独立的回归模型中生成预测，然后对结果进行平均。

图 9-8

压力和积极情感之间关联的线性回归模型的汇总预测

newdat <- data.frame(
  STRESS = seq(from = 0, to = 6, length.out = 100),
  Age = mean(davg$Age),
  EDU = factor("< Uni Graduate", levels = levels(davgmiss$EDU)),
  Female = factor("Female", levels = levels(davgmiss$Female)))

newdat$PosAff <- rowMeans(sapply(1:6, function(i) {
  predict(lm.1$analyses[[i]], newdata = newdat) 

}))

ggplot(newdat, aes(STRESS, PosAff)) +
  geom_line()

并行处理的多重插补

虽然基于回归的插补往往相当快，但随着更多的案例、更多的变量、更多的插补，特别是使用更复杂的插补方法(如随机森林),多重插补可能需要大量的计算和时间。幸运的是，如果可以的话，并行化可以直接利用多个内核。

此示例显示了如何使用一种在 Windows、Linux 和 Mac 操作系统上都适用的方法来并行化多重插补。在 Linux 和 Mac 上有更简单的策略，但是这种方法是最具交叉兼容性的。

首先，我们创建一个包含两个工作进程的本地集群。如果您有更多的内核，您可以将此设置得更高。如果只有两个内核，则设置为两个，以此类推。接下来，我们确保将包加载到每个工作进程中。由于我们将让每个工人只进行一次插补，为了使结果可复制，我们需要为每次插补单独随机种子。最后，我们需要将数据集和随机种子导出到每个工作者进行处理。

cl <- makeCluster(2)
clusterExport(cl, c("book_directory", "checkpoint_directory" ))

clusterEvalQ(cl, {
  library(checkpoint)
  checkpoint("2018-09-28", R.version = "3.5.1",
  project = book_directory,
  checkpointLocation = checkpoint_directory,
  scanForPackages = FALSE,
  scan.rnw.with.knitr = TRUE, use.knitr = TRUE)
  library(mice)
  library(randomForest)
  library(data.table)
})

## [[1]]
##  [1] "data.table"    "randomForest"  "mice"          "lattice"
##  [5] "checkpoint"    "RevoUtils"     "stats"         "graphics"
##  [9] "grDevices"     "utils"         "datasets"      "RevoUtilsMath"
## [13] "methods"       "base"
##
## [[2]]
##  [1] "data.table"    "randomForest"  "mice"          "lattice"
##  [5] "checkpoint"    "RevoUtils"     "stats"         "graphics"
##  [9] "grDevices"     "utils"         "datasets"      

"RevoUtilsMath"
## [13] "methods"       "base"

imputation_seeds <- c(
  403L, 2L, 2118700268L, 1567504751L,
  -161759579L, -1822093220L)

clusterExport(cl, c("davgmiss", "imputation_seeds"))

现在一切都设置好了，我们可以从 1 到 6 循环，得到我们的 6 个多重估算数据集。每一个都被传递给一个实际运行的工作进程，使用parLapplyLB()函数。

system.time(mi.par <- parLapplyLB(cl, 1:6, function(i) {
mice(
  davgmiss,
  m = 1,   maxit = 20,
  defaultMethod = c("norm", "logreg", "polyreg", "polr"),
  seed = imputation_seeds[i])
}))

##    user  system elapsed
##     0.0     0.0     3.5

最后，之前因为我们要求mice()直接做多重插补，所以插补已经合并成一个大对象了。为了分开并行处理，我们每次只要求一个插补，所以我们需要手动合并它们，我们用ibind()函数创建一个具有六个插补的mids类对象。当我们打印结果时，我们可以看到多重插补的数量是 6。

## combine the separate imputations into a single object
mi.par2 <- ibind(mi.par[[1]], mi.par[[2]]) 

for (i in 3:6) {
  mi.par2 <- ibind(mi.par2, mi.par[[i]])
}

mi.par2

## Class: mids
## Number of multiple imputations:  6
## Imputation methods:
##   Female      Age    SES_1      EDU   STRESS  SUPPORT   PosAff
## "logreg"   "norm"   "norm" "logreg"  "norm"  "norm"   "norm"
##   NegAff
##   "norm"
## PredictorMatrix:
##         Female Age SES_1 EDU STRESS SUPPORT PosAff NegAff
## Female       0   1     1   1      1       1      1      1
## Age          1   0     1   1      1       1      1      1
## SES_1        1   1     0   1      1       1      1      1
## EDU          1   1     1   0      1       1      1      1
## STRESS       1   1     1   1      0       1      1      1

## SUPPORT      1   1     1   1      1       0      1      1

使用随机森林的多重插补

使用随机森林输入数据的工作原理与使用更简单的回归方法基本相同。同样，基本函数调用是对mice()的调用。主要区别在于指定了不同的方法“rf”。由于 RFs 可以处理连续变量和分类变量，我们不需要为不同类型的变量指定不同的方法。最后，对于 RFs，我们可以设置一些额外的选项。首先是每片森林中应该有多少棵树。在本例中，我们将其设置为 100。一些估算表明可能需要更少的树，但这将取决于问题，在某些情况下，为了获得良好的预测模型，可能需要 100 多棵树。我们还指定了 RF 的节点大小，在本例中为 10。迭代次数设置为 20，以匹配回归插补的总迭代次数。实际上，并不总是清楚应该使用多少次迭代。模型收敛需要足够的迭代。

首先要注意的是，回归方法只需要几秒钟，而 RF 插补则需要很长时间，甚至是并行运行。该模型似乎已经收敛。最后，在检查诊断时，使用 RF 插补，分布与观察数据更加匹配，并且影响值不会被插补到可能性范围之外(低于 1)。

system.time(mi.rfpar <- parLapplyLB(cl, 1:6, function(i) 

{
  mice(
    davgmiss,
    m = 1, maxit = 30,
    method = "rf",
    seed = imputation_seeds[i],
    ntree = 500, nodesize = 10)
}))

##    user  system elapsed
##    0.16    0.11  850.27

## combine into a single object
mi.rf <- ibind(mi.rfpar[[1]], mi.rfpar[[2]])
for (i in 3:6) {
  mi.rf <- ibind(mi.rf, mi.rfpar[[i]]) 

}

## plot convergence diagnostics
plot(mi.rf, PosAff + NegAff + SUPPORT ~ .it | .ms)

## model diagnostics
densityplot(mi.rf, ~ PosAff + NegAff + SUPPORT + STRESS)

xyplot(mi.rf, NegAff + PosAff ~ STRESS + SUPPORT)

请注意，无论数据是如何多重估算的(GLMs、GAMs、RFs)，一旦它们被估算，我们将以同样的方式对它们进行分析。我们可以比较不同模型的结果来看影响。因为我们生成了缺失，所以在这个例子中，我们不仅可以将插补模型与完整的病例结果进行比较，还可以将它们与“真相”进行比较以下代码运行了一个结果为积极影响的 GLM，模型中包含一些社会人口统计协变量，压力作为焦点预测因子。使用真实数据、仅完整案例、基于链式方程和线性模型的估算数据以及随机森林多重估算来重复该模型。估计值和置信区间绘制在图 9-12 中。在这个虚构的例子中，我们可以看到完全案例方法产生了更大的置信区间和通常不太准确的结果。不同的插补方法具有比“真实”模型更具可比性的置信区间。不幸的是，虽然可以在已知条件下评估不同插补模型的性能，但在实践中，不可能知道真相或确切的缺失数据机制，因此很难判断哪种特定方法最准确。这些结果还表明，根据对发现的解释，许多结果可能被认为是非常相似的。鉴于这种相似性和哪种模型最好的不确定性，考虑与不同方法相关的计算成本是合理的，这可能导致一些人选择比随机森林更简单的插补模型。

图 9-12

压力和积极情感之间关联的线性回归模型的汇总预测

图 9-11

观察值和估算值的情感与压力和社会支持散点图

图 9-10

随机森林模型中观测值和估算值的密度图

图 9-9

随机森林插补模型的收敛诊断

m.true <- lm(PosAff ~ STRESS + Age + EDU + Female, data = davg)
m.cc <- lm(PosAff ~ STRESS + Age + EDU + Female, data = davgmiss)
m.mireg <- summary(pool(with(mi.1,
  lm(PosAff ~ STRESS + Age + EDU + Female))),
  conf.int = TRUE)
m.mirf <- summary(pool(with(mi.rf,
  lm(PosAff ~ STRESS + Age + EDU + Female))),
  conf.int = TRUE)

res.true <- as.data.table(cbind(coef(m.true), confint(m.true)))
res.cc <- as.data.table(cbind(coef(m.cc), confint(m.cc)))
res.mireg <- as.data.table(m.mireg[, c("estimate", "2.5 %", "97.5 %")])
res.mirf <- as.data.table(m.mirf[, c("estimate", "2.5 %", "97.5 %")])
setnames(res.true, c("B", "LL", "UL"))
setnames(res.cc, c("B", "LL", "UL"))
setnames(res.mireg, c("B", "LL", "UL"))
setnames(res.mirf, c("B", "LL", "UL"))

res.compare <- rbind(
  cbind(Type = "Truth", Param = names(coef(m.true)), res.true),
  cbind(Type = "CC", Param = names(coef(m.true)), res.cc),
  cbind(Type = "MI Reg", Param = names(coef(m.true)), res.mireg),
  cbind(Type = "MI RF", Param = names(coef(m.true)), res.mirf))

ggplot(res.compare, aes(factor(""),
   y = B, ymin = LL, ymax = UL, colour = Type)) +
  geom_pointrange(position = position_dodge(.4)) +
  scale_color_viridis(discrete = TRUE) +
  facet_wrap(~Param, scales = "free") +
  theme(
    legend.position = c(1, 0),
    legend.justification = c("right", "bottom"))

## clean up cluster
stopCluster(cl)
rm(cl)

9.3 案例研究:RFs 的多重插补

为了结束这一章，我们将检查一个使用随机森林多重插补的完整工作示例。随机森林多重插补往往非常耗时，因此在几乎所有情况下，您都应该考虑使用并行处理。下面的代码设置了一个本地集群，并加载了所有必需的包。正如我们在介绍中提到的，我们使用checkpoint包来确保我们可以通过指定R的版本和日期来精确地控制和指定使用哪个版本的R包。这有助于使结果具有可重复性，并确保如果您在以后回到您的代码，您确切地知道您使用了什么软件来运行它。

cl <- makeCluster(2)
clusterExport(cl, c("book_directory", "checkpoint_directory" ))

clusterEvalQ(cl, {
  library(checkpoint)
  checkpoint("2018-09-28", R.version = "3.5.1",
  project = book_directory,
  checkpointLocation = checkpoint_directory,
  scanForPackages = FALSE,
  scan.rnw.with.knitr = TRUE, use.knitr = TRUE)
  library(mice)
  library(randomForest)
  library(data.table)
})

## [[1]]
##  [1] "data.table"    "randomForest"  "mice"          "lattice"
##  [5] "checkpoint"    "RevoUtils"     "stats"         "graphics"
##  [9] "grDevices"     "utils"         "datasets"      "RevoUtilsMath"
## [13] "methods"       "base"
##
## [[2]]
##  [1] "data.table"    "randomForest"  "mice"          "lattice"
##  [5] "checkpoint"    "RevoUtils"     "stats"         "graphics"
##  [9] "grDevices"     "utils"         "datasets"      "RevoUtilsMath"
## [13] "methods"       "base"

因为多重插补包含随机成分，所以如果您想在重新运行多重插补模型时获得相同的结果，也需要设置随机种子。获得许多随机种子的一个简单方法是使用R中内置的.Random.seed变量。因为这可能会随着时间的推移而改变，而不是直接依赖于结果，所以使用dput()函数将它们导出为可复制可粘贴的R代码。在下面的代码中，我们展示了一个简单的例子，然后展示了重用种子，以便我们可以确定使用哪些种子。为了在并行处理中使用，我们必须使用clusterExport()将结果导出到本地集群。

## example of how to have R return some seed values
dput(.Random.seed[1:5])

## c(403L, 148L, -1767993668L, 1417792552L, 298386660L)

## random seeds
imputation_seeds <- c(403L, 148L, -1767993668L,
  1417792552L, 298386660L, 1360311820L,
1356573822L, -1472988872L, 1215046494L, 759520201L,
1399305648L, -455288776L, 969619279L, 518793662L,
-383967014L, -1983801345L, -698559309L, 1957301883L,
-1457959076L, 1321574932L, -537238757L,
11573466L, 1466816383L, -2113923363L, 1663041018L)

clusterExport(cl, c("davgmiss", "imputation_seeds"))

实际运行随机森林多重插补的代码相当简单。随机森林模型中的预测因子几乎不需要设置，因为模型容易适应分类和连续预测因子，并且不强加关于关联函数形式的假设。因此，对转换数据的关注不是特别相关。尽管异常值可能会有一些影响，因为随机森林依赖于数据的分割，但异常值或预测值中的极值也往往影响较小。连续结果的异常值可能会带来一些挑战，因此检查这些异常值是一个好主意。以下代码和图 9-13 显示了数据快速诊断的示例。图表显示数据可能足够好，可以继续进行。

图 9-13

插补模型中包含的连续变量的密度图

ggplot(melt(davgmiss[, sapply(davgmiss, is.numeric),
              with = FALSE], measure.vars = 1:6), aes(value)) +
         geom_density() + geom_rug() +
         facet_wrap(~variable, scales = "free")

## Warning: Removed 168 rows containing non-finite values (stat_density) 

.

首先，运行相对简单的模型并生成少量多重估算数据集通常是个好主意。总之，这意味着获得初步结果所需的时间更少，这将有助于识别设置中的任何错误，确保插补模型的结果看起来合理，并可以形成一个基础来估计完成全部插补需要多长时间。下面的代码就是这样一个简单的例子。因为我们设置了一个只有两个处理器的本地集群，所以我们创建了四个估算数据集，所以每个处理器只需要生成两个数据集。我们还将迭代次数限制为五次，这大大提高了速度。我们在开始和结束时都使用了proc.time()函数，这样我们就可以获得总共运行了多长时间的日志。

start.time <- proc.time()

mi.rfpar1 <- parLapplyLB(cl, 1:4, function(i) {
  mice(
    davgmiss,
    m = 1, maxit = 5,
    method = "rf",
    seed = imputation_seeds[i],
    ntree = 100, nodesize = 10)
})
stop.time <- proc.time()

## estimate of how long it took
stop.time - start.time

##    user  system elapsed
##       0       0      15

## combine into a single object
mi.rf1 <- ibind(mi.rfpar1[[1]], mi.rfpar1[[2]])
for (i in 3:4) {
  mi.rf1 <- ibind(mi.rf1, mi.rfpar1[[i]])
}

我们可以看到大约过了 15.5 秒。随着每个核心插补的增加，以及迭代次数的增加，这一数字将大致呈线性增加。我们还应该检查来自这个简单模型的诊断，看看是否有任何异常或可能暗示模型中的问题。我们首先来看看图 9-14 中的一些收敛诊断。

图 9-14

随机森林插补模型的收敛诊断

## plot convergence diagnostics

plot(mi.rf1, NegAff + STRESS + Age ˜ .it | .ms)

接下来，我们将插补与图 9-15 中观察到的分布进行比较。

图 9-15

随机森林模型中观测值和估算值的密度图

## model diagnostics for continuous study variables
densityplot(mi.rf1, ˜ NegAff + STRESS + Age)

现在我们可以拟合我们的目标模型，一个线性回归，从每个估算数据集中提取残差，并检查这些假设、异常值等。因为有多个回归模型(每个估算数据集一个)，我们使用lapply()循环遍历这些模型，组合成标准化残差的单个向量，并绘制。结果如图 9-16 所示。虽然分布有些偏斜，并且有一些相对极端的值，但它们并不太可怕，残差分布也没有严重偏斜。

图 9-16

模型残差的分布图(密度和 Q-Q 偏差)

## fit the models
fit.mirf1 <- with(mi.rf1,
  lm(NegAff ~ STRESS + Age + EDU + Female + SES_1))

testdistr(unlist(lapply(fit.mirf1$analyses, rstandard)))

最后，我们可以汇总和总结结果，然后在表 9-3 中查看它们。这有助于早期识别数据、编码或最终分析模型中的任何明显问题。如果有问题，可以通过这种方式更快地解决，而不是运行整个插补模型，这可能很耗时，而且只有到那时才意识到有问题。有趣的是，我们的经验是，通常在最初几次会发现一些数据问题，直到后来才发现，因为插补依赖于正确的数据，所以数据中的任何问题都需要重新运行插补模型和所有后续分析。

表 9-3

多次估算数据测试运行中汇集的回归结果

|   |

估计

|

Std。错误

|

统计的

|

df

|

p 值

|

2.5%

|

97.5%

|
| — | — | — | — | — | — | — | — |
| (截取) | One point zero five | Zero point three two | Three point three | Seventy-five point five eight | Zero | Zero point four two | One point six nine |
| 强调 | Zero point one nine | Zero point zero two | Ten point four one | One hundred and eighteen point nine seven | Zero | Zero point one five | Zero point two two |
| 年龄 | Zero point zero one | Zero point zero one | Zero point six | One hundred and three point eight eight | Zero point five five | −0.02 | Zero point zero three |
| 我的优势是毕业+ | Zero point zero two | Zero point zero seven | Zero point three four | Twenty-five point three one | Zero point seven three | −0.12 | Zero point one six |
| 女性女性 | −0.01 | Zero point zero six | −0.19 | Twenty-four point eight six | Zero point eight five | −0.14 | Zero point one one |
| SES_1 | −0.02 | Zero point zero two | −0.97 | Sixty point three nine | Zero point three three | −0.06 | Zero point zero two |

## pool results and summarize
m.mirf1 <- summary(pool(fit.mirf1), conf.int = TRUE)

xtable(m.mirf1,
  digits = 2,
  caption = "Regression results pooled across multiply imputed data test run",
  label = "tmd-pooledres2")

一旦我们合理地认为插补模型有效，没有需要解决的数据问题，并且我们的最终分析模型可能有效，那么我们就进入最终插补。在本例中，我们将最大迭代次数从 5 次增加到 30 次，以帮助确保收敛。我们还将估算数据集的数量从 4 个增加到 10 个。在实践中，更常见的是使用 25-100 个估算数据集，但我们保持它的简洁性，以使示例不会运行太长时间。

之前，我们看到每个内核进行 2 次插补需要 5 次迭代，耗时 15.5 秒。当使用 30 次而不是 5 次最大迭代时，我们预计需要大约 6 倍的时间，如果我们进行 10 次插补(每个内核 5 次)，插补次数大约需要 2.5 倍的时间。总之，我们估计完成这个更长的插补大约需要 232.5 秒。如果我们计划 50 个估算数据集，它将会是现在的 5 倍。

start.time2 <- proc.time()
mi.rfpar2 <- parLapplyLB(cl, 1:10, function(i) {
  mice(
    davgmiss,
    m = 1, maxit = 30,
    method = "rf",
    seed = imputation_seeds[i],
    ntree = 100, nodesize = 10)
})
stop.time2 <- proc.time()

## time taken
stop.time2 - start.time2

##    user  system elapsed
##    0.04    0.02  274.58

## combine into a single object
mi.rf2 <- ibind(mi.rfpar2[[1]], mi.rfpar2[[2]])
for (i in 3:10) {
  mi.rf2 <- ibind(mi.rf2, mi.rfpar2[[i]])
}

我们可以看到，与我们预测的 232.5 秒相比，大约过去了 274.6 秒。与更简单的模型一样，我们应该检查诊断。这些如图 9-17 所示。

图 9-17

随机森林插补模型的收敛诊断

## plot convergence diagnostics
plot(mi.rf2, NegAff + STRESS + Age ˜ .it | .ms)

接下来，我们将插补与图 9-18 中观察到的分布进行比较。

图 9-18

随机森林模型中观测值和估算值的密度图

## model diagnostics for continuous study variables
densityplot(mi.rf2, ˜ NegAff + STRESS + Age)

现在我们拟合我们的目标模型，线性回归，并检查图 9-19 中的标准化残差。

图 9-19

模型残差的分布图(密度和 Q-Q 偏差)

## fit the models
fit.mirf2 <- with(mi.rf2,
  lm(NegAff ~ STRESS + Age + EDU + Female + SES_1))

testdistr(unlist(lapply(fit.mirf2$analyses, rstandard)))

最后，我们可以汇总和总结结果，然后在表 9-4 中查看它们。

## pool results and summarize
m.mirf2 <- summary(pool(fit.mirf2), conf.int = TRUE)

xtable(m.mirf2,
  digits = 2,
  caption = "Regression results pooled across multiply imputed data final run",
  label = "tmd-pooledres3")

如果你将表 9-3 与表 9-4 进行比较，你会发现有一些不同。如果我们进行更多的估算，可能会有更多的差异。选择 50 至 100 个插补的一个原因是，数字越大，一组随机插补与另一组插补之间的差异往往越小。由于只有 5 到 10 次插补，随机概率可能会导致相对较大的变化。

表 9-4

多次估算数据最终运行的回归结果汇总

|   |

估计

|

Std。错误

|

统计的

|

df

|

p 值

|

2.5%

|

97.5%

|
| — | — | — | — | — | — | — | — |
| (截取) | One point one one | Zero point three three | Three point three nine | Ninety-five point five two | Zero | Zero point four six | One point seven six |
| 强调 | Zero point one eight | Zero point zero two | Eight point nine four | Sixty-nine point seven four | Zero | Zero point one four | Zero point two two |
| 年龄 | Zero point zero one | Zero point zero two | Zero point four six | Forty-four point zero three | Zero point six five | −0.02 | Zero point zero four |
| 我的优势是毕业+ | Zero point zero three | Zero point zero eight | Zero point four four | Thirty-one point one | Zero point six six | −0.12 | Zero point one nine |
| 女性女性 | −0.01 | Zero point zero seven | −0.12 | Twenty-seven point eight six | Zero point nine one | −0.15 | Zero point one four |
| SES_1 | −0.02 | Zero point zero two | −1.03 | Sixty-two point six five | Zero point three one | −0.07 | Zero point zero two |

9.4 总结

本章介绍了通过链式方程(MICE)的多重插补，这是一种解决缺失数据的灵活技术。与更容易但往往不是最佳的完整病例分析相比，MICE 可以通过使用所有可用数据来提高效率。如果与缺失相关的必要变量可以包括在 MICE 模型中，MICE 还可以减少由于缺失数据而导致的估计偏差。在许多分析中，解决缺失数据通常是关键的第一步。此外，尽管我们在本章中仅展示了一些后续分析，但 MICE 可以用作几乎任何后续分析的第一步，包括本书其他章节中讨论的所有分析。最后，表 9-5 显示了本章中用于解决缺失数据的一些关键函数的简要总结。

表 9-5

本章中描述的关键功能列表及其功能摘要

|

功能

|

它的作用

|
| — | — |
| aggr() | 通过变量可视化数据集中的缺失数据，并检查缺失数据的不同模式 |
| marginplot() | 可视化一个变量的缺失是否依赖于另一个变量 |
| mice() | 使用完全条件规格方法运行多重插补 |
| mice.mids() | 对先前运行的mice中的mids类对象运行额外的迭代，以检查收敛性 |
| densityplot() | 一个mids类对象绘制观察和估算数据密度的方法；有助于了解估算数据的分布是否与观测数据相似或不同 |
| xyplot() | 一个mids类对象创建两个变量散点图的方法，分别显示观察值和估算值 |
| with() | 一个mids类对象对每个估算数据集分别运行指定分析的方法 |
| pool() | 汇集对不同估算数据集分别运行的分析结果 |
| ibind() | 用所有多重插补将单独的mids对象合并成一个对象 |
| complete() | 从一个mids类对象中提取完整的数据集 |

十、GLMMs：简介

广义线性混合模型(glmm)扩展了前面章节中介绍的广义线性模型(GLMs ),以统计方式说明聚集的数据(例如，学校中的儿童、特定医院诊所中的个人、对同一个人的重复测量)并使这些非独立观察值有条件地独立。

考虑观察值之间的关系是至关重要的，无论是由于阶层(如班级、学校、工作场所)内的聚类，还是由于家庭内的相似性或重复测量，因为许多统计模型假设观察值是独立的，或者至少是有条件独立的。在许多情况下，这可能是一个站得住脚的假设，特别是当观察的基本单位(例如，个人等。)每个人都对我们的数据贡献了一个观察值。本章介绍了观测值可能不独立的情况。当处理观察值可能不独立的数据时，不仅需要调整最终模型和统计推断，还必须调整数据可视化、探索和描述性统计。因此，本章涵盖了处理非独立数据的各个方面，后续章节涵盖了此类数据的统计建模。

library(checkpoint)
checkpoint("2018-09-28", R.version = "3.5.1",
  project = book_directory,
  checkpointLocation = checkpoint_directory,
  scanForPackages = FALSE,
  scan.rnw.with.knitr = TRUE, use.knitr = TRUE)

library(knitr)
library(ggplot2)
library(cowplot)
library(viridis)
library(JWileymisc)
library(data.table)
library(lme4)
library(lmerTest)
library(chron)
library(zoo)
library(pander)
library(texreg)

options(width = 70, digits = 2)

10.1 多级数据

观察值不独立的两种常见情况是对每个个体进行重复测量，如纵向研究或重复测量实验，以及对个体进行聚类或分组，如家庭、学校、公司等。，然后将个体参与者聚集或嵌套在这些高阶单元中。例如，如果家族被招募，来自同一家族的兄弟姐妹的观察结果可能比来自其他家族的兄弟姐妹的观察结果更相似，这是由于共有的遗传因素和共有的环境影响。类似地，如果对个人进行一段时间的跟踪，并在一周内每天进行评估，来自同一个人的观察结果自然更有可能彼此相关，而不是与其他人的观察结果相关。

在这一章中，我们将关注一个重复测量数据的例子。尽管重复测量数据可能与家庭内兄弟姐妹或公司内多名员工的数据有很大不同，但它们带来了许多相同的挑战，并有共同的解决方案。主要区别在于，重复测量数据自然地按时间排序，而聚类数据通常是无序的(例如，如果从 100 家公司中抽取 10 名员工，通常没有自然的方法对员工进行排序或排序)。

多级数据的第一个区别是有两种常见的数据组织方式。一种方法，有时称为宽数据集，类似于单级数据。每一行代表一个观察单位(例如，一个人，一所学校)，并添加额外的变量(列)用于单位内的重复测量。例如，假设高血压患者在研究开始时(??)、6 个月后(??)和 1 年随访时(??)纵向测量血压。表 10-1 给出了这些数据的一个示例。在这个例子中，ID 3 错过了最后的时间点。

表 10-1

在三个时间点测量的收缩压(SBP)的示例宽数据集

|

身份

|

SBPT1 型

|

SBPT2 型

|

SBPT3 型

|
| — | — | — | — |
| 1 | One hundred and thirty-five | One hundred and thirty | One hundred and twenty-five |
| 2 | One hundred and twenty | One hundred and twenty-five | One hundred and twenty-one |
| 3 | One hundred and twenty-one | One hundred and twenty-five | 。 |

虽然宽数据在某些情况下很方便，但它通常不是多级数据的理想结构。首先，对于纵向数据，比如表 10-1 中的例子，部分信息实际上编码在变量名中:??，??，??。一般来说，变量名最好只描述变量或测量值，而附加信息(如时间)要被捕获并编码到另一个变量(如时间点)中。第二，如果在一个单元内有不同数量的评估，如不同长度的纵向随访或学校内不同数量的学生，宽格式变得非常低效。例如，想象一个数据集，其中从一所大型学校招聘了 600 名学生，而从一所较小的学校招聘了 15 名学生。要将这些数据以宽格式放置，每个测量需要 600 个变量(最大学校的每个学生一个变量)。在这个小学校里，这 600 个变量中除了 15 个之外，其他都将丢失。

构建多级数据集的另一种方法有时称为长格式。在长格式中，多行可能属于任何一个单元，这由一个标识(ID)变量来表示。表 10-2 显示了一个使用相同的血压假设纵向研究的长数据示例。在这种情况下，如果一个特定的单元错过了一个观测值或者比其他单元的观测值少(例如，一所小学校对一所大学校)，则可以省略这些行。例如，在表 10-2 中，ID 3 没有时间 3 血压读数，因此该行完全不在数据中。

表 10-2

示例在三个时间点测量的收缩压(SBP)的长数据集

|

身份

|

自发性腹膜炎

|

时间

|
| — | — | — |
| 1 | One hundred and thirty-five | one |
| 1 | One hundred and thirty | Two |
| 1 | One hundred and twenty-five | three |
| 2 | One hundred and twenty | one |
| 2 | One hundred and twenty-five | Two |
| 2 | One hundred and twenty-one | three |
| 3 | One hundred and twenty-one | one |
| 3 | One hundred and twenty-five | Two |

重塑数据

如果数据是以一种格式存储的，就有可能，有时也有必要将其转换成另一种格式。下面的代码创建一个宽数据集，然后使用reshape()函数将其从宽格式转换为长格式。varying参数表示随时间变化的宽数据中的变量。v.names参数表示长数据集中每组可变变量的名称。在这个例子中，我们只有 SBP，但如果有更多的变量，这些可以添加。方向参数指示数据是应该由宽到长，还是由长到宽。

ex.wide <- data.table(
  ID = c(1, 2, 3),
  SBPT1 = c(135, 120, 121),
  SBPT2 = c(130, 125, 125),
  SBPT3 = c(125, 121, NA))

print(ex.wide)

##    ID SBPT1 SBPT2 SBPT3
## 1:  1   135   130   125
## 2:  2   120   125   121
## 3:  3   121   125    NA

reshape(
  data = ex.wide,
  varying = list(paste0("SBPT", 1:3)),
  v.names = c("SBP"),
  idvar = "ID",
  direction = "long")

##    ID time SBP
## 1:  1    1 135
## 2:  2    1 120
## 3:  3    1 121
## 4:  1    2 130
## 5:  2    2 125
## 6:  3    2 125
## 7:  1    3 125
## 8:  2    3 121
## 9:  3    3  NA

相反，如果数据已经是长格式，就可以将它们转换成宽格式。下面的代码创建了一个长数据，并再次使用reshape()函数将其转换为宽数据。sep参数指示如何使用宽变量名。在这种情况下，基数是 SBP，然后是分隔符 T，后面是时间:1、2、3。

ex.long <- data.table(
  ID = c(1, 1, 1, 2, 2, 2, 3, 3),
  SBP = c(135, 130, 125, 120, 125, 121, 121, 125),
  Time = c(1, 2, 3, 1, 2, 3, 1, 2))

print(ex.long)

##    ID SBP Time
## 1:  1 135    1
## 2:  1 130    2
## 3:  1 125    3
## 4:  2 120    1
## 5:  2 125    2
## 6:  2 121    3
## 7:  3 121    1
## 8:  3 125    2

reshape(
  data = ex.long,
  v.names = "SBP",
  timevar = "Time",
  sep = "T",
  idvar = "ID",
  direction = "wide")

##    ID SBPT1 SBPT2 SBPT3
## 1:  1   135   130   125
## 2:  2   120   125   121
## 3:  3   121   125    NA

每日数据集

在开始之前，我们将介绍一个新的数据集。这些数据来自 2017 年在莫纳什大学进行的每日日记研究，年轻成年人在大约 12 天的时间里每天完成多达三次(早上、下午和晚上)的测量。因此，每个参与者对数据集贡献了大约 36 次观察。为了保护参与者的保密性和匿名性，这里使用的数据是从原始数据模拟的，但是以这样的方式保留了变量之间的关系和原始数据的大部分特征。

表 10-3 中列出了变量名和每个变量名的简要描述。

模拟数据作为JWileymisc包的一部分，可以使用data()功能加载。

表 10-3

每日日记研究数据中的变量名称列表

|

变量名

|

描述

|
| — | — |
| UserID | 每个人的唯一标识符 |
| SurveyDay | 每次观察发生的日期 |
| SurveyInteger | 调查编码为整数(1 =上午，2 =下午，3 =晚上) |
| SurveyStartTimec11 | 调查开始时间，从上午 11:00 开始以小时为中心 |
| Female | 0 或 1 变量，其中 1 =女性，0 =男性 |
| Age | 参与者的年龄，最高编码为 25 岁 |
| BornAUS | 一个 0 或 1 变量，其中 1 =出生在澳大利亚，0 =出生在澳大利亚以外 |
| SES_1 | 参与者的主观 SES，底部编码为 4，顶部编码为 8 |
| EDU | 教育水平，其中 1 =大学毕业或以上，0 =大学毕业以下 |
| SOLs | 自我报告的睡眠开始潜伏期(分钟),仅限早晨调查 |
| WASONs | 自我报告的睡眠开始后醒来的次数，最高编码为 4，仅早晨调查 |
| STRESS | 总体压力等级为 0-10，每天重复 3 次 |
| SUPPORT | 总体社会支持评分为 0-10 分，每天重复 3 次 |
| PosAff | 正面影响评分为 1-5 分，每天重复 3 次 |
| NegAff | 负面影响评分为 1-5 分，每天重复 3 次 |
| COPEPrb | 问题集中在 1-4 级的应对上，在晚上的调查中每天重复 1 次 |
| COPEPrc | 情绪处理应对等级为 1-4，在晚间调查中每天重复 1 次 |
| COPEExp | 情绪表达应对等级为 1-4，在晚间调查中每天重复 1 次 |
| COPEDis | 精神解脱应对在 1-4 的范围内，在晚上的调查中每天重复 1 次 |

data(aces_daily)
str(aces_daily)

## 'data.frame':        6599 obs. of  19 variables:
##  $ UserID            : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ SurveyDay         : Date, format: "2017-02-24" ...
##  $ SurveyInteger     : int  2 3 1 2 3 1 2 3 1 2 ...
##  $ SurveyStartTimec11: num  1.93e-01 4.86e-01 1.16e-05 1.93e-01 4.06e-01 ...
##  $ Female            : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Age               : num  21 21 21 21 21 21 21 21 21 21 ...
##  $ BornAUS           : int  0 0 0 0 0 0 0 0 0 0 ... 

##  $ SES_1             : num  5 5 5 5 5 5 5 5 5 5 ...
##  $ EDU               : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ SOLs              : num  NA 0 NA NA 6.92 ...
##  $ WASONs            : num  NA 0 NA NA 0 NA NA 1 NA NA ...
##  $ STRESS            : num  5 1 1 2 0 0 3 1 0 3 ...
##  $ SUPPORT           : num  NA 7.02 NA NA 6.15 ...
##  $ PosAff            : num  1.52 1.51 1.56 1.56 1.13 ...
##  $ NegAff            : num  1.67 1 NA 1.36 1 ...
##  $ COPEPrb           : num  NA 2.26 NA NA NA ...
##  $ COPEPrc           : num  NA 2.38 NA NA NA ...
##  $ COPEExp           : num  NA 2.41 NA NA 2.03 ...
##  $ COPEDis           : num  NA 2.18 NA NA NA ... 

数据是长格式的，因此每一行代表一个人在某一天的一次调查中的观察结果。在长格式中，每个人贡献大约 36 行数据。虽然长格式是存储具有许多重复测量值的数据的有效方式，但它会使识别缺失数据变得困难，因为缺失的调查不会被注册为一行缺失数据，而是整行缺失。第一步有助于确定缺失数据的比率，并且对未来的几项分析非常有用，这就是添加任何缺失的调查，只需将观测值设置为缺失即可。为此，我们可以创建一个临时数据集，其中包含从个人的第一次调查/日期到最后一次调查/日期的所有调查和日期，然后对原始数据执行完全连接(或合并，保留所有行)。

下面的代码通过查找最小和最大日期以及这些日期中每一天的最早和最晚调查来创建这样一个临时数据集。然后，使用这些信息，我们可以创建一个“完整”的数据集，包含第一天的第一次调查和最后一天的最后一次调查之间的所有调查和所有天数。最后，我们可以将两个数据集合并在一起，保留所有将根据需要填充缺失值的行。

draw <- as.data.table(aces_daily)
draw <- draw[order(UserID, SurveyDay, SurveyInteger)]
draw[, UserID := factor(UserID)]

tmpdata <- draw[!is.na(SurveyDay) & !is.na(SurveyInteger)][, .(
  MinD = min(SurveyDay),
  MinS = min(SurveyInteger[SurveyDay == min(SurveyDay)]),
  MaxD = max(SurveyDay),
  MaxS = max(SurveyInteger[SurveyDay == max(SurveyDay)])),
  by = UserID]

tmpdata <- tmpdata[, .( 

  SurveyInteger = c(
    MinS:3L, #first day
    rep(1L:3L, times = MaxD - MinD - 1), #all days between first/last
    1L:MaxS), #last day
  SurveyDay = as.Date(rep(MinD:MaxD, c(
      4L - MinS, #first day
      rep(3, MaxD - MinD - 1), #all days between first/last
      MaxS)), origin = "1970-01-01")), #lastday
  by = UserID]

d <- merge(draw, tmpdata, by = c("UserID", "SurveyDay", "SurveyInteger"),
           all = TRUE)

nrow(draw)

## [1] 6599

nrow(d)

## [1] 6927

nrow(draw)/nrow(d)

## [1] 0.95

在添加丢失的行后，数据中的行数从原始数据开始增加，反映了丢失的数据。

10.2 描述性统计

对于非独立数据，基本描述性统计也可以用不同的方法计算。为了开始理解这些差异，以及多级结构意味着什么，检查图 10-1 中的两个图。该图显示了在不同条件下对四个不同的人进行十次评估的虚构数据。实线表示每个人的平均值，圆点表示观察到的数据。

对于长格式的多级数据，如果我们计算一个变量y的平均值和方差，这将是所有人和所有时间的平均值。方差将包括人与人之间的差异(线与线之间的距离)和人与人之间的方差(数据点围绕每个人的平均值变化的程度)。相反，如果我们首先对一个人在不同时间的观察进行平均，那么平均值将是四条线的平均值，方差将只是个体平均值之间的可变性，如图 10-1 中的线。

此外，计算所有数据点的描述符与计算人平均数据的描述符会导致参与者的权重不同。当计算所有观察的汇总时，有十个观察的参与者将得到只有一个观察的参与者的十倍权重(例如，由于缺失数据)。加权的问题往往不太重要，因为所有聚类的大小都差不多，并且如果所有聚类都相同(例如，每个人正好有 10 个观察值)，那么加权也没有什么区别。

set.seed(1234)
ex.data.1 <- data.table(
  ID = factor(rep(1:4, each = 10)),
  time = rep(1:10, times = 4),
  y = rnorm(40, rep(1:4, each = 10), .2))

ex.data.2 <- data.table(
  ID = factor(rep(1:4, each = 10)),
  time = rep(1:10, times = 4),
  y = rnorm(40, 2.5, 1))

plot_grid(
 ggplot(ex.data.1,
        aes(time, y, colour = ID, shape = ID)) +
  stat_smooth(method = "lm", formula = y ~ 1, se=FALSE) +
  geom_point() +
  scale_color_viridis(discrete = TRUE),
 ggplot(ex.data.2,
        aes(time, y, colour = ID, shape = ID)) +
  stat_smooth(method = "lm", formula = y ~ 1, se=FALSE) +
  geom_point() +
  scale_color_viridis(discrete = TRUE),
 ncol = 1,
 labels = c(
   "High Between Variance",
   "Low Between Variance"),
 align = "hv")

对于多级数据，计算描述性统计数据没有简单的正确或错误的方法，但理解差异并准确描述所使用的方法是很重要的。一般来说，三种常见的方法如下:

• 忽略结构并计算所有观察值的描述符。如果它们具有不同数量的观察值，这可能会对单元进行不同的加权。它还将提供一个方差估计值，该值结合了人与人之间以及人与人之间的差异。也就是说，可变性是变量的总可变性。

• 首先平均(或合并)一个单元内的观察值，然后计算描述符。这使得每个单元重量相等。方差估计只会捕获变量在人与人之间的方差，因此它可能更有助于描述样本特征，而不是变量的总可变性。

• 仅计算第一个时间点的描述性统计数据。这将仍然倾向于包括一些个人之间和个人内部的可变性，因为个人内部的可变性没有被平均。如果随着时间的推移出现有意义的变化，平均值可能不代表研究的整体平均值。这种方法只对纵向数据有意义。对于其他多级结构(例如，嵌套在教室中的学生)，没有合理的方法来选择使用哪个学生。

基本描述

计算描述性统计数据的每种方法的例子在我们前面为正面影响加载的每日数据示例中显示。使用data.table()直接在数据集中动态地执行必要的数据管理，而不是创建新的变量或新的数据集。为了通过 ID 获得第一个观察，我们按照 ID、星期几和调查(上午、下午、晚上)进行排序，然后通过 ID 选择第一个观察。

图 10-1

显示高方差和低方差假设数据的图。在高方差中，一个人的观察值变化很小，但个体差异很大。在低间方差中，个体差异不多，但每个人内部的变异性很大。

## mean and SD on all observations
egltable("PosAff", data = d)

##                M (SD)
## 1: PosAff 2.68 (1.07)

## mean and SD first averaging within ID
egltable("PosAff",
  data = d[, .(
    PosAff = mean(PosAff, na.rm = TRUE)) 

,
    by = UserID])

##                M (SD)
## 1: PosAff 2.68 (0.80)

## mean and SD on first observations
egltable("PosAff", data = d[
  order(UserID, SurveyDay, SurveyInteger)][,
    .(PosAff = PosAff[1]), by = UserID])

##                M (SD)
## 1: PosAff 2.71 (1.02) 

对于长数据集，为了计算时不变变量的汇总度量，而不是在参与者中求平均值，我们必须首先对数据进行子集化以删除重复的行，以便它回到单级结构。这可以通过要求数据删除任何重复的 id 来实现。

tab <- egltable(c("Female", "Age", "BornAUS", "SES_1", "EDU"),
                data = d[!duplicated(UserID)],
                strict = FALSE)
tab

##             M (SD)/N (%)
##  1:  Female
##  2:       0    28 (40.6)
##  3:       1    41 (59.4)
##  4:     Age 21.91 (2.38)
##  5: BornAUS
##  6:       0    41 (60.3)
##  7:       1    27 (39.7)
##  8:   SES_1  6.05 (1.21)
##  9:     EDU
## 10:       0    45 (66.2)
## 11:       1    23 (33.8)

绘图是一种显示描述性统计数据的有用方式，例如按不同的组。图 10-2 显示了女性和男性的平均应对能力。平均值标绘为点，并添加了标度中的锚点，以便于解释。

图 10-2

显示男女平均应对评级的图表

## create a dataset of the means and labels by gender
copeplotdata <- d[!is.na(Female), .(
  M = c(
    mean(COPEPrb, na.rm = TRUE),
    mean(COPEPrc, na.rm = TRUE),
    mean(COPEExp, na.rm = TRUE),
    mean(COPEDis, na.rm = TRUE)),
  Var = 1:4,
  Low = sprintf("I usually don’t do this at all\n[%s]",
                c("Problem Focused", "Emotional Processing",
                  "Emotional Expression", "Disengagement")),
  High = sprintf("I usually do this a lot\n[%s]",
                 c("Problem Focused", "Emotional Processing",
                   "Emotional Expression", "Disengagement"))),
  by = Female]

## coded 0/1 but for plotting, R needs to know
## it is discrete not a continuous number
copeplotdata[, Female := factor(Female)]

## create a plot

gglikert(x = "M", y = "Var", leftLab = "Low", rightLab = "High",
         data = copeplotdata, colour = "Female",
  xlim = c(1, 4), title = "Average Coping") +
  scale_colour_manual(values =
    c("1" = "grey70", "0" = "grey30"))

描述性统计也可以被其他变量分解。例如，以下代码根据当时报告的压力水平计算并绘制了图 10-3 中积极和消极影响的平均水平。请注意，由于这是通过调查得出的，因此将在调查而非人员层面进行解释。也就是说，在调查中，人们给自己的压力评分在 5 分以上，他们平均给自己的影响评分是多少？平均来说，它并没有告诉我们压力大或小的人的平均影响。事实上，如果同一个人有时报告的压力高于 5，而有时报告的压力低于 5，那么他可能会将一些调查贡献给高压力的平均值，而将一些调查贡献给低压力的平均值。

## create a dataset of the means and labels by stress
afplotdata <- d[!is.na(STRESS), .(
  M = c(
    mean(PosAff, na.rm = TRUE),
    mean(NegAff, na.rm = TRUE)),
  Var = 1:2,
  Low = sprintf("Very Slightly or\nNot at all\n[%s]",
                c("Positive Affect", "Negative Affect")),
  High = sprintf("Extremely\n\n[%s]",
                c("Positive Affect", "Negative Affect"))),
  by = .(Stress = STRESS > 5)]

## add labels to understand stress
afplotdata[, Stress := factor(Stress, levels = c(FALSE, TRUE),
                              labels = c("<= 5", "> 5"))]

## create a plot
gglikert(x = "M", y = "Var", leftLab = "Low", rightLab = "High",
         data = afplotdata, colour = "Stress",
  xlim = c(1, 5), title = "Affect by Stress") +
  scale_colour_manual(values =
    c("<= 5" = "grey70", "> 5" = "grey30"))

当观察值可以有意义地排序时，比如按一天中的时间排序，我们也可能希望按时间点分别计算描述性统计数据。在下面的代码中，我们为带有更好标签的调查创建了一个新变量，然后对参与者的回答进行平均，但根据调查分别进行，最后计算描述性统计数据。为了得到一个好的分组摘要，我们可以将Survey作为分组变量传递。但是，请注意，在这种情况下，组间差异的统计测试将是不准确的，因为它假设了独立的组。我们忽略测试，专注于仍然准确的描述性统计。这个例子还展示了一次获取多个变量的描述性统计数据。

图 10-3

显示男女平均应对评级的图表

d[, Survey := factor(SurveyInteger, levels = 1:3,
    labels = c("Morning", "Afternoon", "Evening"))]

egltable(c("PosAff", "NegAff", "STRESS"), g = "Survey",
  data = d[, .(
    PosAff = mean(PosAff, na.rm = TRUE),
    NegAff = mean(NegAff, na.rm = TRUE),
    STRESS = mean(STRESS, na.rm = TRUE)
    ), by = .(UserID, Survey)])

##           Morning M (SD) Afternoon M (SD) Evening M (SD)
## 1: PosAff    2.67 (0.84)      2.69 (0.81)    2.67 (0.81)
## 2: NegAff    1.53 (0.46)      1.57 (0.49)    1.56 (0.49)
## 3: STRESS    2.14 (1.47)      2.52 (1.60)    2.39 (1.56)
##                          Test
## 1: F(2, 570) = 0.05, p = .947
## 2: F(2, 570) = 0.33, p = .720
## 3: F(2, 570) = 3.01, p = .050

附加选项(参与者平均与否，按时间点分别报告等。)的出现是因为对于重复测量或非独立数据，观测值可以分解到不同的层次。具体来说，我们可以想象具体的积极情绪观察是参与者的典型或平均积极情绪加上特定日期或时间的影响的组合。

我们可以通过计算每个参与者的平均值(即，中间部分)，然后取观察值和每个参与者自己的平均值之间的差值(即，内部部分)，将变量分解为“中间”和“内部”两个方面。下面的代码展示了一个有积极影响的例子。在计算了这两个分量之后，我们还可以分别得到每个分量的描述性统计数据。注意，对于每个参与者的平均值，我们应该首先删除重复的值；否则，它将更倾向于具有更多重复测量的参与者，这对于级别间或时间不变性变量是不合适的。

d[, BPosAff := mean(PosAff, na.rm = TRUE), by = UserID]
d[, WPosAff := PosAff - BPosAff]

egltable("BPosAff", data = d[!duplicated(UserID)])

##                 M (SD)
## 1: BPosAff 2.68 (0.80)

egltable("WPosAff", data = d)

##                 M (SD)
## 1: WPosAff 0.00 (0.72)

分解变量通常不仅对描述性统计有用，而且对分析人们希望了解变量在人与人之间以及人与人之间的关系也有用。为了便于这样的工作，我们可以对所有时变变量进行变量间和变量内转换。为了减少我们编写的代码量，我们可以定义一个函数，bwmean(),来计算平均值和与平均值的偏差，然后使用data.table()通过 ID 来应用这个函数。

## define a new function
bwmean <- function(x, na.rm = TRUE) {
  m <- mean(x, na.rm = na.rm)
  list(m, x - m)
}

## apply it to affect, support, and stress, by ID
d[, c("BNegAff", "WNegAff") := bwmean(NegAff), by = UserID]
d[, c("BSUPPORT", "WSUPPORT") := bwmean(SUPPORT), by = UserID]
d[, c("BSTRESS", "WSTRESS") := bwmean(STRESS), by = UserID] 

我们还可以为睡眠和应对措施制定一个介于和内部变量。然而，它们的工作方式略有不同。睡眠只在早上测量，应对只在晚上测量，作为当天的整体应对。如果我们通过调查统计非缺失观测值的数量，这一点很容易看出，首先使用is.na()函数根据是否缺失将每个值转换为 0 或 1，然后使用sum()求和，所有这些都通过调查分解。注意is.na()返回 TRUE 或 1 表示缺失，返回 FALSE 或 0 表示不缺失。为了扭转这种情况，我们添加了感叹号，这意味着代码已被读取，并通过调查对非缺失值求和。

d[, .(
  NCope = sum(!is.na(COPEPrb)),
  NSOLs = sum(!is.na(SOLs))),
  by = Survey]

##       Survey NCope NSOLs
## 1: Afternoon     0     0
## 2:   Evening  2090  2097
## 3:   Morning     0     0

尽管睡眠和应对能力只在特定的调查中进行评估，但它们适用于一整天，所以我们可以在其他调查中填写它们。对于参与者之间的变量，很容易做到这一点:我们只需填写所有调查的平均值。填写其他调查时间点更复杂。关键是传递一个值，这个值R将根据需要回收。我们可以通过让data.table()应用 ID 和调查日的操作来实现这一点，然后忽略任何缺失的应对值(每天最多留下一个)，并减去人与人之间的应对变量。因为我们已经填写了，人与人之间的应对变量在每次调查时都会有值，所以内部应对变量也会有值。

d[, BCOPEPrb := mean(COPEPrb, na.rm = TRUE), by = UserID]
d[, WCOPEPrb := na.omit(COPEPrb) - BCOPEPrb,
  by = .(UserID, SurveyDay)]
d[, BCOPEPrc := mean(COPEPrc, na.rm = TRUE), by = UserID]
d[, WCOPEPrc := na.omit(COPEPrc) - BCOPEPrc,
  by = .(UserID, SurveyDay)]
d[, BCOPEExp := mean(COPEExp, na.rm = TRUE), by = UserID]
d[, WCOPEExp := na.omit(COPEExp) - BCOPEExp,
  by = .(UserID, SurveyDay)]
d[, BCOPEDis := mean(COPEDis, na.rm = TRUE), by = UserID]
d[, WCOPEDis := na.omit(COPEDis) - BCOPEDis,
  by = .(UserID, SurveyDay)]

d[, BSOLs := mean(SOLs, na.rm = TRUE), by = UserID]
d[, WSOLs := na.omit(SOLs) - BSOLs,
  by = .(UserID, SurveyDay)]
d[, BWASONs := mean(WASONs, na.rm = TRUE), by = UserID]
d[, WWASONs := na.omit(WASONs) - BWASONs,
  by = .(UserID, SurveyDay)] 

组内相关系数

另一种对多水平模型有用的描述性统计称为组内相关系数或 ICC。在多层次的背景下，ICC 基于将可变性分解为两个来源:个体之间的可变性和个体内部的可变性。计算变量的 ICC 的可靠方法是使用最简单类型的多级模型:只有随机截距的模型，即只包含截距但允许截距随 ID 随机变化的模型。随机截距的方差是两者之间的方差，因为它本质上是个体均值的方差，而剩余方差是个体内部的方差，它不能用个体自身的均值来解释。这两个方差来源共同构成总方差。

)

(10.1)

使用这两个方差来源，我们可以计算个体之间的变异性与总变异性的比率。该比率在 0 和 1 之间变化，并提供关于个体之间发生的总可变性的多少的信息。值为 0 表示所有个体平均值彼此相等，因此所有可变性都发生在个体内部。相反，值为 1 表示在个体中所有值都是相同的，个体之间存在所有可变性。该比率被称为组内相关系数(ICC ),可通过下式计算:

)

(10.2)

ICC 可以通过拟合随机截距模型手动计算，或者更方便地使用iccMixed()函数。iccMixed()需要变量名、ID 变量(或者多个 ID 的变量)的名称和数据集。它返回每个级别的估计方差，称为适马，以及每个级别的方差与总方差的比率，即 ICC。虽然 ICC 通常用于两级结构，但该函数可以推广到更高阶的结构，如学生内部的观察和有班级的学生。如果学生和班级 id 可用，则可以计算每个级别的方差，ICC 将是每个级别的方差与总方差的比率。

iccMixed("NegAff", "UserID", d)

##         Var Sigma  ICC
## 1:   UserID  0.21 0.44
## 2: Residual  0.27 0.56

iccMixed("PosAff", "UserID", d)

##         Var Sigma  ICC
## 1:   UserID  0.63 0.54
## 2: Residual  0.53 0.46

除了提供个体之间或个体内部变异程度的指标外，ICC 还用于计算“有效”样本量。“有效”样本量是对数据提供的独立样本数量的近似估计。例如，如果连续 10 天每天收集 10 个人的数据，则总共有 100 个观察值，但是这些数据不太可能提供与 100 个人测量一次(独立样本)相同的有效信息。

为了更好地理解有效样本量，通常称为 NEffective，考虑两个极端。首先，假设在一个个体中，每个观察都是相同的。比如，想象一下测量一个成年人的身高。一旦在第一天测量了身高，另外 9 天的评估不太可能提供任何进一步的有用信息。在这个例子中，ICC 将是 1:所有的可变性发生在个体之间，个体内部没有可变性。也就是说，成年人有许多不同的身高(可变性之间)，但同一成年人每天的身高基本相同(可变性之内没有)。

在另一个极端，一些变量可能每天都在变化，就像人与人之间一样。想象一下评估同一城市不同成年人的通勤时间。忽略不同的路线，就目前而言，变异可能只是因为不同的日子和交通状况而存在。因此，平均而言，他们的所有通勤时间可能是相同的(可变性之间没有差异)，所有可变性都发生在日常基础上(可变性之内)。在后一个例子中，对 10 个人进行 10 天的测量提供了 100 个观察结果和等同于在一天中对 100 个人进行评估的信息。

计算 NEffective 试图提供关于等效独立样本大小的估计。该公式取决于参与者或真正独立单元的数量，N；个人(单位)人均考核次数，k；和 ICC，如下所示:

)

(10.3)

也可以使用nEffective()功能在R中计算有效样本量。这也有助于凸显国际刑事法院的影响力有多大。ICC 越高，无效性越低。下面的R代码显示了消极和积极影响的 NEffective。尽管观察次数相似，但由于 ICC 的差异，消极和积极影响的效果却大不相同，ICC 由nEffective()函数自动计算。

## number of units
n <- length(unique(d$UserID))

## average observations per unit
k <- nrow(d[!is.na(NegAff)])/n

## effective sample size
nEffective(n, k, dv = "NegAff", id = "UserID", data = d)

##                     Type    N
## 1: Effective Sample Size  420
## 2:     Independent Units  191
## 3:    Total Observations 6389

k <- nrow(d[!is.na(PosAff)])/n
nEffective(n, k, dv = "PosAff", id = "UserID", data = d)

##                     Type    N
## 1: Effective Sample Size  343
## 2:     Independent Units  191
## 3:    Total Observations 6399

10.3 探索和假设

分布和异常值

在本书的开始，我们检查了多种方法来可视化单变量和多变量数据，为分析做准备。对于非独立或重复的测量数据，类似的方法适用，除了它们可以应用于不同的水平或不同的单位(例如，观察值、人与人之间的平均值等)。).为了便于分解和图形化检查，我们可以使用meanDecompose()函数。它使用一个公式界面，左侧是主要变量，右侧是要分解的 id 或其他变量。它的工作方式类似于我们在变量之间和变量内部创建的方式，只是它创建了单独的数据集，因此在变量之间的级别上，没有重复。

tmp <- meanDecompose(PosAff ~ UserID, data = d)
str(tmp, max.level = 1)

## List of 2
##  $ PosAff by UserID  :Classes 'data.table' and 'data.frame': 191 obs. of 2 variables:
##   ..- attr(*, "sorted")= chr "UserID"
##   ..- attr(*, ".internal.selfref")=<externalptr>
##  $ PosAff by residual:Classes 'data.table' and 'data.frame': 6927 obs. of  1 variable:
##   ..- attr(*, ".internal.selfref")=<externalptr>

单独的数据集存储在一个列表中，它们由变量和级别命名。这里的“用户标识”和“剩余”对应于级别之间和级别之内。如前几章所述，我们可以对照正态(或其他)分布检查每个变量。在下面的代码中，我们只关注人与人之间的积极影响值，如图 10-4 所示。

图 10-4

与正态分布相反的人与人之间的积极影响

testdistr(tmp[[1]]$X, varlab = names(tmp)[1],
          extremevalues = "theoretical", robust=TRUE)

我们可以让R遍历所有级别的数据并绘制它们，而不是为每个级别编写代码。我们不单独绘制它们，而是在最后使用plot_grid()将它们组合起来，如图 10-5 所示。

图 10-5

与正态分布相反的个人之间和个人内部的积极影响

plots <- lapply(names(tmp), function(x) {
  testdistr(tmp[[x]]$X, plot = FALSE, varlab = x,
            extremevalues = "theoretical", robust=TRUE)[1:2]
})

do.call(plot_grid, c(unlist(plots, FALSE), ncol = 2))

meanDecompose()功能的另一个特点是我们可以添加更多的级别。例如，我们可以查看参与者之间、参与者内部的每日差异，以及参与者和日期的最终残差。图 10-6 显示了一个例子。

图 10-6

对正态分布的不同程度的负面影响

tmp <- meanDecompose(NegAff ~ UserID + SurveyDay, data = d)
do.call(plot_grid, c(unlist(lapply(names(tmp), function(x) {
  testdistr(tmp[[x]]$X, plot = FALSE, varlab = x,
            extremevalues = "theoretical", robust=TRUE)[1:2]
}), FALSE), ncol = 2))

这些图表(图 10-6 )显示，尽管参与者和日或残差的负面影响是对称的，但人与人之间的负面影响是相当偏斜的。我们可以尝试对数转换。结果如图 10-7 所示。尽管参与者之间负面情绪的偏差有所改善，但它仍然存在。按 ID 和日期划分的级别大致呈正态分布。然而，残差的分布虽然是对称的，但却是稀疏的。这些结果表明，我们应该谨慎假设负面情绪的常态。

图 10-7

不同水平的自然对数对正态分布的负面影响

d[, logNegAff := log(NegAff)]
tmp <- meanDecompose(logNegAff ~ UserID + SurveyDay, data = d)
do.call(plot_grid, c(unlist(lapply(names(tmp), function(x) {
  testdistr(tmp[[x]]$X, plot = FALSE, varlab = x,
            extremevalues = "theoretical", robust=TRUE)[1:2]
}), FALSE), ncol = 2))

时间趋势

除了探索分布和异常值，纵向重复测量数据还有其他有用的诊断。如果数据中没有预期的时间趋势，或者时间趋势不是研究的重点，那么根据经验证明没有时间趋势是有用的。这一点很重要，因为许多分析都假设过程是稳定的(即过程不会随着时间的推移而发生实质性变化)。

一个简单的开始方式是绘制不同时间的方法。为此，我们将使用参与者内部变量，以忽略仅由不同参与者驱动的任何潜在差异。为了一次绘制许多变量，我们可以将数据按变量融合成一个长数据集。结果如图 10-8 所示。

图 10-8

gam 平滑下变量随时间变化的趋势

dt <- d[, .(
  WPosAff = mean(WPosAff, na.rm = TRUE),
  WNegAff = mean(WNegAff, na.rm = TRUE),
  WSTRESS = mean(WSTRESS, na.rm = TRUE),
  WSUPPORT = mean(WSUPPORT, na.rm = TRUE),
  WSOLs = mean(WSOLs, na.rm = TRUE),
  WWASONs = mean(WWASONs, na.rm = TRUE)) , by = SurveyDay]
dt <- melt(dt, id.var = "SurveyDay")

ggplot(dt, aes(SurveyDay, value)) +
  geom_point() +
  stat_smooth(method = "gam", formula = y ~ s(x, k = 10)) +
  facet_wrap(~ variable, scales = "free")

除了系统的时间趋势之外，我们还可以看看一周中不同的日子，或者最常见的工作日与周末的差异。我们可以使用weekdays()函数将日期转换成星期几，然后测试这些日期是否与星期六或星期天相匹配，从而得到一个逻辑比较，显示今天是否是周末。结果如图 10-9 所示。

图 10-9

gam 平滑下变量随时间变化的趋势

dt[, Weekend := weekdays(SurveyDay) %in% c("Saturday", "Sunday")]
ggplot(dt, aes(Weekend, value)) +
  stat_summary(fun.data = mean_cl_boot) +
  facet_wrap(~ variable, scales = "free")

## Warning: Removed 2 rows containing non-finite values (stat_summary).

总的来说，时间趋势和一些工作日与周末的差异都出现了。这些初步结果表明，未来的分析应该调整这种差异，或者我们可以计算新的变量，这些变量是去除任何时间趋势后的残差。

自相关

除了时间趋势之外，观察变量之间的相关程度也很有帮助。这被称为自相关。使用默认自相关工具时，通常假设观测值在时间上间隔相等，并且不会有任何缺失值。现在，作为一种快速而粗略的探索方法，我们将分两步填充缺失的数据。首先，我们将估计任何遗漏的调查开始时间，作为每个人在每个调查时间(上午、下午、晚上)的平均时间。然后我们可以将日期和时间结合起来，得到一个日期和时间变量，如下面的代码所示。

d[, StartTimec11Alt := ifelse(is.na(SurveyStartTimec11),
                              mean(SurveyStartTimec11, na.rm = TRUE),
                              SurveyStartTimec11),
  by = .(UserID, Survey)]
d[, StartDayTimec11Alt := chron(
      dates. = format(SurveyDay, "%m/%d/%Y"),
     times. = StartTimec11Alt)]

大多数自相关函数是为单个时间序列设计的，因此我们将一次操作一个参与者。为了看一个例子，我们从单个参与者的图开始。我们先把数据做成时间序列对象，使用zoo()函数(其中“zoo”代表 Z 的有序观测值，“Z”是作者姓氏的第一个字母)。接下来，我们使用na.approx()函数通过插值来填充缺失值。最后，我们准备使用acf()函数计算自相关。结果如图 10-10 所示，显示积极影响在滞后 0 时(即同一时间点)完全相关，并且显示自相关性在较晚滞后时减小。由于我们每天有三次调查，也有可能滞后 3 的观测值具有更高的自相关性，因为这些观测值代表不同日期的相同时间。然而，数据并不支持这一点，相反，这表明绝对时间的流逝可能是最显著的因素。

图 10-10

一个参与者的自相关

tmpd <- d[UserID == 1]
acf(na.approx(zoo(tmpd$PosAff,
    order.by = tmpd$StartDayTimec11Alt)),
    lag.max = 10)

前面的代码显示了一个参与者的结果，但是我们有很多参与者。接下来的代码使用了一个扩展，acfByID()用来计算 ID 的自相关。我们可以用它来生成一个数据集，用于从 0 到 10 的积极影响和每个 ID 的自相关。我们可以对负面影响和压力重复这个过程，然后使用箱线图可视化结果，如图 10-11 所示。箱线图显示了每个滞后的自相关 id 的分布。在 0(无自相关)和0.5处添加线条，作为相对较强自相关的粗略指示。

acf.posaff <- acfByID("PosAff", "StartDayTimec11Alt",
                      "UserID", d)

print(acf.posaff)

##       UserID Variable Lag AutoCorrelation
##    1:      1   PosAff   0          1.0000
##    2:      1   PosAff   1          0.0016
##    3:      1   PosAff   2          0.1249
##   ---
## 2099:    191   PosAff   8         -0.0199
## 2100:    191   PosAff   9          0.1337
## 2101:    191   PosAff  10          0.1828

## make for other measures
acf.negaff <- acfByID("NegAff", "StartDayTimec11Alt",
                      "UserID", d)
acf.stress <- acfByID("STRESS", "StartDayTimec11Alt",
                      "UserID", d)

## put into one dataset for plotting a panel
acf.all <- rbind(
  acf.posaff, acf.negaff,
  acf.stress)

ggplot(acf.all,
    aes(factor(Lag), y = AutoCorrelation)) +
  geom_hline(yintercept = 0, colour = "grey50", size = 1) +
  geom_hline(yintercept = c(-.5, .5),
             linetype = 2, colour = "grey50", size = 1) +
  geom_boxplot() + ylab("Auto Correlation") +
  facet_wrap(~ Variable, ncol = 1)

图 10-11 的结果表明，滞后 1 秒后，自相关性相当小。因此，对于分析，我们将探索 lag1 预测。为此，我们需要创建包含 lag 1 值的变量，我们在下面的代码中就是这样做的。首先，因为滞后 1 是一个调查差异，我们需要一个新的衡量标准，对每个参与者的调查从开始到结束进行排序。请注意，以下代码仅适用于没有遗漏调查或遗漏调查已经添加的情况，就像我们在本章开始时所做的那样。然后，我们可以用它来计算滞后值。对于应对和睡眠，我们只能滞后一天，因为它们是每天测量的。

最后，我们可以将处理后的数据压缩并保存为 RDS 文件，以便在后面的章节中使用。

图 10-11

所有参与者积极和消极影响和压力的自相关

## ensure data ordered by ID, date, and time
d <- d[order(UserID, SurveyDay, SurveyInteger)]
## calculate a number for the survey from 1 to total
d[, USURVEYID := 1:.N, by = .(UserID)]

d[,
  c("NegAffLag1", "WNegAffLag1",
    "PosAffLag1", "WPosAffLag1",
    "STRESSLag1", "WSTRESSLag1") :=
    .SD[.(UserID = UserID, USURVEYID = USURVEYID - 1),
    .(NegAff, WNegAff,
      PosAff, WPosAff,
      STRESS, WSTRESS),
      on = c("UserID", "USURVEYID")]]

d[,
  c("WCOPEPrbLag1", "WCOPEPrcLag1",
    "WCOPEExpLag1", "WCOPEDisLag1",
    "WSOLsLag1", "WWASONsLag1") :=
  .SD[.(UserID = UserID, Survey = Survey, SurveyDay = SurveyDay - 1),
      .(WCOPEPrb, WCOPEPrc, WCOPEExp, WCOPEDis,
        WSOLs, WWASONs),
      on = c("UserID", "Survey", "SurveyDay")]]

## save data after processing, with compression
## for use in subsequent chapters
saveRDS(d, file = "aces_daily_sim_processed.RDS",
        compress = "xz")

假设

GLMMs 与其他回归模型有相似的假设。他们认为

• 用于随机效果的单位是独立的

• 在预测因素和结果之间有一个线性关系

• 对于正态分布的结果，剩余方差是齐次的

• 随机效应遵循(多元)正态分布

• 结果来自预期分布(正态、泊松等。)

虽然没有很好的独立性测试，但其余的假设可以通过视觉进行评估。对于具有非正态结果的 GLMMs，评估结果的残差和分布可能更棘手。对于具有(至少假设的)正态分布结果的 GLMMs，可以检查残差和拟合与残差值的标准图，以检查残差方差是否同质，以及是否满足所有分布假设。在plotDiagnosticsLMER()功能中，这些诊断图被捆绑在一起。为了说明这些，我们拟合了一个模型，预测人与人之间和人与人之间压力的负面影响，包括随机截距和斜率。我们将等到以后再解释模型本身。目前，重点是检查假设。图 10-12 中的曲线显示残差过于集中，不符合正态分布。

图 10-12

负面影响混合效应模型诊断图，显示残差分布(左上)、残差与用于评估方差同质性的拟合值(右上)、随机截距分布(左中)、随机斜率分布(右中)以及随机效应是否为多元正态(左下)。

m.negaff <- lmer(NegAff ~ 1 + BSTRESS + WSTRESS +
            (1 + WSTRESS | UserID), data = d)
assumptiontests <- plotDiagnosticsLMER(m.negaff, plot = FALSE)
do.call(plot_grid, c(
  assumptiontests[c("ResPlot", "ResFittedPlot")],
  assumptiontests$RanefPlot, ncol = 2))

图 10-13 中的曲线显示，总体而言，正态性假设大致得到满足。确实有一个明显的多元异常值。

图 10-13

正面影响混合效应模型诊断图，显示残差分布(左上)、残差与用于评估方差同质性的拟合值(右上)、随机截距分布(左中)、随机斜率分布(右中)以及随机效应是否为多元正态(左下)。

m.posaff <- lmer(PosAff ~ 1 + BSTRESS + WSTRESS +
            (1 + WSTRESS | UserID), data = d)

assumptiontests <- plotDiagnosticsLMER(m.posaff, plot = FALSE)
do.call(plot_grid, c(
  assumptiontests[c("ResPlot", "ResFittedPlot")],
  assumptiontests$RanefPlot, ncol = 2))

为了进一步探索多元异常值，我们可以检查极值。这些按类型分开:Residuals, Random Effect UserID : WSTRESS和Multivariate Random Effect UserID。根据这些图，我们可能不太担心残差，因为我们有如此大的样本量。我们可以检查出多元异常值。在这种情况下，我们看到 id57 和 123 似乎是罪魁祸首。要查看没有这些多元异常值影响的结果，我们可以删除它们并重新估计模型。图 10-14 中的曲线显示，总体而言，正态假设大致得到满足，不再有任何单个随机效应的明显异常值，也不再有多元异常值。

图 10-14

正面影响混合效应模型诊断图，显示残差分布(左上)、残差与用于评估方差同质性的拟合值(右上)、随机截距分布(左中)、随机斜率分布(右中)以及随机效应是否为多元正态(左下)。结果去除两个多元异常值，IDs 57 和 123。

assumptiontests$ExtremeValues[
  EffectType == "Multivariate Random Effect UserID"]

##     PosAff UserID                        EffectType
##  1:    4.7    123 Multivariate Random Effect UserID
##  2:    3.9    123 Multivariate Random Effect UserID
##  3:    3.8    123 Multivariate Random Effect UserID
## ---
## 20:    3.7    123 Multivariate Random Effect UserID
## 21:    4.9    123 Multivariate Random Effect UserID
## 22:    4.6    123 Multivariate Random Effect UserID

m.posaff <- lmer(PosAff ~ 1 + BSTRESS + WSTRESS +
            (1 + WSTRESS | UserID),
            data = d[!UserID %in% c(57, 123)])

assumptiontests <- plotDiagnosticsLMER(m.posaff, plot = FALSE)
do.call(plot_grid, c(

  assumptiontests[c("ResPlot", "ResFittedPlot")],
  assumptiontests$RanefPlot, ncol = 2))

为了更好地理解这些多元异常值的不同寻常之处，我们还可以将他们的数据与其他人在积极情感和压力方面的数据进行比较。一种方法是为数据集中的每个人绘制压力和积极影响关系的斜率，然后单独突出显示极端情况下的关系。我们也可以绘制数据点，只是为了极端的 id，以确保单个极端的观察不会驱动效果。这在图 10-15 中完成，该图显示 IDs 57 和 123 确实是极端的，尽管观察值似乎与估计的斜率一致。

图 10-15

积极的情感和压力联系突出了极端的情况。

ggplot() +
  stat_smooth(aes(WSTRESS, PosAff, group = UserID),
    data = d[!UserID %in% c(123)], method = "lm",
   se = FALSE, colour = "grey50") +
  stat_smooth(aes(WSTRESS, PosAff, group = UserID),
    data = d[UserID %in% c(123)], method = "lm",
   se = FALSE, colour = "blue", size = 2) +
  geom_point(aes(WSTRESS, PosAff),
    data = d[UserID %in% c(123)], colour = "blue", size = 2) +
  stat_smooth(aes(WSTRESS, PosAff, group = UserID),
  data = d[UserID %in% c(57)], method = "lm",
 se = FALSE, colour = "orange", size = 2) +
geom_point(aes(WSTRESS, PosAff),
  data = d[UserID %in% c(57)], colour = "orange", size = 2)

10.4 总结

本章介绍了多级数据结构以及它们通常是如何格式化或存储的(宽或长)。它展示了在处理多级数据和单级数据时，描述性统计的基本数据探索和报告会发生怎样的变化。它还展示了如何通过一次绘制一个级别的数据来可视化和评估广义线性混合模型(GLMMs)的一些常见假设。最后，它涵盖了纵向数据的具体情况，在这种情况下，额外的诊断，如是否有一致的时间趋势和变量本身随时间的自相关是重要的。表 10-4 总结了所使用的功能。

表 10-4

本章中描述的关键功能列表及其功能摘要

|

功能

|

它的作用

|
| — | — |
| [] | 在data.table对象中执行数据管理，以计算缺失的数据，并对所有数据执行操作，或首先跨 id 合并，然后计算结果；处理多级数据时常用的操作 |
| acfByID() | 按 ID 计算各种滞后的自相关系数，并返回适用于汇总或绘图的数据集 |
| egltable() | 计算描述性统计数据，可选地按另一个变量的分组级别 |
| gglikert() | 显示描述性统计数据，如平均值，响应锚点位于图表的左侧和右侧 |
| iccMixed() | 计算变量的组内相关系数 |
| meanDecompose() | 获取一个变量，并将其分解为不同级别的平均值和残差，以便在多级数据的每个级别快速绘图 |
| nEffective() | 计算多级数据中变量的有效样本量 |
| plotDiagnosticsLMER() | 为线性混合模型创建各种诊断图 |
| reshape() | 将数据从宽格式调整为长格式或从长格式调整为宽格式 |