R-统计编程和数据建模高级教程-一-R 统计编程和数据建模高级教程（一）一、单变量数据可视化本书其余部分讨论的大多数

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

原文：Advanced R Statistical Programming and Data Models

协议：CC BY-NC-SA 4.0

一、单变量数据可视化

本书其余部分讨论的大多数统计模型对数据和最佳模型做出假设。作为数据分析师，我们经常必须指定我们假设数据来自的分布。异常值，也称为极端值或异常值，也可能对许多统计模型的结果产生不适当的影响。在这一章中，我们研究了一次探索一个变量(即单变量)的分布和异常值的视觉和图形方法。本章的目标并不是专门创建漂亮的或出版物质量的图表，也不是显示结果，而是使用图表来理解变量的分布并识别异常值。本章重点关注单变量数据可视化，下一章将采用一些相同的概念，但应用于多变量分布，并涵盖如何评估变量之间的关系。

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(MASS)
library(JWileymisc)
library(data.table)

options(width = 70, digits = 2)

ggplot2包[109]创建了优雅的图形，而cowplot包是一个使图形更整洁的插件[117]。MASS包提供了测试不同分布如何拟合数据的函数。JWileymisc包是由本文作者之一维护的，它提供了各种各样的功能，让我们可以专注于本章中的图形。这个data.table包将会被大量用于数据管理。

1.1 分销

可视化观察到的分布

许多统计模型要求指定变量的分布。直方图使用条形来绘制分布图，可能是最常用的显示单个变量分布的图形。虽然相对较少，但堆积点图是另一种方法，它提供了一种精确的方式来可视化显示各个数据点的数据分布。最后，密度图也很常见，它用一条线来表示在任何给定值下的近似密度或数据量。

堆积点图和直方图

点状图的工作原理是为每个观察到的数据值绘制一个点，如果两个点落在另一个点上，它们就被堆叠起来[118]。与直方图或密度图相比，点阵图的独特之处在于，它们实际绘制的是原始的单个数据点，而不是汇总或汇总它们。这使得点状图成为观察变量分布或扩散的一个很好的起点，尤其是当你的观察值相对较少的时候。

绘制单个数据点的粒度方法也有散点图的局限性。即使是中等规模的数据集(例如几百个)，绘制单个值也是不切实际的。对于数千或数百万的观察值，点状图在可视化总体分布方面甚至更不有效。

我们可以使用ggplot2轻松创建一个图，结果如图 1-1 所示。

img/439480_1_En_1_Fig1_HTML.png

图 1-1

旧车每加仑英里数的堆积点图

ggplot(mtcars, aes(mpg)) +
  geom_dotplot()

## 'stat_bindot()' using 'bins = 30'. Pick better value with 'binwidth'.

简单地说，ggplot2的大部分代码遵循以下代码片段所示的格式。在我们的例子中，我们想要一个点状图，所以几何对象，或“geom”，是一个点状图(geom_dotplot())。有许多优秀的在线教程和书籍可以学习如何使用ggplot2包来制作图形，所以我们在这里不再对ggplot2做更多的介绍。特别是，开发ggplot2的 Hadley Wickham 最近更新了关于该包的书， ggplot2:用于数据分析的优雅图形【109】，这是一个极好的指南。对于那些喜欢更少的概念背景和更多的食谱的人，我们推荐温斯顿·张的 R 图形食谱。

  ggplot(the-data, aes(variable-to-plot)) +
    geom_type-of-graph()

与绘制原始数据的点图不同，直方图是一个条形图，其中条形的高度是在条形宽度指定的范围内的值的计数。您可以改变条形的宽度，以控制在一个条形中聚合和计数的相邻值的数量。较窄的条形聚合了较少的数据点，提供了更精细的视图。较宽的条形聚集更多，并提供更宽的视图。图 1-2 显示了著名的 iris 数据集中花的萼片长度分布的直方图。

img/439480_1_En_1_Fig2_HTML.png

图 1-2

来自虹膜数据的萼片长度直方图

ggplot(iris, aes(Sepal.Length)) +
  geom_histogram()

## 'stat_bin()' using 'bins = 30'. Pick better value with 'binwidth'.

如果您知道分布的形状(例如，正态分布)，您可以检查变量的直方图是否看起来像您认识的分布的形状。在萼片长度数据的情况下，它们看起来近似正态分布，尽管它们显然不是完美的。

如果数据看起来不符合我们希望的分布(如正态分布)，通常会对原始数据进行转换。同样，直方图是检查变换后分布情况的有用方法。图 1-3 显示了一年一度的加拿大猞猁捕获的直方图。从图中我们可以看到变量是正偏的(有一个长的右尾巴)。

img/439480_1_En_1_Fig3_HTML.png

图 1-3

加拿大猞猁年度捕获直方图

ggplot(data.table(lynx = as.vector(lynx)), aes(lynx)) +
  geom_histogram()

## 'stat_bin()' using 'bins = 30'. Pick better value with 'binwidth'.

对于正偏差，平方根或对数变换有助于减少正偏差，并使变量更接近正态分布(假设没有负值)。图 1-4 显示了自然对数变换后的猞猁诱捕直方图。

img/439480_1_En_1_Fig4_HTML.png

图 1-4

自然对数变换后的加拿大猞猁年度捕获直方图

ggplot(data.table(lynx = as.vector(lynx)), aes(log(lynx))) +
  geom_histogram()

## 'stat_bin()' using 'bins = 30'. Pick better value with 'binwidth'.

密度图

将观察到的数据分布可视化的另一个常用工具是绘制经验密度图。除了用geom_density()代替geom_histogram()外，ggplot2的代码与直方图的代码相同。代码如下，结果如图 1-5 所示。

img/439480_1_En_1_Fig5_HTML.png

图 1-5

这是我们萼片长度的密度图

ggplot(iris, aes(Sepal.Length)) +
  geom_density()

经验密度图包括一定程度的平滑，因为对于连续变量，在任何特定值都不会有很多观察值(例如，可能没有观察值为 3.286，即使有值 3.281 和 3.292)。经验密度图通过应用某种程度的平滑来显示分布的总体形状。有时，调整平滑度会有助于查看更粗糙(更接近原始数据)或更平滑(更接近“分布”)的图表。使用adjust参数在ggplot2中控制平滑。我们在图 1-5 中看到的默认值是adjust = 1。小于 1 的值“更嘈杂”或平滑度较低，而大于 1 的值会增加平滑度。我们比较和对比了图 1-6 中的噪声和图 1-7 中的非常平滑。

img/439480_1_En_1_Fig7_HTML.png

图 1-7

非常平滑的密度图

img/439480_1_En_1_Fig6_HTML.png

图 1-6

噪声密度图

ggplot(iris, aes(Sepal.Length)) +
  geom_density(adjust = .5)

ggplot(iris, aes(Sepal.Length)) +
  geom_density(adjust = 5)

将观察分布与预期分布进行比较

虽然检查观察到的数据分布是有帮助的，但是我们经常检查分布，看它是否满足我们希望应用的统计分析的假设。例如，线性回归假设数据(有条件地)是正态分布的。如果经验分布非常不正态或更接近不同的分布，那么使用正态线性回归可能是不合适的。

Q-Q 图

为了评估数据是否符合或接近特定的预期分布，我们可以使用分位数-分位数或 Q-Q 图。Q-Q 图绘制了观察到的数据分位数与来自预期分布(例如，正态分布、贝塔分布等)的理论分位数之间的关系。).Q-Q 图可以用来检查数据是否来自几乎任何分布。因为正态或高斯分布是最常见的，所以ggplot2中制作 Q-Q 图的默认函数默认为正态分布。在 Q-Q 图中，如果数据完全符合预期分布，那么这些点将落在一条直线上。以下代码创建了图 1-8 。

img/439480_1_En_1_Fig8_HTML.png

图 1-8

正常数据看起来像一条直线。萼片。长度似乎相当正常

ggplot(iris, aes(sample = Sepal.Length)) +
  geom_qq()

为了更好地理解geom_qq()是如何工作的，我们可以自己做一个。R内置了许多概率分布的基本函数。这允许人们生成随机数(例如，rnorm())，获得来自给定分布的观察值落在某个值之上或之下的概率(例如，pnorm())，计算来自分布的分位数(例如，qnorm())，以及获得特定值的分布密度(例如，dnorm())。按照惯例，它们被命名为r、p、q、d，后面是发行版名称(或其缩写，如“norm”代表 normal)。应用这一知识，我们使用qnorm()获得以下分位数，其中均值= 0 且标准差为 1 的正态分布值的 10% (.10)将位于该分位数中:

qnorm(p = .1, mean = 0, sd = 1)

## [1] -1.3

如何将此应用于具有不同均值或标准差的正态分布是很简单的。假设我们有三个数据点。如果它们是正态分布的，我们可能会期望中间点以 0.5 的概率落在正态分布上，而另外两个点大约一半落在 0 和 0.5 或 0.5 和 1 之间(即 0.25 和 0.75)。我们可以很容易地获得这些概率的正态分位数。

qnorm(p = c(.25, .50, .75), mean = 0, sd = 1)

## [1] -0.67 0.00 0.67

为了帮助得出适当间隔的概率，我们可以使用ppoints()函数。

ppoints(n = 3, a = 0)

## [1] 0.25 0.50 0.75

ppoints()默认为小调整，而不是完全间隔。对于十个或十个以下的数据点，(1:N - 3/8)/(n + 1 - 2 * 3/8)，对于十个以上的数据点，(1:N - 1/2)/(n + 1 - 2 * 1/2)。无论哪种方式，想法都是一样的。

ppoints(n = 3)

## [1] 0.19 0.50 0.81

剩下的工作就是对我们的数据进行分类，并绘制出理论上的正态分位数。添加平均值和标准偏差在技术上是不必要的；它们是线性调整。无论哪种方式，这些点应该落在一条直线上，但使用它们会使理论分位数具有与我们的原始数据相同的均值和规模。

这里我们使用ggplot2中的qplot()函数来绘图。注意，qplot()中的q代表“快速”,因为它是使用用于更高级图形的ggplot()功能的较长规范的简写。所有这些都是说q与分位数无关。最后，为了帮助可视化，我们使用geom_abline()添加一条斜率为 1、截距为 0 的直线。该函数因一条直线的普通方程而得名，作为 x 的函数:

$b\ast x+a$

(1.1)

其中参数名为截距(a)和斜率(b)。我们在图 1-9 中显示了结果。

img/439480_1_En_1_Fig9_HTML.png

图 1-9

在 x 轴上显示理论标准(基于平均值和标准偏差的预测)

qplot(
  x = qnorm(
    p = ppoints(length(iris$Sepal.Length)),
    mean = mean(iris$Sepal.Length),
    sd = sd(iris$Sepal.Length)),
  y = sort(iris$Sepal.Length),
  xlab = "Theoretical Normal Quantiles",
  ylab = "Sepal Length") +
  geom_abline(slope = 1, intercept = 0)

在这种情况下，我们可以看到数据呈合理的正态分布，因为所有点都非常接近法线，并且围绕该线大致对称。

虽然测试数据是否符合正态分布是常见的，但实际数据可能更接近许多其他分布。我们可以使用geom_qq()通过指定期望分布的分位数函数来绘制 Q-Q 图。例如，回到猞猁诱捕数据，图 1-10 用对数正态分布(qlnorm())评估原始猞猁数据的拟合度。

img/439480_1_En_1_Fig10_HTML.png

图 1-10

测试 lynx 数据是否符合对数正态分布

当使用不常用的分布时，有时ggplot2不知道默认情况下如何选择分布的参数。如果需要，我们可以将预期分布的参数作为命名列表传递给dparams参数。以下示例测试 lynx 数据是否符合图 1-11 中的泊松分布。

img/439480_1_En_1_Fig11_HTML.png

图 1-11

测试 lynx 数据是否符合泊松分布

ggplot(data.table(lynx = as.vector(lynx)), aes(sample = lynx)) +
  geom_qq(distribution = qlnorm)

ggplot(data.table(lynx = as.vector(lynx)), aes(sample = lynx)) +
  geom_qq(distribution = qpois, dparams = list(lambda = mean(lynx)))

密度图

检验观察分布是否与预期分布一致的另一种方法是绘制经验密度与预期分布密度的对比图。为此，我们可以使用geom_density()来绘制经验密度，然后使用stat_function()函数，这是绘制任何函数的通用方法。如果我们绘制函数dnorm()，它将绘制一个正常密度。我们只需要指定正态分布的均值和标准差应该基于我们的数据。结果如图 1-12 所示。同样，数据似乎接近正态分布，尽管并不完美。

img/439480_1_En_1_Fig12_HTML.png

图 1-12

正常曲线和我们的密度图(默认平滑度为 1)

ggplot(iris, aes(Sepal.Length)) +
  geom_density() +
  stat_function(fun = dnorm,
                args = list(
                  mean = mean(iris$Sepal.Length),
                  sd = sd(iris$Sepal.Length)),
                colour = "blue")

虽然绘制经验密度和预期密度不会提供 Q-Q 图未捕捉到的任何信息，但更直观的方法是“看到”彼此上下的分布，而不是看到彼此相对绘制的两个分布。

拟合更多分布

通过 Q-Q 图或观察到的和预期的密度图，我们可以评估许多不同的分布。然而，对于正态分布之外的分布，通常需要指定它们的参数，以获得分位数或密度。可以手动计算参数，并将其传递给适当的分位数或密度函数，但是使用MASS包中的fitdistr()函数，我们可以拟合许多分布，并让 R 通过指定分布的名称来估计参数。目前，fitdistr()支持以下发行版:

贝塔
柯西
卡方检验
指数的
F
微克
几何学的
对数正态
符号逻辑的
负二项式
标准
泊松
t
（统计学家）威伯尔（或韦布尔）

尽管这远不是统计分布的详尽列表，但对于几乎所有使用的统计数据来说，这已经足够了，并且涵盖了本书讨论的统计分析中使用的所有分布。

为了了解如何使用fitdistr()，我们从 beta 发行版中虚构了一些随机数据。贝塔分布对比例很有用，因为贝塔分布以 0 和 1 为界。我们使用set.seed()使我们的例子可重复。

set.seed(1234)
y <- rbeta(150, 1, 4)
head(y)

## [1] 0.138 0.039 0.111 0.099 0.377 0.384

fitdistr()函数将数据、表示分布名称的单变量字符串和分布参数的初始值作为一个列表。

y.fit <- fitdistr(y, densfun = "beta",
                  start = list(shape1 = .5, shape2 = .5))

从fitdistr()我们可以得到分布的估计参数。我们明确地提取它们。

y.fit

##   shape1    shape2
##   1.08      4.28
##  (0.11)    (0.52)

y.fit$estimate["shape1"]

## shape1
##    1.1

y.fit$estimate["shape2"]

## shape2
##    4.3

我们还可以提取对数似然性(通常缩写为 LL ),它告诉我们数据有多大可能来自具有这些参数的分布。可能性越高，表示测试的分布和数据之间的匹配越接近。需要注意的是，除了对数似然或 LL 之外，通常还会报告–2 *LL，通常简称为 2LL。最后，更复杂的模型通常(尽管不总是)至少能稍微更好地拟合数据。为了说明这一点，您可以评估用于给定可能性的自由度。logLik()函数提取对数似然和自由度。

logLik(y.fit)

## 'log˽Lik.' 95 (df=2)

虽然可能性值不容易单独解释，但它们对于比较非常有用。如果我们拟合两个分布，提供较高(对数)可能性的分布更适合数据。我们可以再次使用fitdistr()来拟合正态分布，然后将 beta 分布的 LL 与正态分布的 LL 进行比较。

y.fit2 <- fitdistr(y, densfun = "normal")
logLik(y.fit2)

## 'log˽Lik.' 67 (df=2)

在自由度相同的情况下，β(LL = 95.4)的 LL 高于正态分布(LL = 67.3)。这些结果表明我们应该为这些数据选择 beta 分布。

JWileymisc包提供了一种自动拟合多种分布的方法，并通过testdistr()功能自动查看密度或 Q-Q 图。对于正态分布，只需要很少的R代码就可以完成，结果如图 1-13 所示。

img/439480_1_En_1_Fig13_HTML.png

图 1-13

叠加正态分布的密度图和正态 Q-Q 图

testdistr(y)

为了比较正态分布和贝塔分布的拟合程度，我们可以拟合两者，并使用cowplot包中的plot_grid()函数将两个图形绘制成一个图形面板。图 1-14 中的结果显示了数据与贝塔分布的良好一致性(尽管这并不奇怪，因为我们是从贝塔分布中产生的数据！)以及与正常人的不匹配。请注意，来自密度函数的警告消息很常见，在这种情况下不必担心。

img/439480_1_En_1_Fig14_HTML.png

图 1-14

显示了叠加 beta 或正态分布的密度图以及 Q-Q 图拟合

test.beta <- testdistr(y, "beta",
                       starts = list(shape1 = .5, shape2 = .5),
                       varlab = "Y", plot = FALSE)

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

test.normal <- testdistr(y, "normal", varlab = "Y", plot = FALSE)

plot_grid(
    test.beta$DensityPlot, test.beta$QQPlot,
    test.normal$DensityPlot, test.normal$QQPlot,
    ncol = 2)

对于离散分布，如计数，testdistr()绘制了一种略有不同的图，旨在更好地显示观察到的比例与理论分布的概率质量函数。具体来说，密度值是观察到的比例，然后是给定分布中每个值的预期概率。

举例来说，首先我们模拟负二项分布的一些数据，然后在图 1-15 中绘制假设泊松分布的结果，在图 1-16 中绘制假设负二项分布的结果。这种比较显示，就对数似然性和与期望值的偏差而言，负二项式比泊松更适合数据。

img/439480_1_En_1_Fig16_HTML.png

图 1-16

观察到的离散比例，负二项分布的理论概率用蓝色标出

img/439480_1_En_1_Fig15_HTML.png

图 1-15

观察到的离散比例和泊松分布的理论概率用蓝色标出

set.seed(1234)
y <- rnbinom(500, mu = 5, size = 2)
testdistr(y, "poisson")

testdistr(y, "nbinom")

1.2 异常值

异常值是不同于其他值的值，或者在某些方面不标准或不典型。异常值通常也称为异常值或极值。很难精确定义什么是异常值，但通常它们是以某种方式与大多数人不一致的数据点，通常是以一种相对极端的方式。

对于来自正态分布的数据，异常值的常见阈值是 z 分数为 3 之外的任何值。这些阈值基于假设正态分布的概率。具体来说，如果数据呈正态分布，大约 0.10%的数据将低于 z 值 3，大约 0.10%的数据将高于 z 值+3。使用pnorm()函数，确切的概率如下。

pnorm(c(-3, 3))

## [1] 0.0013 0.9987

因为这些阈值基于正态分布，所以它们不一定有意义地应用于非正态分布的数据。虽然许多统计分析，如线性回归，假设结果是(有条件的)正态分布，很少关于预测的分布假设。然而，预测变量的异常值，特别是当它们处于极端时，会强烈影响结果。

视觉识别异常值通常比使用定量标准定义它们更容易。例如，图 1-17 显示了两个图表。两个图都有三个相对异常的点，值为 5。然而，这些点在画面 A 中可能显得更不合适，在画面 A 中，所有其他数据点形成或多或少连续的组，并且与异常点之间存在相对较大的间隙。即使图 B 中的异常点也有间隙，因为数据点中还有其他间隙，但有几个数据点与其他数据点分离并不奇怪——分离似乎是图 B 中的一种模式，而图 a 中没有。

img/439480_1_En_1_Fig17_HTML.png

图 1-17

显示带有异常值的堆积点图的面板图

set.seed(1234)
d <- data.table(
  y1 = rnorm(200, 0, 1),
  y2 = rnorm(200, 0, .2) + rep(c(-3, -1, 1, 3), each = 50))

plot_grid(
  qplot(c(d$y1, rep(5, 3)), geom = "dotplot", binwidth = .1),
  qplot(c(d$y2, rep(5, 3)), geom = "dotplot", binwidth = .1),
  ncol = 1, labels = c("A", "B"))

根据分布的形状，定义异常值也很困难。图 1-18 显示了两种分布。图 A 是伽马分布，显示了伽马分布的特征性长右尾。即使只有少数几个是相当极端的，数据中也没有明显的“间断”,这是一种典型的连续、长右尾分布类型。此外，从数据中很容易看出，有一种减少频率但相当极端的正值的模式。相反，面板 B 中的正态分布更加对称，没有证据表明存在如此长的尾巴。添加到 B 图中的一两个极端正值可能确实看起来“不正常”

img/439480_1_En_1_Fig18_HTML.png

图 1-18

面板图显示从伽玛和正态分布中随机生成的(没有添加异常值)数据

set.seed(1234)
d2 <- data.table(
  y1 = rgamma(200, 1, .1),
  y2 = rnorm(200, 10, 10))

plot_grid(
  qplot(d2$y1, geom = "dotplot", binwidth = 1),
  qplot(d2$y2, geom = "dotplot", binwidth = 1),
  ncol = 1, labels = c("A", "B"))

这些不同的例子突出了准确定义什么是或不是异常值的困难。虽然我们不能提供任何单一的规则来遵循，但是有一些额外的工具内置在testdistr()函数中来帮助做出这些判断。extremevalues参数可用于指定是否应突出显示低于或高于指定分布的经验数据或理论百分点的值。图 1-19 显示了一个示例，根据经验位置突出显示了最低的 1%和最高的 1%的点。地毯中的密度图(密度曲线下方的线条)通过将它们涂成纯黑色来突出显示。在 Q-Q 图中，极值点也是实心黑色，而不是灰色。如果没有点被涂成纯黑色，这将表明没有点落在第 1 个和第 99 个经验百分位数之外。当使用经验值时，除非有大型数据集可用，否则仔细检查顶部和底部 1%这样的值可能是合理的。然而，在可能的情况下，通常需要一个比顶部和底部 1%更极端的阈值，例如顶部和底部 0.10%，以确保这些值确实不太可能仅由于偶然因素而出现。

img/439480_1_En_1_Fig19_HTML.png

图 1-19

显示突出显示极值的图表

testdistr(d$y1, extremevalues = "empirical",
          ev.perc = .01)

除了根据经验定义极值之外，还可以根据特定的理论分布来定义极值。也就是说，我们可以查看是否有任何值超出了正态分布的底部或顶部 0.10%(图 1-20 )，或者对于相同的数据，是否有任何点是基于伽马分布的百分位数的极值(图 1-21 )。这些图表显示，对于相同的数据，如果我们预期数据遵循正态分布，则一些点可能被视为异常或异常值，但是如果我们预期数据遵循伽玛分布，则这些点可能不是异常的。在实践中，作为数据分析师，我们必须判断任何数据值在多大程度上极端或明显异常，以及如何最终对其进行分析。

img/439480_1_En_1_Fig21_HTML.png

图 1-21

基于伽马分布的理论百分位数突出显示极值的图表

img/439480_1_En_1_Fig20_HTML.png

图 1-20

基于正态分布的理论百分位数突出显示极值的图表

testdistr(d2$y1, "normal", extremevalues = "theoretical",
          ev.perc = .001)

testdistr(d2$y1, "gamma", extremevalues = "theoretical",
          ev.perc = .001)

可能有几个异常值。然而，当使用理论分布时，如果一个异常值比另一个异常值更极端，则不太极端的值可能会被“掩盖”,因为参数估计会受到更极端的值的影响。图 1-22 显示了一个例子，其中有两个异常值:100 和 1000。值 100 被值 1，000 所掩盖，因为理论正态分布的均值和方差被值 1，000 拉高了很多，以至于值 100 不再是异常的。

img/439480_1_En_1_Fig22_HTML.png

图 1-22

图表显示异常值 100 被更极端的异常值 1000 所掩盖

testdistr(c(d2$y2, 100, 1000), "normal",
          extremevalues = "theoretical",
          ev.perc = .001)

如果有多个异常值，可以使用迭代过程，处理最极端的值，然后重新检查，直到不再有异常值出现。然而，通过使用稳健的方法，可以在一定程度上简化该过程。当均值和(共)方差是感兴趣的参数时，一种这样的稳健方法是最小协方差行列式(MCD)估计器。简而言之，MCD 估计器找到原始病例的子集，该子集具有其样本协方差矩阵的最低行列式[82]。在单变量情况下，这相当于具有较低方差的原始数据情况的子集。testdistr()函数有一个可选的robust参数，可用于正态分布。当robust = TRUE时，testdistr()使用robustbase包【59】中的covMcd()函数，该函数使用【83】提出的快速算法来计算(近似)MCD。使用稳健估计器的结果如图 1-23 所示。使用稳健估计器，甚至在去除更极端的异常值之前，两个异常值都被识别。

img/439480_1_En_1_Fig23_HTML.png

图 1-23

基于稳健估计的突出显示极值的图形

testdistr(c(d2$y2, 100, 1000), "normal",
          robust = TRUE,
          extremevalues = "theoretical",
          ev.perc = .001)

最后，如果发现异常值，有几种方法可以解决它们。如果可能的话，最好先检查这些值是否准确。异常值通常是由于编码或数据输入错误引起的。如果不存在错误或者无法检查，最简单的选择是排除具有异常值的情况。当异常值的事例很少，并且在排除这些异常值后，大型数据集仍有许多事例时，排除或移除这些事例可能特别有效。在每个案例都很重要的较小数据集中，排除异常值可能会导致太多数据丢失。这也可能发生在更大的数据集中，其中更多的情况是异常的。

排除病例的另一种方法是 winsorizing，以查尔斯·温索尔命名。Winsorizing 采用异常值并用最接近的非异常值替换它们[92，第 14-20 页]。自动完成这一任务的一种方法是计算所需的经验分位数，并将这些值之外的任何值设置为计算的百分位数。即使异常值只存在于分布的一个尾部，这个过程也同样适用于较低和较高的尾部。调整两个尾部有助于确保程序本身不会使分布位置变得更低或更高。Winsorizing 很容易在R中使用JWileymisc包中的winsorizor()函数来完成。唯一需要的参数是每个尾部 winsorize 的比例。winsorizor()函数的另一个特性是，除了返回 winsorized 变量之外，它还添加了注意用于 winsorizing 的阈值和百分点的属性。

winsorizor(1:10, .1)

##  [1] 1.9 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 9.1
## attr(,"winsorizedValues")
##   low high percentile
## 1 1.9  9.1        0.1

图 1-24 比较了我们之前看到的伽马分布变量(在图 A 中)和之后(在图 B 中)的上下 1%。图 B 显示了 winsorizing 的特征“变平”,原始值现在只有 49.2，而不是 58.4。

img/439480_1_En_1_Fig24_HTML.png

图 1-24

面板图比较了(A)和(B)将(经验)底部和顶部的 1%进行 winsorizing 之前和之后的数据

plot_grid(
  testdistr(d2$y1, "gamma", extremevalues = "theoretical",
            ev.perc = .005, plot=FALSE)$QQPlot,
  testdistr(winsorizor(d2$y1, .01), "gamma", extremevalues = "theoretical",
            ev.perc = .005, plot=FALSE)$QQPlot,
  ncol = 2, labels = c("A", "B"), align = "hv")

1.3 摘要

在这一章中，我们学习了使用R可视化原始或混合格式数据的各种方法。此外，除了图形探索性数据分析，我们还学习了一些方法来量化我们的数据与各种分布的拟合。关于本章使用的关键功能的总结，请参见表 1-1 。

表 1-1

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

功能

它的作用

| | --- | --- | | ggplot() | 使用 ggplot2 包创建新地块 | | qplot() | 具有更简单(和更少细微差别)语法的“快速”绘图 | | geom_dotplot() | 创建点状图几何对象-显示所有原始数据 | | geom_histogram() | 创建直方图几何对象 | | geom_density() | 创建密度分布，本质上是平滑的直方图 | | geom_qq() | Q-Q 图将观察到的数据分位数与理论分位数进行对比 | | testdistr() | 自动查看密度或 Q-Q 图 | | winsorizor() | 用最接近的非异常值替换异常值 | | plot_grid() | 将多个图放入定义的网格中 | | stat_function() | 在当前图形的顶部绘制一个函数 | | fitdistr() | 获取数据、分布和分布的参数 | | logLik() | LL 是数据来自 fitdistr()的可能性，越大越好 |

二、多元数据可视化

前一章介绍了单变量数据可视化的方法。这一章延续了那个主题，但是从可视化单个变量转移到一次可视化多个变量。除了像前一章一样检查分布和异常值，我们还将讲述如何可视化变量之间的关系。可视化变量之间的关系尤其有助于更传统的统计模型，其中数据分析师必须指定函数形式(例如，线性、二次等)。).在后面的章节中，我们还将介绍机器学习模型，它采用算法来学习数据中的函数形式，而不需要分析师指定它。

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(MASS)
library(mvtnorm)
library(mgcv)
library(quantreg)
library(JWileymisc)
library(data.table)

options(width = 70, digits = 2)

2.1 分销

尽管评估单变量分布相对容易，但多变量分布更具挑战性。多个个体变量的联合分布由它们的个体分布组成。因此，除了所有可能的单变量分布，单变量分布的不同组合被组合以形成正在研究的联合分布。从可视化的角度来看，也很难在通常用于绘制数据的两个维度中可视化多个维度的数据。

在本书中，我们将涉及的唯一多元分布是多元正态(MVN)分布。在实践中，这不是过分的限制，因为大多数分析一次只关注一个结果，并且只需要知道结果的分布。因子分析和结构方程模型是两种同时模拟多种结果的常见分析类型，通常假设数据是多元正态的。如果一个人可以评估数据是否是多元正态的，我们的经验涵盖了大多数做出多元分布假设的分析。

正态分布由两个参数决定，均值 μ 和标准差 σ ，分别控制分布的位置和规模。 p 维的多元正态分布由两个矩阵控制，一个是均值的 px 1 矩阵 μ ，另一个是 pxp 协方差矩阵σ。

在二元情况下，其中 p = 2，绘制多元分布图是简单的。我们使用mvtnorm()包中的rmvnorm()函数来模拟一些多元正态数据。然后可以使用geom_density2d()函数和ggplot2绘制经验密度，如图 2-1 所示。

img/439480_1_En_2_Fig1_HTML.png

图 2-1

多元正态数据的 2D 经验密度图

mu <- c(0, 0)
sigma <- matrix(c(1, .5, .5, 1), 2)

set.seed(1234)
d <- as.data.table(rmvnorm(500, mean = mu, sigma = sigma))
setnames(d, names(d), c("x", "y"))

ggplot(d, aes(x, y)) +
  geom_point(colour = "grey60") +
  geom_density2d(size = 1, colour = "black") +
  theme_cowplot()

检查经验密度是有帮助的，但是不能将它们与多元正态分布下的预期进行精确的比较。为了将获得的经验密度与多元正态分布的预期密度进行对比，我们基于数据范围生成了成对的 x 和 y 值的网格，然后使用dmvnorm()函数计算每个 x，y 对的密度。多元正态分布的参数 μ 和σ由观测数据计算得出。结果如图 2-2 所示，为了清晰起见，这次删除了原始数据点。

testd <- as.data.table(expand.grid(
  x = seq(from = min(d$x), to = max(d$x), length.out = 50),
  y = seq(from = min(d$y), to = max(d$y), length.out = 50)))
testd[, Density := dmvnorm(cbind(x, y), mean = colMeans(d), sigma = cov(d))]

ggplot(d, aes(x, y)) +
  geom_contour(aes(x, y, z = Density), data = testd,
               colour = "blue", size = 1, linetype = 2) +
  geom_density2d(size = 1, colour = "black") +
  theme_cowplot()

在图 2-2 中，经验密度和正态密度非常接近，我们可以得出结论，假设数据是多元正态是合理的。接下来，我们模拟两个非多元正态的正态分布变量。图 2-3 显示了每个变量的单变量密度图，表明它们呈正态分布。

img/439480_1_En_2_Fig3_HTML.png

图 2-3

显示模拟变量为单变量正态的单变量密度图

img/439480_1_En_2_Fig2_HTML.png

图 2-2

2D 经验密度与多元正态密度图

set.seed(1234)
d2 <- data.table(x = rnorm(500))
d2[, y := ifelse(abs(x) > 1, x, -x)]

plot_grid(
  testdistr(d2$x, plot = FALSE)$Density,
  testdistr(d2$y, plot = FALSE, varlab = "Y")$Density,
  ncol = 2)

尽管变量是单个正态的，但这并不保证它们是多元正态的。这是一个关键点，并强调了评估多元分布的重要性，如果使用一种分析技术，对多元分布作出假设(例如，验证性因素分析)。使用我们之前使用的相同代码，我们绘制了图 2-4 中的经验和多元正态密度。

testd2 <- as.data.table(expand.grid(
  x = seq(from = min(d2$x), to = max(d2$x), length.out = 50),
  y = seq(from = min(d2$y), to = max(d2$y), length.out = 50)))
testd2[, Density := dmvnorm(cbind(x, y), mean = colMeans(d2), sigma = cov(d2))]

ggplot(d2, aes(x, y)) +
  geom_contour(aes(x, y, z = Density), data = testd2,
               colour = "blue", size = 1, linetype = 2) +
  geom_density2d(size = 1, colour = "black") +
  theme_cowplot()

图 2-4 清楚地显示了这两个变量不是多元正态的。尽管这是一个极端的例子，但它强调了单变量评估可能会产生多大的误导。

img/439480_1_En_2_Fig4_HTML.png

图 2-4

显示非多元正态数据的 2D 密度图

当 p > 2 时，直接可视化一个多元正态分布是困难的。相反，我们可以利用 Mahalanobis [60]的工作，他开发了一种方法来计算数据与“中心”的距离，数据和中心可以是多维的。距离测量被命名为 Mahalanobis 距离[60]，并将其用于我们数据中的每种情况，我们可以计算其与(多变量)中心的距离。假设数据是多元正态的，这些距离将分布为具有 p 自由度的卡方变量。

我们之前看到的testdistr()函数还包括一个基于马氏距离绘制多元正态数据的选项。我们模拟的双变量正态数据的结果如图 2-5 所示，我们模拟的单变量正态而非多变量正态数据如图 2-6 所示。

img/439480_1_En_2_Fig6_HTML.png

图 2-6

叠加多元正态分布的密度图和 Q-Q 图显示非多元正态的数据

img/439480_1_En_2_Fig5_HTML.png

图 2-5

叠加多元正态分布的密度图和显示多元正态数据的 Q-Q 图

testdistr(d, "mvnorm", ncol = 2)

testdistr(d2, "mvnorm", ncol = 2)

即使当变量的数量增加时，只要数据正在进行多元正态性检验，图表也是相似的。例如，我们可以测试mtcars数据集，它包含 32 辆汽车不同方面的 11 个变量。结果如图 2-7 所示，表明数据近似为多元正态，尽管与 Q-Q 图上直线的偏差表明了一些非正态性。

img/439480_1_En_2_Fig7_HTML.png

图 2-7

mtcars 数据的密度图叠加多元正态分布和 Q-Q 图

testdistr(mtcars, "mvnorm", ncol = 2)

2.2 异常值

对于多元数据，任何给定变量的值都可能是异常的，或者在多元空间中是异常的。如果一个值是单变量反常的，它也将是多变量反常的。然而，正如我们看到的，单变量正常变量不一定是多变量正常的，同样，单变量异常值是多变量异常的，但单变量非异常值并不保证是多变量非异常的。

在下面的代码中，我们模拟了强正相关的多元正态数据，并向V3添加了两个异常值。单独检查V3(即单变量)，一个异常值清晰可见，移除后，没有其他异常值出现(图 2-8 )。

img/439480_1_En_2_Fig8_HTML.png

图 2-8

具有异常值(图 A)和移除异常值(图 B)的数据的密度图叠加正态分布

mu <- c(0, 0, 0)
sigma <- matrix(.7, 3, 3)
diag(sigma) <- 1

set.seed(12345)
d <- as.data.table(rmvnorm(200, mean = mu, sigma = sigma))[order(V1)]
d[c(1, 200), V3 := c(2.2, 50)]

plot_grid(
  testdistr(d$V3, extremevalues = "theoretical", plot=FALSE)$Density,
  testdistr(d[V3 < 40]$V3, extremevalues = "theoretical", plot=FALSE)$Density, ncol = 2, labels = c(“A”, “B”))

接下来，我们可以再次使用testdistr函数来寻找多元异常值。使用所有情况的结果如图 2-9 所示，清楚地显示了一个异常值，没有其他点出现异常。然而，在排除一个已识别的异常情况后，第二个多元异常情况出现(图 2-10 )。当更极端的值被移除时，这种额外的异常情况被揭露，导致变量V3的均值和方差减小。

img/439480_1_En_2_Fig10_HTML.png

图 2-10

移除一个极端异常值后的多元正态和(多元)异常值的图表

img/439480_1_En_2_Fig9_HTML.png

图 2-9

多元正态和(多元)异常值的图表

testdistr(d, "mvnorm", ncol = 2, extremevalues = "theoretical")

testdistr(d[V3 < 40], "mvnorm", ncol = 2, extremevalues = "theoretical")

通过使用多元正态分布参数的稳健估计器:均值和协方差矩阵，可以直接识别多个异常情况，而不是迭代地去除异常情况。当robust选项与testdistr()的mvnorm分布结合使用时，使用robustbase包中的快速最小协方差行列式(MCD)估计器【82，83】。结果如图 2-11 所示。使用多元正态分布和稳健估计量，单变量和多变量异常情况都可以在一次通过中识别。去除这些数据后，数据看起来接近多元正态分布(图 2-12 )。

img/439480_1_En_2_Fig12_HTML.png

图 2-12

移除了两个异常值的多元正态和(多元)异常值的图表

img/439480_1_En_2_Fig11_HTML.png

图 2-11

使用稳健估计量的多元正态和(多元)异常值的图表

testdistr(d, "mvnorm", ncol = 2, robust = TRUE, extremevalues = "theoretical")

testdistr(d[-c(1,200)], "mvnorm", ncol = 2, extremevalues = "theoretical")

2.3 变量之间的关系

对于连续变量，大多数模型假定某种函数形式，变量之间通常是线性关系。有许多方法来检查这一点，但一个快速的方法是将预测值 x 切割成 k 个离散的容器。如果有有意义的断点，可以使用这些断点，但通常使用五分位数、四分位数或三分位数，这取决于数据的数量。在大型数据集中，即使是细粒度的分割，如十分位数，也可能有意义。

将预测因子分成离散的组后，可以绘制箱线图或均值图来显示趋势的大致形状。在下面的代码中，我们模拟了预测值 x 和结果 y 之间的二次关系。使用%+%操作符，我们可以重复使用同一个图，并简单地更新数据以使用具有不同分界点集合(四分位数、五分位数、十分位数)的数据。结果如图 2-13 所示。尽管在具有三等分切割的面板 A 中趋势有些不清楚，但是在具有十分位数的面板 D 中，在这少量的数据中，结果变得更加嘈杂。

img/439480_1_En_2_Fig13_HTML.png

图 2-13

将连续变量切割成四分位数的箱线图，显示非线性关系

set.seed(12345)
d2 <- data.table(x = rnorm(100))
d2[, y := rnorm(100, mean = 2 + x + 2 * xˆ2, sd = 3)]

p.cut3 <- ggplot(
  data = d2[, .(y,
    xcut = cut(x, quantile(x,
      probs = seq(0, 1, by = 1/3)), include.lowest = TRUE))],
  aes(xcut, y)) +
  geom_boxplot(width=.25) +
  theme(axis.text.x = element_text(
          angle = 45, hjust = 1, vjust = 1)) +
  xlab("")

p.cut4 <- p.cut3 %+% d2[, .(y,
    xcut = cut(x, quantile(x,
      probs = seq(0, 1, by = 1/4)), include.lowest = TRUE))]

p.cut5 <- p.cut3 %+% d2[, .(y,
    xcut = cut(x, quantile(x,
      probs = seq(0, 1, by = 1/5)), include.lowest = TRUE))]

p.cut10 <- p.cut3 %+% d2[, .(y,
    xcut = cut(x, quantile(x,
      probs = seq(0, 1, by= 1/10)), include.lowest = TRUE))]

plot_grid(
  p.cut3, p.cut4,
  p.cut5, p.cut10,
  ncol = 2,
  labels = c("A", "B", "C", "D"),
  align = "hv")

确定两个变量之间关系的函数形式的另一种方法是使用局部加权回归(黄土)线[21]。黄土线背后的思想类似于拟合最佳拟合的直线，但是加权到附近的数据点。以下代码创建一个散点图，并覆盖黄土线和直线回归直线(图 2-14 )。黄土线很容易识别被直线遗漏的二次趋势。

img/439480_1_En_2_Fig14_HTML.png

图 2-14

黄土最佳拟合线呈现非线性关系

ggplot(d2, aes(x, y)) +
  geom_point(colour="grey50") +
  stat_smooth(method = "loess", colour = "black") +
  stat_smooth(method = "lm", colour = "blue", linetype = 2)

一旦我们知道了近似的函数形式，我们就可以尝试用参数函数来近似它。在下面的代码中，我们修改了线性模型 smooth，以包含一个自定义公式，该公式指示应该在 x 和x?? 2 上回归 y 。结果得到了很大的改善，显示黄土线和二次线之间非常一致(图 2-15 )。

img/439480_1_En_2_Fig15_HTML.png

图 2-15

黄土线和二次线

ggplot(d2, aes(x, y)) +
  geom_point(colour="grey50") +
  stat_smooth(method = "loess", colour = "black") +
  stat_smooth(method = "lm",
              formula = y ~ x + I(xˆ2),
              colour = "blue", linetype = 2)

尽管黄土线有其优点，但也不是绝对可靠的。一个选择是平滑度。平滑由span参数控制，该参数传递给负责估计直线的loess()函数。在下一个示例中，绘制了两条黄土线，一条跨度低，一条跨度高。跨度越大，线条越平滑。图 2-16 为跨度为 0.2 和 2.0 的两条黄土线。

img/439480_1_En_2_Fig16_HTML.png

图 2-16

平滑程度不同的黄土线

ggplot(d2, aes(x, y)) +
  geom_point(colour="grey50") +
  stat_smooth(method = "loess", span = .2,
              colour = "black") +
  stat_smooth(method = "loess", span = 2,
              colour = "black", linetype = 2)

接下来，我们将扩展到两个变量之外，并研究许多变量的可视化方式。以下代码模拟的结果是三个预测变量(两个连续预测变量和一个分类预测变量)的函数。x和y之间关系的初始双变量图表明关系相对较弱，可能有一些异常值(图 2-17 )。

img/439480_1_En_2_Fig17_HTML.png

图 2-17

多元数据中二元关系的黄土线

set.seed(1234)
d3 <- data.table(
  x = rnorm(500),
  w = rnorm(500),
  z = rbinom(500, 1, .4))
d3[, y := rnorm(500, mean = 3 +
       ifelse(x < 0 & w < 0, -2, 0) * x +
       ifelse(x < 0, 0, 2) * w * xˆ2 + 4 * z * w,
    sd = 1)]
d3[, z := factor(z)]

ggplot(d3, aes(x, y)) +
  geom_point(colour="grey50") +
  stat_smooth(method = "loess", colour = "black")

尽管对于完全连续的变量，数据可视化通常仅限于二维，但我们可以通过图形的形状、颜色和面板来添加额外的维度。下面的代码再次检查了结果y的预测值，这次使用了所有的变量。二进制的z用于给点和线着色，我们使用前面的技巧将连续的w切割成四分之一。虽然剪切w不会完美地捕捉它的连续关系，但它足以暗示正在发生什么，并表明是否有必要与w互动。图 2-17 中提到的几个相对极端的值通过 winsorizing 底部和顶部 1%的数据来解决。最终结果如图 2-18 所示。

img/439480_1_En_2_Fig18_HTML.png

图 2-18

多元数据的黄土线

ggplot(d3, aes(x, winsorizor(y, .01), colour = z)) +
  geom_point() +
  stat_smooth(method = "loess") +
  scale_colour_manual(values = c("1" = "black", "0" = "grey40")) +
  facet_wrap(~ cut(w, quantile(w), include.lowest = TRUE))

图 2-18 显示，当w为低时，对于低于 0 的x值，在y和x之间往往有一个负的、大致线性的斜率，但是当w为高时，没有关系。根据w的电平，z会向下或向上移动数值。

除黄土外，广义加性模型(gam)是另一种拟合直线的灵活方法[40]。我们将在后面讨论它们的统计性质和用途。现在的优势在于拟合灵活的非线性平滑项，包括连续变量和分类变量之间的相互作用。具体来说，我们利用了mgcv包【122】中的gam()函数。

在黄土中，我们使用span参数指定平滑程度。在 GAMs 中，我们通过k参数使用近似自由度来指定允许生产线有多灵活。目前，我们并不关注 gam 的统计属性或正式的统计推断，而是将它们作为一种工具来拟合数据并生成预测，以图形化和可视化我们数据中的(平滑)模式。te()功能用于允许x和w之间的非线性交互，并且我们还允许这些对于不同级别的z是分开的。拟合的模型存储在对象m中。接下来，我们生成一个用于预测的网格值。与密度一样，我们在x和w上为z的所有级别选择等距点。生成预测后，我们可以绘制如图 2-19 所示的结果。注意图 2-19 颜色最好看。

对于图 2-19 中的一条轮廓线，所有预测的y值都是相同的。因此，通过追踪一条线，你可以检查预测因素与结果的关系。例如，当z为 0 时(图 2-19 的左图)，对于低于 0 的x值，w是否在-3 和 0 之间移动对预测的y值影响不大。一般来说，任何平行于特定维度的等高线都表示在该维度上该点的预测变化很小。光栅背景也有助于显示预测的y值。

img/439480_1_En_2_Fig19_HTML.png

图 2-19

显示 x 和 w 的不同组合的预测 y 值的等值线图，由 z 面板显示

m <- gam(winsorizor(y, .01) ~ z + te(x, w, k = 7, by = z), data = d3)

newdat <- expand.grid(
  x = seq(min(d3$x), max(d3$x), length.out = 100),
  w = seq(min(d3$w), max(d3$w), length.out = 100),
  z = factor(0:1, levels = levels(d3$z)))

newdat$yhat <- predict(m, newdata = newdat)

ggplot(newdat, aes(x = x, y = w, z = yhat)) +
  geom_raster(aes(fill = yhat)) +
  geom_contour(colour = "white", binwidth = 1, alpha = .5) +
  facet_wrap(~ z)

评估方差的同质性

方差的同质性或同质性是指各组之间或连续预测值之间存在相同的有限方差，结果在预测值的各个级别上具有相同的残差方差。例如，在单向方差分析中，对于解释变量(通常称为自变量)的每个水平，结果的方差应该大致相等。

为了评估方差的同质性，我们为 iris 数据制作了一个数据表，并通过简单地计算物种的方差来进行第一次检查。

diris <- as.data.table(iris)
diris[, .(V = var(Sepal.Length)), by = Species]

##       Species    V
## 1:     setosa 0.12
## 2: versicolor 0.27
## 3:  virginica 0.40

可视化数据和分布也很有帮助。箱线图或盒须图可能是一种有用的方法，因为“盒”部分覆盖了四分位数之间的范围，即第 25 至 75 个百分点。这是一种稳健的分布范围测量方法。我们检查盒子的分布,以了解不同物种之间的可变性是否大致相等。箱线图也比仅仅按组计算方差更能提供信息，因为它们可以显示数据中是否有任何异常值。

即使物种间的传播是可比的，中位数或位置也可能不同，这使得比较传播更加困难。如果您不想检查位置，并且希望只关注分布，那么在制图之前确定中间值是很有用的。图 2-20 中显示了无(A 图)和有(B 图)中间居中的箱线图。

img/439480_1_En_2_Fig20_HTML.png

图 2-20

各种萼片长度的盒须图。异常值显示为点

plot_grid(
  ggplot(diris, aes(Species, Sepal.Length)) +
    geom_boxplot() +
    xlab(""),
  ggplot(diris[, .(Sepal.Length = Sepal.Length -
                              median(Sepal.Length)), by = Species],
         aes(Species, Sepal.Length)) +
    geom_boxplot() +
    xlab(""),
  ncol = 2, labels = c("A", "B"), align = "hv")

我们已经看到了密度图，并用它来评估变量的分布；如果密度图被反射形成镜像，则称为小提琴图。我们可以用geom_violin()来制作；它提供了与箱线图相似的信息。然而，violin 图提供了关于分布的更多信息，因为箱线图仅显示了中位数(第 50 个百分点)、第 25 个百分点和第 75 个百分点，以及数据的范围(如果有异常值，则略小于整个范围)。但是，由于查看中位数和四分位间距范围也很有用，因此将宽度较窄的箱线图叠加到 violin 图上会很有帮助。从小提琴图中，我们可以看到，对于具有较高分布均值(位置)的物种，图 2-21 中的扩散稳步增加。

img/439480_1_En_2_Fig21_HTML.png

图 2-21

中间有盒须图的小提琴图

ggplot(diris, aes(Species, Sepal.Length)) +
  geom_violin() +
  geom_boxplot(width = .1) +
  xlab("")

正如我们在研究变量之间的关系时所做的那样，violin 图和 box 图可以扩展到使用颜色和面板来处理多个变量。下面的代码使用我们研究变量之间关系的数据，并删除所有连续的预测因子，以便于检查预测因子水平上结果的分布和传播。结果如图 2-22 所示。请注意，我们有目的地限制了 y 轴的范围，以便更容易看到图形，并且不太强调极值。

img/439480_1_En_2_Fig22_HTML.png

图 2-22

小提琴图，中间用 x 的四分位数表示盒须图，用 z 着色

## create cuts
d3[, xquartile := cut(x, quantile(x), include.lowest = TRUE)]
d3[, wquartile := cut(w, quantile(w), include.lowest = TRUE)]
d3[, yclean := winsorizor(y, .01)]

## median center y by group to facilitate comparison
d3[, yclean := yclean - median(yclean),
   by = .(xquartile, wquartile, z)]

p <- position_dodge(.5)

ggplot(d3, aes(xquartile, yclean, colour = z)) +
  geom_violin(position = p) +
  geom_boxplot(position = p, width = .1) +
  scale_colour_manual(values = c("1" = "black", "0" = "grey40")) +
  facet_wrap(~ wquartile) +
  theme(axis.text.x = element_text(angle = 45, hjust=1, vjust=1)) +
  coord_cartesian(ylim = c(-5, 5), expand = FALSE)

到目前为止，我们已经研究了评估连续结果和离散解释变量的方差齐性，或者将连续解释变量分成离散类别。现在我们转向如何对连续变量进行可视化处理。第一步是创建一个散点图；然而，为了看到“传播”,我们需要对方差或四分位范围等进行一些连续的估计。这可以通过分位数回归来实现。我们不会在这里详细讨论这个过程，但可以说分位数回归可以估计作为一个或多个解释变量的函数的分布的分位数。我们将首先模拟一些同质和异质数据。

set.seed(1234)
d4 <- data.table(x = runif(500, 0, 5))
d4[, y1 := rnorm(500, mean = 2 + x, sd = 1)]
d4[, y2 := rnorm(500, mean = 2 + x, sd = .25 + x)]

图 A 是同质可变性的例子，图 B 是异质可变性的例子。我们在图 2-23 中展示了这两个数据集的视觉对比。分位数回归[51，50，107]。

img/439480_1_En_2_Fig23_HTML.png

图 2-23

同异方差与异方差

plot_grid(
  ggplot(d4, aes(x, y1)) +
    geom_point(colour = "grey70") +
    geom_quantile(quantiles = .5, colour = 'black') +
    geom_quantile(quantiles = c(.25, .75),
                  colour = 'blue', linetype = 2) +
    geom_quantile(quantiles = c(.05, .95),
                  colour = 'black', linetype = 3),
  ggplot(d4, aes(x, y2)) +
    geom_point(colour = "grey70") +
    geom_quantile(quantiles = .5, colour = 'black') +
    geom_quantile(quantiles = c(.25, .75),
                  colour = 'blue', linetype = 2) +
    geom_quantile(quantiles = c(.05, .95),
                  colour = 'black', linetype = 3),
  ncol = 2, labels = c("A", "B"))

## Smoothing formula not specified. Using: y ~ x
## Smoothing formula not specified. Using: y ~ x
## Smoothing formula not specified. Using: y ~ x
## Smoothing formula not specified. Using: y ~ x
## Smoothing formula not specified. Using: y ~ x
## Smoothing formula not specified. Using: y ~ x

2.4 总结

本章演示了可视化多元数据的技术，以了解它如何与多元正态分布进行比较(表 2-1 )。

表 2-1

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

功能

它的作用

| | --- | --- | | geom_quantile() | 以给定的分位数拟合直线(默认为四分位数) | | geom_density2d() | 为多元正态数据点创建 2D 密度等值线图 | | geom_contour() | 创建 2D 等高线图 | | geom_violin() | 镜像密度图 | | stat_smooth() | 将曲线拟合到数据点 |

三、GLM 1

广义线性模型(GLMs)是一大类模型，包括回归分析和方差分析(ANOVA ),是常用于指 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(data.table)
library(ggplot2)
library(visreg)
library(ez)
library(emmeans)
library(rms)
library(ipw)
library(JWileymisc)
library(RcppEigen)
library(texreg)

options(
  width = 70,
  stringsAsFactors = FALSE,
  datatable.print.nrows = 20,
  datatable.print.topn = 3,
  digits = 2)

3.1 概念背景

广义线性模型提供了一个通用的框架和符号，可以适应许多特定类型的模型和分析。本章涵盖了几种特定类型的 GLMs。如果您不熟悉一些更基本的 GLMs 类型，如方差分析(ANOVA)或线性回归，并且您希望更好地了解 GLMs 的背景和概念框架，那么阅读 ANOVA 和线性回归可能是值得的。一个易于使用的免费资源是在线统计教育:一门多媒体学习课程( http://onlinestatbook.com/ ，项目负责人:莱斯大学的大卫·m·莱恩)，其中有一节是关于方差分析和回归的。关于 GLMs 的全面数理统计背景，我们推荐 McCullagh 和 Nelder 的经典著作[61]，其完整参考资料在参考资料部分。

首先，我们介绍一些统计和 GLMs 中常用的函数，包括:

E ( x )表示一个变量的期望值，或者它的(可能有条件的)均值。
Var ( x )表示一个变量的方差，或者说它的离差。
g ( x )表示链接函数，该函数获取原始结果变量并将其转换为由 GLM 预测的线性标度。
g^-1(x)表示反向链接函数，它采用 GLM 预测的线性标度，并将其转换回原始标度。
exp ( x )表示指数函数。
ln ( x )表示自然对数函数。

glm 的结构是这样的，有一个结果或因变量，我们称之为 y。一个或多个预测值或解释变量存储在一个矩阵中，该矩阵有 n 行(数据点的数量)和 k 列(预测值或解释变量的数量)，我们称之为 x。回归系数，即要估计的参数，形成一个长度为 k (预测值的数量)的向量，称为 β 。结局的期望值是 E (y) = μ 。加粗的希腊文小写字母 μ 是书写预期结果值的较短方式。这些将总是在原始或原始数据范围内。线性标尺上的期望值是另一个称为 η 的向量。

有了这些约定，我们就可以像内尔德和威德伯恩的开创性文章[71]那样定义广义线性模型的构建模块。我们使用与内尔德和威德伯恩[71]稍有不同的符号来反映应用统计学中更常见的近期实践。每个 GLM 都有一个结果或因变量， y 。在线性预测因子上进行调节后，假设结果遵循来自指数族的概率分布。此外，每个 GLM 具有一组 k 预测变量，x ₁ ，...，x _k ，总体称为 X，以及一个期望的线性结果:

$\eta = X\beta$

(3.1)

最后，每个 GLM 都有一个链接功能

$\eta =g\left(\mu \right)=g\left(E(y)\right)$

(3.2)

以及类似的反向链接功能

$E(y)=\mu ={g}^{-1}\left(\eta \right)$

(3.3)

因为参数 β 总是在线性预测标度 η 上估计，所以不管结果的分布如何，GLMs 的估计都是相似的。改变的是为结果和链接以及反向链接函数假设的分布。我们涵盖了三种常见类型的结果变量的分布和链接函数(以及如何分析它们):连续、二元结果和计数结果。表 3-1 显示了我们在本章中涉及的最常用的(规范的)链接函数。请注意，列出的分布并不是给定结果类型的所有可能分布。例如，对于连续结果，根据其形状或界限，有许多其他分布(例如，连续但以 0 和 1 为界限的数据的 Beta 分布)。

表 3-1

结果类型、分布和相应的链接函数

结果类型

分布

链接功能

反向链接功能

| | --- | --- | --- | --- | | 连续的(实数) | 正常(高斯) | =【g】()=** | =【g】^-1()=** | | 二进制(0/1) | 伯努利多项式 | $\eta =g\left(\mu \right)=\mathit{\ln}\left(\frac{\mu }{1-\mu}\right)$ | $\mu ={g}^{-1}\Big(\eta =\frac{1}{1+\mathit{\exp}\left(-\eta \right)}$ | | 计数(正整数) | 泊松，负二项式 | =【g】()=(** | =【g】^-1()=【exp】** |

基于结果分布，有可能写出一个似然函数。我们在本书中不涉及似然函数的细节，因为没有必要知道什么是应用数据分析的似然函数。然而，基本思想是，似然函数(取决于所选择的分布)量化了数据从具有给定参数集的分布中出现的可能性；因此，一般似然函数被写成 L (y， θ )。注意，这里我们将参数称为 θ ，因为对于许多(但不是所有)分布，必须估计回归系数之外的附加参数；具体来说，经常需要估计色散参数。例如，正态分布的参数是均值和方差或标准差。另一个注意事项是，出于实用的原因，可能性通常被报告为对数可能性。其他地方描述了进一步的细节[71，61]。

似然函数有几个有用的目的。首先，它们是最大似然估计的基础。也就是说，估计 GLM(以及许多其他模型)的参数，以最大化数据的似然函数。此外，来自给定模型的总体似然性可以与来自另一个模型的似然性进行比较，作为哪个模型提供更好的数据拟合的相对比较。接下来，我们研究如何将这个通用模型应用到具体的案例中。

3.2 分类预测器和虚拟编码

两级分类预测器

在 GLMs 中检查分类预测因子是很常见的。例如，人们可能希望测试结果中是否存在性别差异。但是，在这样做之前，需要一个系统将“女人”和“男人”转换成数值。将离散类别编码成数字的最常见系统称为虚拟编码。虚拟编码包括对代表一个特定类别的一系列二进制 0/1 变量进行编码。对于性别，我们可以编写两个虚拟变量，一个代表女性，一个代表男性。如表 3-2 所示。

表 3-2

虚拟编码性别示例

性

| | --- | --- | --- | | 妇女 | one | Zero | | 男人 | Zero | one | | 男人 | Zero | one | | 妇女 | one | Zero | | 妇女 | one | Zero | | 男人 | Zero | one |

在这种情况下，伪代码 D1 和 D2 是完全负相关的，表明它们捕获相同的信息，但是被颠倒。因此，我们将只把两个伪码中的一个输入实际的 GLM。总的原则是，对于一个 k- 级变量，可以生成 k 个伪代码变量，但只包括伪代码变量的k1。省略的变量成为参考组。为了理解为什么会这样，我们可以检查一个简单的设计矩阵 X，(表 3-3 )，它对应于一个只测试性别差异的模型。作为 GLMs 的标准，我们有一个常数列，它是截距，当所有其他预测值都等于零时结果的期望值。我们还有两个虚拟编码变量之一，D2。在这种情况下，当参与者是男性时，D2 为 1，当参与者是女性时，为 0。相应的系数向量将具有两个元素，一个用于截距，一个用于伪码 D2。

表 3-3

具有截距和一个虚拟变量的性别差异设计矩阵

拦截

| | --- | --- | | one | Zero | | one | one | | one | one | | one | Zero | | one | Zero | | one | one |

截距系数将是女性的期望值，因为只有当参与者是女性时，所有其他变量(在这种情况下，只有 D2)才为零。D2 系数是 D2 一个单位变化结果的预期变化。因为我们用 0/1 来编码 D2，在 D2 一个单位的变化正好是从女性(0)到男性(1)的转变。因此，D2 系数可以更简单地认为是男女之间的预期差异。

三层或更多层分类预测因子

具有三个或更多级别的虚拟编码分类变量的工作方式类似于基本的两级方法。同样，我们创建一组 0/1 变量，对变量的每个特定级别进行编码。表 3-4 显示了锻炼类型的三级变量示例。

表 3-4

虚拟编码性别示例

锻炼

| | --- | --- | --- | --- | | 奔跑 | one | Zero | Zero | | 游泳 | Zero | one | Zero | | 自行车 | Zero | Zero | one | | 奔跑 | one | Zero | Zero | | 自行车 | Zero | Zero | one | | 自行车 | Zero | Zero | one | | 奔跑 | one | Zero | Zero | | 游泳 | Zero | one | Zero | | 游泳 | Zero | one | Zero |

为了进行分析，我们将在 GLM 中包括任何k1 = 31 = 2 的伪码变量，同样，被排除的组将是参考组。然而，组间比较变得稍微复杂一些。假设我们忽略了 D1，运行的虚拟代码。游泳和自行车虚拟代码的系数将捕捉跑步和游泳之间以及跑步和自行车之间的预期差异。然而，没有一个系数会直接测试游泳和骑自行车之间的差异。总的原则是，对于一个k水平变量，有 $\left(\begin{array}{l}k\\ {}2\end{array}\right)$ 可能的成对比较；对于三级变量，三次成对比较；对于四级变量，六次成对比较。要测试所有成对组合，有几种选择。一种选择是，尽管效率不高，但只要有伪代码变量，就运行 GLM，每次迭代只省略一个伪代码。另一种方法是使用从模型估计的系数矩阵和参数协方差矩阵来测试特定的对比。在我们的练习示例中，测试游泳与跑步的系数是否不同于自行车与跑步的系数将提供游泳与自行车是否不同的测试。

考虑如何测试变量的整体效果也很重要。对于单个连续变量或两级分类变量，单个系数反映了变量的整体效应。然而，在输入多个虚拟代码变量的三级或更多级分类变量的情况下，任何特定虚拟代码变量系数的测试都不能提供该变量显著性的整体测试。相反，我们需要一个综合测试。对于具有正态分布结果的 GLMs，这可能是一个综合 f 检验。对于无法计算自由度的其他 GLM，标准测试是基于测试所有虚拟变量系数共同为零的 Wald 测试，或似然比测试，这需要在没有问题变量的情况下重新调整 GLM，并测试最终对数似然变化的程度，最终对数似然变化将作为自由度等于排除的虚拟变量数量的 χ ² 分布。

3.3 相互作用和调节效应

当两个或多个解释变量之间的关系和结果依赖于彼此的值时，就会发生交互作用。例如，Wiley 和他的同事们[115]发现了压力源的可控程度和人们积极感受的应对方式之间的相互作用。对于人们几乎无法控制的压力源，无论人们试图积极思考如何改善问题还是避免思考问题，对他们的积极情绪水平没有影响。然而，对于可控的压力源，那些避免思考如何解决问题的人积极情绪水平较低，而那些积极思考如何解决问题的人积极情绪水平较高。

两个变量之间的相互作用称为双向相互作用。三个变量之间的相互作用称为三向相互作用，等等。交互作用可以在 GLMs 中进行测试，方法是将单个变量单独添加到模型中，并添加需要交互作用的变量的乘积项。这可以通过在原始数据集中创建新变量来实现。例如，对于 x ₁ 和 x ₂ ,“新”交互变量将被创建为int = x₁⋅x₂。R和其他统计软件包还提供了指定应该在模型中测试两个(或更多)变量的交互作用的能力，在这种情况下，产品术语是自动创建的，而无需在数据集中创建额外的变量。无论哪种情况，最终结果都是一样的，表 3-5 中显示了一个设计矩阵示例。

对于两个以上变量的相互作用，还有几个附加项。如果有三个解释变量( x ₁ ， x ₂ ， x ₃ )，则有三个双向交互(x₁⋅x₂，x₁⋅x) x2⋅x3 以及除了三个个体变量之外可以考虑的一个三向交互(x1⋅x2⋅x3)。即使主要关注的是三方交互，标准做法也是包括所有低阶项。因此，对于三向互动，通常所有可能的双向互动和所有单个变量也将包括在分析中。

因为相互作用涉及不止一个变量，所以总是有可能用不止一种方式来解释它们。例如，如果两个变量 a 和 b 之间存在相互作用，则每个变量与结果 y 之间的关系取决于另一个变量的水平。因此 a 和 y 之间的关系取决于 b 的级别，同样 b 和 y 之间的关系取决于 a 。这一事实改变了对单个变量回归系数的解释。查看表 3-5 ，当 x2 = 0 时，x1 的系数将被解释为 x1 一个单位变化的预期 y 变化。同样，当 x1 = 0 时，x2 的系数将被解释为 x2 变化一个单位时 y 的预期变化。最后，x1x2 的系数可以解释为(1)x2 变化一个单位时 x1 系数的预期变化，或者(2)x1 变化一个单位时 x2 系数的预期变化。

表 3-5

双向互动设计矩阵示例

拦截

x1x2

| | --- | --- | --- | --- | | one | one | Two | Two | | one | Two | Two | four | | one | three | one | three | | one | three | three | nine | | one | one | one | one | | one | Two | Zero | Zero |

当用标准代数而不是矩阵代数写出时，这种解释更有意义，如下所示:

${y}_i={b}_0+{b}_1\cdot {x}_{1i}+{b}_2\cdot {x}_{2i}+{b}_3\cdot \left({x}_{1i}\cdot {x}_{2i}\right)$

(3.4)

这可以分解如下，以强调相互作用如何最终导致 x1 和 x2 之间的关系，其中 y 取决于相互作用中的其他变量:

$\begin{array}{l}\kern2.78em {y}_i={b}_0+\\ {}\left({b}_1+{b}_3\cdot {x}_{2i}\right)\cdot {x}_{1i}+\\ {}\kern0.5em \left({b}_2+{b}_3\cdot {x}_{1i}\right)\cdot {x}_{2i}\end{array}$

(3.5)

类似的逻辑也适用于三方互动，只是依赖于另外两个变量。因为所有较低阶的双向交互也是标准的，所以模型的复杂性(参数的数量)急剧增加。

$\begin{array}{l}{y}_i={b}_0\\ {}\kern0.75em +{b}_1\cdot {x}_{1i}+{b}_2\cdot {x}_{2i}+{b}_3\cdot {x}_{3i}\\ {}\kern0.75em +{b}_4\cdot \left({x}_{1i}\cdot {x}_{2i}\right)\\ {}\kern0.75em +{b}_5\cdot \left({x}_{1i}\cdot {x}_{3i}\right)\\ {}\kern0.75em +{b}_6\cdot \left({x}_{2i}\cdot {x}_{3i}\right)\\ {}\kern0.75em +{b}_7\cdot \left({x}_{1i}\cdot {x}_{2i}\cdot {x}_{3i}\right)\end{array}$

(3.6)

这可以分解如下，以突出每个变量对其他两个变量的依赖性:

$\begin{array}{l}{y}_i={b}_0\\ {}\kern0.75em +\left({b}_1+{b}_4\cdot {x}_{2i}+{b}_5\cdot {x}_{3i}+{b}_7\cdot \left({x}_{2i}\cdot {x}_{3i}\right)\right)\cdot {x}_{1i}\\ {}\kern0.75em +\left({b}_2+{b}_4\cdot {x}_{1i}+{b}_6\cdot {x}_{3i}+{b}_7\cdot \left({x}_{1i}\cdot {x}_{3i}\right)\right)\cdot {x}_{2i}\\ {}\kern0.75em +\left({b}_3+{b}_5\cdot {x}_{1i}+{b}_6\cdot {x}_{3i}+{b}_7\cdot \left({x}_{1i}\cdot {x}_{2i}\right)\right)\cdot {x}_{3i}\end{array}$

(3.7)

3.4 公式界面

在R中，许多模型和几乎所有 glm 都是使用公式接口指定的。公式是一种指定简单到复杂模型的灵活方式，它由两部分组成，中间用波浪号(∩)隔开。左手边(LHS)位于波浪号的左侧，右手边(RHS)位于波浪号的右侧。基本形式是

outcome ∼predictor1 + predictor2.

“+”操作符将变量添加到模型中。R的公式接口是一种灵活的指定模型的方式。主要操作符有“+”、-、:、*，它们分别是加法项、-和乘法项。使用update()功能可以修改现有公式。

除了从数据中检查变量的个别影响，GLMs 通常包括两个(或更多)变量的乘积。这是一个非常常见的任务，以至于 formula 接口有一种特殊的方式来指示应该包含两个项的乘积，即“:”运算符。例如，y∾x1 + x2 + x1:x2包括x1、x2，以及它们的交互作用(乘积项)作为y的预测因子。其中一个很好的特性是它正确地处理了连续变量和分类变量。如果x1和x2是连续测度，那么x1:x2将是正则代数积。如果其中一个或两个都是虚拟代码，那么产品将针对虚拟代码进行适当扩展。

当包含交互项时，主要效应和每个变量的单独效应几乎总是包含在内。因为个体效应几乎总是包含在交互作用中，所以可以使用“*”运算符来表示两个变量的交互作用和个体效应:y∩x1 * x2扩展为y∩x1 + x2 + x1:x2。多个操作符可以链接在一起，以便y∩x1 * x2 * x3扩展为y∩x1 + x2 + x3 + x1:x2 + x1:x3 + x2:x3 + x1:x2:x3。有时，一个变量可能会调节与三个或更多其他预测因素的相互作用，但三向或四向相互作用是不可取的。括号可用于对术语集进行分组，以便将运算符分配给组中的所有术语。于是，y∾x1 * (x2 + x3)展开为y∾x1 + x2 + x3 + x1:x2 + x1:x3。

这涵盖了公式中最常用的运算符。另外两个细节有时会有所帮助，尤其是在修改现有公式时。一个点，“.”，可以作为指代一切事物的简称。最后，可以使用“-”操作符删除术语。当使用update()函数更新现有的、通常存储的公式时，这些是最常用的。update()函数将一个现有的公式对象作为其第一个参数，然后是所需的修改。以下代码显示了可能的不同类型的公式更新示例。注意，其他运算符“*”和“:”也可以与“.”一起使用。如果点被完全省略，那么旧公式的这一部分根本不会被重用。我们将在R中使用 formula 接口来构建大部分模型，因此值得花时间彻底学习。

f1 <- y ~ x1 + x2 + x1:x2

update(f1, . ~ .)

## y ~ x1 + x2 + x1:x2

update(f1, w ~ .)

## w ~ x1 + x2 + x1:x2

update(f1, . ~ . + x3)

## y ~ x1 + x2 + x3 + x1:x2

update(f1, . ~ . - x1:x2)

## y ~ x1 + x2

3.5 差异分析

概念背景

方差分析(ANOVA)是一种统计技术，用于划分不同因素导致的结果差异。ANOVAs 是 GLM 的一个特例，具有连续、正态分布的结果和离散/分类解释变量，如性别(女性、男性)或随机实验研究中的状况(如治疗 A、治疗 B 或对照)。由于这些限制，ANOVAs 可以被概念化为测试结果的平均值是否在每个组中是相等的。也就是说，ANOVAs 测试是否:

${\mu}_{TreatmentA}={\mu}_{TreatmentB}={\mu}_{Control}$

(3.8)

传统上，ANOVAs 被用作零假设统计检验(NHST)的一部分。从本质上讲，NHST 设立了零假设，并问道，假设零假设在总体中为真，在这个数据样本中获得观察结果的概率是多少？零假设的反面是替代假设。在方差分析的情况下，另一个假设是至少一个组的平均值不等于其余的平均值。例如，治疗 A 可能比治疗 B 或对照组具有更低或更高的平均值。

为了将方差分析参数化为 GLM，我们将前面的等式写成一系列的差，例如:

$\begin{array}{l}{\mu}_{Control}-{\mu}_{TreatmentA}\\ {}{\mu}_{Control}-{\mu}_{TreatmentB}\end{array}$

(3.9)

这些差异被编码到保存预测变量 x 的设计矩阵中。如我们之前所述，使用虚拟编码将组转换为预测变量。

默认情况下，R将使用因子的第一级作为参考组来创建虚拟代码。参考组被省略，但是设计矩阵增加了截距，截距是一个包含 1 的常数列。然后，GLM 将估计参数， β ，设计矩阵的每一列一个参数。因为截距在模型中，所以治疗 A 的系数(R标记为ConditionA)将是治疗 A 的平均值和对照组(参考水平)之间的差值。同样，治疗 B 的系数将是治疗 B 的平均值和对照组之间的差值。截距将获取对照组的平均值。

为了看到这一点，我们可以使用函数lm()估计R中的回归参数，并使用函数coef()提取系数。我们用一个公式接口:outcome predictor写出我们希望 R 拟合的模型。

set.seed(1234)
example <- data.table(
  y = rnorm(9),
  Condition = factor(rep(c("A", "B", "Control"), each = 3),
                     levels = c("Control", "A", "B")))

coef(lm(y ~ Condition, data = example))

## (Intercept)  ConditionA  ConditionB
##      -0.562       0.614       0.092

通过计算每组的平均值，很容易检查这些是否与组平均值的差异相匹配。

example[, .(M = mean(y)), by = Condition]

##    Condition      M
## 1:         A  0.052
## 2:         B -0.470
## 3:   Control -0.562

我们可以立即看到截距与对照组的平均值相同。处理 A 的系数等于对照组和处理 A 组之间的差值。平均值:0.61 = 0.05 - (-0.56).

如果我们要抑制截距，那么设计矩阵将包含每个条件的虚拟代码，回归系数将正好是组平均值，我们在下面的代码中可以看到。通过在公式中添加 0 来抑制截距。但是请注意，这通常只有在模型中包含伪代码时才有意义。在具有连续解释变量的 GLM 中，抑制截距会迫使截距精确为零，这很少是明智的。

model.matrix(~ 0 + Condition, data = example)

##   ConditionControl ConditionA ConditionB
## 1                0          1          0
## 2                0          1          0
## 3                0          1          0
## 4                0          0          1
## 5                0          0          1
## 6                0          0          1
## 7                1          0          0
## 8                1          0          0
## 9                1          0          0
## attr(,"assign")
## [1] 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$Condition
## [1] "contr.treatment"

coef(lm(y ~ 0 + Condition, data = example))

## ConditionControl       ConditionA       ConditionB
##           -0.562            0.052           -0.470

标准方差分析检验了总体结果中有多少可变性，其中有多少可变性在组均值之间，以及在考虑组均值后还有多少可变性。为了检验是否有任何差异，组均值之间的变异量与组内变异量之比(剩余方差)形成。该比率称为 F 比率，可用于获得 p 值，因为 F 分布的比例比观察到的 F 比率更极端。

F 比率基于均方比率，即平方和(SS)除以自由度(DF)。总是有两个，一个是分子(利益的影响)，一个是分母(误差或残差)，具体来说

$F=\frac{M{S}_{model}}{M{S}_{residual}}=\frac{S{S}_{model}/D{F}_{model}}{S{S}_{residual}/D{F}_{residual}}$

(3.10)

自由度也用于使用pf()功能查找 F 分布的 F 比率。

pf(.72, df1 = 1, df2 = 6, lower.tail = FALSE)

## [1] 0.43

为了使 ANOVAs 中的统计检验有效，必须满足几个假设。第一，观测必须是独立的(独立性)。例如，如果对几个参与者进行重复测量，使得观察结果集中在一个人内(或在一所学校内，或任何其他分组单位内)，则违背了这一假设。第二，以解释变量为条件，结果必须是连续的和正态分布的(正态)。第三，所有解释变量的每个水平内的方差必须相等(方差齐性)。最后一个假设是必需的，因为单个剩余方差用于估计所有组的不确定性。

前面的示例在模型中使用了一个分组因子。可以同时测试多个独立变量或分组变量。包含多个解释变量也允许测试交互作用的效果——一个变量的效果是否依赖于另一个变量的水平。在某种程度上，这是对单个解释变量的巨大进步，因为变量的数量增加了一倍，增加了一个额外的交互(调节)项。然而，从 GLM 的角度来看，这是一个很小的变化。设计矩阵增加了额外的列，其中一些是伪代码的产物，而不是来自单个变量的伪代码。这些都是比较表面的区别。下面的代码修改了我们之前看到的 mtcars 数据集，将一些变量转换为因子，并为两个变量的主要影响及其相互作用创建了一个设计矩阵(如果两个变量都是连续的，则只是两个变量的乘积，如果它们是离散的，则是它们的伪码的乘积)。只有当其他两个伪码都为 1 时，交互作用(标为 vs1:am1)才为 1，它反映了当 vs = 1 和 am = 1 时，均值的差异有多大，超出了它们的平均效应预期。

mtcars <- as.data.table(mtcars)
mtcars[, ID := factor(1:.N)]
mtcars[, vs := factor(vs)]
mtcars[, am := factor(am)]

head(model.matrix(~ vs * am, data = mtcars))

##   (Intercept) vs1 am1 vs1:am1
## 1           1   0   1       0
## 2           1   0   1       0
## 3           1   1   1       1
## 4           1   1   0       0
## 5           1   0   0       0
## 6           1   1   0       0

R 中的方差分析

注意

方差分析检验正态分布结果的组均值是相等还是不同。ez包中的ezANOVA()函数可以运行独立测量、重复测量和混合模型方差分析，并提供假设测试。

为了在R中运行 ANOVA，我们使用了ez包中的ezANOVA()函数，该函数适合各种类型的 ANOVA，并提供 ANOVA 通常报告的附加信息。要使用它，我们需要为ezANOVA()添加一个 ID 变量。ezANOVA()函数接受一个数据集、结果变量(dv)、主题 ID 变量(wid)、主题变量之间的变量(between)。其余参数是可选的，它们控制组不平衡时方差的计算方式以及要打印的输出量。在接下来的代码中，我们使用我们的小型示例数据集，测试条件的整体效果。它还打印一个统计测试，无论方差齐性假设是否满足，Levene 的测试。小的 F 比率(高 p 值)表明很少有证据表明条件之间的差异显著。

example[, ID := factor(1:.N)]

print(ezANOVA(
  data = example,
  dv = y,
  wid = ID,
  between = Condition,
  type = 3,
  detailed = TRUE))

## Coefficient covariances computed by hccm()

## $ANOVA
##        Effect DFn DFd  SSn SSd    F    p p<.05   ges
## 1 (Intercept)   1   6 0.96   8 0.72 0.43       0.108
## 2   Condition   2   6 0.66   8 0.25 0.79       0.076
##
## $'Levene's Test for Homogeneity of Variance'
##   DFn DFd SSn SSd    F    p p<.05
## 1   2   6 1.5 6.1 0.73 0.52

接下来，我们使用包含交互作用的 mtcars 数据集进行方差分析。我们在 subjects 变量和 interaction 变量之间添加了一个加法，但是其余的保持不变。我们首先看到的是关于每种情况下样本大小不等的警告。当样本大小在不同组之间不平衡时，不同的平方和计算方法会产生不同的结果，这是一个有争议的问题[55]。这些结果表明相互作用项没有显著影响。vs 和 am 的两个主要效应(即其本身)在统计上是显著的，并且具有大的效应大小，表明它们与每加仑英里数有关。同样，没有证据表明违反了方差齐性假设。与 GLM 框架相比，F 比率测试对回归系数的类似影响，但方式略有不同。当一个因子有两个以上的级别时，F 比率有多个分子自由度，相当于测试多个回归系数是否同时为零，而不是一次测试一个系数。以下代码输出的最后一列ges显示了广义 eta 平方的效果大小度量。然而，除了测试是如何构建的，ANOVAs 使用的潜在线性模型是 GLMs 允许的模型的子集。

print(ezANOVA(
  data = mtcars,
  dv = mpg,
  wid = ID,
  between = vs * am,
  type = 3,
  detailed = TRUE))

## Warning: Data is unbalanced (unequal N per group). Make sure you specified a well-considered value for the type argument to ezANOVA().

## Coefficient covariances computed by hccm()

## $ANOVA
##        Effect DFn DFd   SSn SSd      F       p p<.05   ges
## 1 (Intercept)   1  28 13144 337 1090.6 5.7e-24     * 0.975
## 2          vs   1  28   382 337   31.7 4.9e-06     * 0.531
## 3          am   1  28   284 337   23.5 4.2e-05     * 0.457
## 4       vs:am   1  28    16 337    1.3 2.6e-01       0.045
##
## $'Levene's Test for Homogeneity of Variance'
##   DFn DFd SSn  SSd    F    p p<.05
## 1   3  28  15 156 0.88 0.46

我们可以使用 Tukey 的诚实显著性差异(HSD)来测试细胞之间的成对差异。下面的代码创建一个新变量，它是 vs 和 am 的组合，然后创建一个显示平均值和 95%置信区间(解释为以相同方式进行的 95%的区间将包括真实总体参数)的图形。任何共享一个字母的细胞在统计学上彼此没有显著差异。不共享一个字母的细胞在统计学上有显著差异。我们使用JWileymisc包中的TukeyHSDgg()进行绘图，调整轴标签的角度并删除 x 轴标题。图 3-1 显示了结果。

img/439480_1_En_3_Fig1_HTML.png

图 3-1

具有置信区间的单元均值图。基于 Tukey 的诚实显著差异，共享字母的细胞在统计上没有显著差异。

mtcars[, Cells := factor(sprintf("vs=%s, am=%s", vs, am))]
TukeyHSDgg("Cells", "hp", mtcars) +
  theme(axis.text.x = element_text(angle=45, hjust=1, vjust=1)) +
  xlab("")

虽然简短，但希望对方差分析的介绍有助于突出方差分析如何用于检验 R 的组均值差异，以及方差分析如何只是 GLMs 的一个特例。最终，ANOVAs 是 GLM 的一个非常有限的特例，因为它们不允许包含连续的解释变量。接下来，我们将线性回归作为 GLM 的一个更灵活的特例进行检验，该特例适用于连续正态分布的结果，同时允许离散和连续的解释变量。

3.6 线性回归

注意

线性回归是连续、正态分布结果变量的 GLM 的特例。与方差分析不同，线性回归适用于离散和连续的解释变量/预测值。rms包中的ols()函数可以运行线性回归并打印综合汇总输出信息。

概念背景

线性回归是 GLMs 的另一种特殊情况，其中链接和反向链接函数只是恒等函数，即 η = g(μ) = μ 和μ= g^-1(η)=η并且结果被假设为正态分布。具体来说，分布假设写成 y ~ N ( μ ，σ)。这里的关键信息是 y 是一个向量，因为 μ 是一个常数。分布的平均值通常被称为它的位置，分散参数或标准偏差被称为它的尺度。另一类模型，位置比例模型，允许位置和比例参数作为数据的函数而变化，但对于 GLMs，我们将假设比例是恒定的。分布的另一种常见写法是残差分布，∑~N(0，σ)，其中∑= yμ。这样写强调了这样一个事实，即原始数据不需要遵循正态分布，它们只需要围绕期望值正态分布。它还强调了离差参数 σ 是残差的离差。也就是说， σ 捕捉期望值附近的离差。在最简单的 GLM 中，唯一的预测因子是一个常数项(截距)， σ 将与 y 的标准差相同，但如果回归可以解释 y 的一些或所有变化，那么 σ 将趋向于零。

基于中心极限定理，当样本大小收敛到无穷大时，回归系数与其标准误差之比的参数分布将收敛到正态分布。在线性回归中，我们可以考虑这样一个事实，即我们通常使用有限的样本来测试参数，而不是针对正态分布来测试参数(针对单个回归系数)。当有无限个自由度时，t 分布收敛到正态分布，而当有有限个自由度时，t 分布的尾部稍微重一些。在线性回归中，根据样本大小和参数个数(df=N–k_参数 )计算自由度。在后面的章节中，当我们检查其他类型的广义线性模型时，自由度不能很容易地计算出来，所以个体回归系数是根据标准正态分布而不是 t 分布来测试的。

除了测试单个回归系数之外，可以使用似然比测试来测试整个模型，该测试将我们拟合的模型的似然比与仅包含截距作为预测值的零模型进行比较。可能性有用的一个原因是，它们可以很容易地进行比较和测试，从而提供对模型的多个变量或其他限制的准确测试。

在线性回归中，一个常见的效应大小是模型(或单个预测因子)解释的结果中方差的百分比。所占的百分比方差称为R?? 2。在我们能够计算出R2 之前，我们需要几个棋子。我们将偏差平方和(SS) total 定义为结果与其总体预期或均值(SSTotal)的偏差，SSRegression 定义为我们的模型预测结果与结果总体预期的偏差，SSResidual 定义为结果与我们的模型预测值之间差异的 SS。这在以下等式中更正式地示出:

$S{S}_{Total}=\sum \limits_{i=1}^N{\left({y}_i-E(y)\right)}²$

(3.11)

$S{S}_{Regression}=\sum \limits_{i=1}^N{\left({\mu}_i-E(y)\right)}²$

(3.12)

$S{S}_{Residual}=\sum \limits_{i=1}^N{\left({y}_i-{\mu}_i\right)}²$

(3.13)

给定这些定义，对于具有正态分布结果的线性模型， R ² 可计算如下:

${R}²=\frac{S{S}_{Regression}}{S{S}_{Total}}=1-\frac{S{S}_{Residual}}{S{S}_{Total}}= cor{\left(\mathrm{y},\mu \right)}²$

(3.14)

除非存在无限的样本量，否则计算样本数据中 R ² 的公式是对总体 R ² 的有偏估计。因此，报告调整后的 R ² 也很常见，它考虑了模型自由度，以提供总体 R ² 的无偏估计。当我们训练模型并使用相同的数据对其进行测试时，对自由度的这种调整会调整人口中对方差的过于乐观的估计。当我们在R中讨论机器学习时，我们将讨论过拟合的概念以及使用单独的数据集进行模型估计和测试。在线性回归中，偏差趋于最小，因为与样本量相比，参数相对较少。随着观测值的预测值/参数数量的增加，过度拟合带来的问题和偏差变得更加棘手。 R ² 可用于估计整体模型的预测精度，但也可用于通过比较每次增加或减少一个预测器时模型 R ² 的变化量来量化单个预测器增加的量。

尽管 R ² 是迄今为止最常见的线性回归拟合或判别指数，但也存在其他指数。另一种选择是基于基尼指数的 g 指数[25]。g 指数和 R ² 的一个关键区别是，g 指数不是标准化的，所以它取决于结果的规模和预测因素。然而，与R2 一样，更高表示更好的区分度。

R 中的线性回归

随着我们转向应用和实际数据分析，我们将开始使用真实数据。在这一章中，我们将使用美国人不断变化的生活[45]研究数据。“数据设置”一节介绍了数据的读取和准备从技术上讲，数据具有采样权重，但为了简单起见，我们忽略这些权重。没有加权，分析仍然是正确的；它们只是不能反映抽样人口。

虽然不是对回归假设的直接检验，但了解结果的大致分布是有用的，即生活满意度。图 3-2 显示了正常曲线重叠的密度图。使用adjust = 2到testdistr()的参数，原始密度比默认值更加平滑。从图 3-2 我们可以看到，生活满意度近似正态分布，没有大的异常值。

img/439480_1_En_3_Fig2_HTML.png

图 3-2

正常密度覆盖(蓝线)的生活满意度密度图(黑线)

acl <- readRDS("advancedr_acl_data.RDS")

testdistr(acl$SWL_W1, "normal",
          varlab = "Satisfaction with Life", plot = FALSE,
          extremevalues = "theoretical",
          adjust = 2)$DensityPlot

R有内置函数来拟合线性回归，但是我们使用rms包中的ols()函数，因为它提供了方便的特性和更全面的默认输出。ols这个名字来源于线性回归的另一个名字:普通最小二乘法。这个名称是基于使用最小平方偏差作为估计回归系数的标准。

模型输出首先回显用于拟合模型的公式。它显示了观察值的数量、剩余标准偏差的估计值、*、*和整个模型的自由度。似然比检验同时提供了对所有预测值的统计显著性的检验，检验了至少一个系数显著不同于零的假设。鉴别指数包括 R ² 和调整后的 R ² 值以及 g 指数。构造残差的均值为零，但由于偏差或异常值，中值可能会有很大不同。检查最小和最大残差对于识别残差异常值也很有用。最后，表中显示了回归系数以及相应的标准误差、t 值和 p 值。

m.ols <- ols(SWL_W1 ~ Sex + AGE_W1 + SESCategory, data = acl, x = TRUE)
m.ols

## Linear Regression Model
##
##  ols(formula = SWL_W1 ~ Sex + AGE_W1 + SESCategory, data = acl,
##      x = TRUE)
##
##                  Model Likelihood    Discrimination
##                        Ratio Test           Indexes
##  Obs    3617    LR chi2    118.62    R2       0.032
##  sigma1.0355    d.f.            5    R2 adj   0.031
##  d.f.   3611    Pr(> chi2) 0.0000    g        0.213
##
##  Residuals
##
##       Min       1Q   Median       3Q      Max
##  -3.44270 -0.67206  0.01543  0.75504  2.36635
##
##
##                 Coef    S.E.   t     Pr(>|t|)
##  Intercept      -0.7057 0.0755 -9.35 <0.0001
##  Sex=(2) FEMALE  0.0308 0.0360  0.86 0.3921
##  AGE_W1          0.0103 0.0011  9.75 <0.0001
##  SESCategory=2  -0.0133 0.0447 -0.30 0.7654
##  SESCategory=3   0.2558 0.0482  5.31 <0.0001
##  SESCategory=4   0.2654 0.0635  4.18 <0.0001
##

使用texreg包，我们可以自动创建格式良好的表格。我们可以使用screenreg()函数为屏幕输出创建表格，使用htmlreg()函数为 HTML 输出创建表格，或者使用texreg()函数为 LATEX 创建表格。

下面的例子展示了如何制作 LATEX 表，在表 3-6 中给出。

表 3-6

统计模型

| |

模型 1

| | --- | --- | | 拦截 | -0.71(0.08)□ | | 性别=(2)女性 | 0.03 (0.04) | | 年龄 _W1 | 0.01(0.00)□ | | SESCategory=2 | −0.01 (0.04) | | SESCategory=3 | 0.26(0.05)□ | | SESCategory=4 | 0.27(0.06)□ | | 编号 obs。 | Three thousand six hundred and seventeen | | R ² | Zero point zero three | | 调整 R ² | Zero point zero three | | L.R | One hundred and eighteen point six two |

<【0.001】【t】**

texreg(m.ols, single.row = TRUE, label = "tglm1-olstex")

我们还可以研究关于模型的几个诊断。首先，我们可以通过使用方差膨胀因子(VIF)来探索任何共线性(解释变量之间的高度相关性)的影响。VIF 值接近 1 表示共线性的影响很小。非常高的 VIF 值可能表明包含高度相关的解释变量会增大参数协方差矩阵的方差，从而导致非常大的标准误差和置信区间。当包含两个非常相似的解释变量时，这种情况最常见。

vif(m.ols)

## Sex=(2) FEMALE         AGE_W1  SESCategory=2  SESCategory=3
##            1.0            1.2            1.4            1.5
##  SESCategory=4
##            1.3

接下来，我们用拟合值和残差创建一个数据表，并用它在图 3-3 中绘制残差图，以检查正态性。图 3-4 使用分位数回归【51，50，107】，正如我们在多元数据可视化章节中介绍的，通过绘制第 5、25、50、75 和 95 个百分点的分位数回归线来评估异方差性。这些线相对平坦，表明几乎没有异方差的证据。

img/439480_1_En_3_Fig4_HTML.png

图 3-4

用分位数回归检验残差与拟合值，以探索异方差性

img/439480_1_En_3_Fig3_HTML.png

图 3-3

绘制残差图以评估正态性

diagnostic.data <- data.table(
  fitted = fitted(m.ols),
  resid = residuals(m.ols))

testdistr(diagnostic.data$resid,
          "normal",
          varlab = "Satisfaction with Life Residuals", plot = FALSE,
          extremevalues = "theoretical",
          adjust = 2)$DensityPlot

ggplot(diagnostic.data, aes(fitted, resid)) +
  geom_point(alpha = .2, colour = "grey50") +
  geom_quantile(quantiles = .5, colour = 'black', size = 1) +
  geom_quantile(quantiles = c(.25, .75),
                colour = 'blue', linetype = 2, size = 1) +
  geom_quantile(quantiles = c(.05, .95),
                colour = 'black', linetype = 3, size = 1)

## Smoothing formula not specified. Using: y ~ x
## Smoothing formula not specified. Using: y ~ x
## Smoothing formula not specified. Using: y ~ x

我们用矩阵符号介绍了 GLM:

$\mu ={g}^{-1}\left(\eta \right)={g}^{-1}\left(\mathrm{X}\beta \right)$

(3.15)

为了更好地理解如何解释系数，使用正则代数将它写出来会很有帮助。在线性回归的情况下，我们可以通过去除反向链接函数来进一步简化，因为链接和反向链接(通常)是相同的函数。还通常对每个变量使用下标 i 来编写下面的等式，以指示对第 i 个个体执行操作。我们用粗体表示变量是向量，包含多个个体的数据。

$\mu ={\beta}_0+{\beta}_1\mathrm{Sex}+{\beta}_2\mathrm{Age}+{\beta}_3{\mathrm{SES}}_2+{\beta}_4{\mathrm{SES}}_3+{\beta}_5{\mathrm{SES}}_4$

(3.16)

每个系数都反映了预测因子中一个单位变化的结果的预期变化。一个单位的含义取决于每个预测器的规模。例如，年龄是用年来编码的，所以一个单位意味着一年。性别是虚拟编码的，所以一个单位代表了男女之间的区别。社会经济地位(SES)被编码为四分位数，参考类别(省略)是 SES 的最低四分位数。因此，举例来说，年长一岁意味着对生活的满意度提高 0.01%。

我们注意到线性回归是 GLM 的特例。R有内置函数glm()，适合 GLMs。glm()功能的好处是您可以对许多特定类型的 glm 使用相同的功能。公式界面与ols()相同，但glm()允许用户指定不同的分布和链接函数。如果我们简单地打印存储的 GLM 对象，我们得到的输出很少。为了得到一个好的摘要，我们需要使用summary()函数。虽然这在计算上稍微更有效，但是根据我们的经验，大多数用户倾向于想要这个输出，所以ols()函数默认提供一个好的输出摘要是很方便的。glm()函数的另一个缺点是它不显示R2 值。这是因为方差并不适用于所有类型的 GLMs。与独立模型相比，没有对整体模型的默认测试，也没有关于每个变量缺失值数量的信息。

m.glm <- glm(SWL_W1 ~ Sex + AGE_W1 + SESCategory,
             data=acl, family = gaussian(link="identity"))
m.glm

##
## Call:  glm(formula = SWL_W1 ~ Sex + AGE_W1 + SESCategory, family = gaussian(link = "identity"),
##     data = acl)
##
## Coefficients:
##   (Intercept)  Sex(2) FEMALE         AGE_W1   SESCategory2
##       -0.7057         0.0308         0.0103        -0.0133
##  SESCategory3   SESCategory4
##        0.2558         0.2654
##
## Degrees of Freedom: 3616 Total (i.e. Null);  3611 Residual
## Null Deviance:     4000
## Residual Deviance: 3870      AIC: 10500

summary(m.glm)

##
## Call:
## glm(formula = SWL_W1 ~ Sex + AGE_W1 + SESCategory, family = gaussian(link = "identity"),
##     data = acl)
##
## Deviance Residuals:
##    Min      1Q  Median      3Q     Max
## -3.443  -0.672   0.015   0.755   2.366
##
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)   -0.70565    0.07549   -9.35  < 2e-16 ***
## Sex(2) FEMALE  0.03083    0.03602    0.86     0.39
## AGE_W1         0.01030    0.00106    9.75  < 2e-16 ***
## SESCategory2  -0.01333    0.04467   -0.30     0.77
## SESCategory3   0.25579    0.04819    5.31  1.2e-07 ***
## SESCategory4   0.26544    0.06353    4.18  3.0e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1
##
## (Dispersion parameter for gaussian family taken to be 1.1)
##
##     Null deviance: 4000.7  on 3616  degrees of freedom
## Residual deviance: 3871.6  on 3611  degrees of freedom
## AIC: 10525
##
## Number of Fisher Scoring iterations: 2

输出与R的内置lm()函数非常相似，该函数专门用于线性回归模型。它确实增加了对R2 的估计，但是同样需要使用summary()并且没有显示每个变量的缺失值。

通常，会评估多个相关的模型。例如，可能有几个焦点预测值，目标是当一些额外的变量(协变量)在模型之内或之外时，检查它们的影响如何变化。其他时候，某些预测值可能会被包括在内，而其他预测值只有在它们具有统计显著性或改善模型性能时才会被包括在内。尝试其他功能形式也很常见。例如，变量可以作为线性效果输入，或者作为线性效果与相同变量的平方(即，二次)一起添加。在本章的前面，我们展示了如何更新R中的公式。update()函数也有许多模型的方法。它的工作方式与公式上的update()类似，除了改变模型公式之外，它还重新拟合模型。为了更新和查看存储在对象中的结果，我们可以将整个调用放在括号中，这将强制打印。在下面的代码中，我们更新了基本模型，并向模型中添加了雇佣状态。

(m.ols2 <- update(m.ols, . ~ . + Employment_W1))

## Linear Regression Model
##
##  ols(formula = SWL_W1 ~ Sex + AGE_W1 + SESCategory + Employment_W1,
##      data = acl, x = TRUE)
##
##                  Model Likelihood    Discrimination
##                        Ratio Test           Indexes
##  Obs    3617    LR chi2    173.43    R2       0.047
##  sigma1.0286    d.f.           12    R2 adj   0.044
##  d.f.   3604    Pr(> chi2) 0.0000    g        0.252
##
##  Residuals
##
##       Min       1Q   Median       3Q      Max
##  -3.50191 -0.66381  0.03265  0.74125  2.55188
##
##
##                             Coef    S.E.   t     Pr(>|t|)
##  Intercept                  -1.1197 0.1244 -9.00 <0.0001
##  Sex=(2) FEMALE              0.0109 0.0387  0.28 0.7776
##  AGE_W1                      0.0092 0.0013  6.83 <0.0001
##  SESCategory=2              -0.0253 0.0451 -0.56 0.5746
##  SESCategory=3               0.2179 0.0498  4.37 <0.0001
##  SESCategory=4               0.2174 0.0655  3.32 0.0009
##  Employment_W1=(2) 2500+HRS  0.5832 0.1098  5.31 <0.0001
##  Employment_W1=(3) 15002499  0.4675 0.0985  4.75 <0.0001
##  Employment_W1=(4) 500-1499  0.5497 0.1085  5.07 <0.0001
##  Employment_W1=(5) 1-499HRS  0.6135 0.1250  4.91 <0.0001
##  Employment_W1=(6) RETIRED   0.5345 0.0962  5.55 <0.0001
##  Employment_W1=(7) UNEMPLOY  0.2498 0.1233  2.03 0.0428
##  Employment_W1=(8) KEEP HS   0.6218 0.0991  6.28 <0.0001
##

为了测试 SES 或就业总体上是否重要，我们需要同时测试所有的伪代码。这可以通过比较两个模型或调用拟合模型上的anova()功能来完成。由此产生的方差分析表显示了对社会经济地位的三自由度测试和对就业的七自由度测试，两者总体上具有统计学意义。

anova(m.ols2)

##                 Analysis of Variance          Response: SWL_W1
##
##  Factor        d.f. Partial SS MS     F     P
##  Sex              1 8.4e-02     0.084  0.08 0.78
##  AGE_W1           1 4.9e+01    49.364 46.65 <.0001
##  SESCategory      3 3.8e+01    12.631 11.94 <.0001
##  Employment_W1    7 5.8e+01     8.318  7.86 <.0001
##  REGRESSION      12 1.9e+02    15.608 14.75 <.0001
##  ERROR         3604 3.8e+03     1.058

没有性别差异，所以我们可以考虑放弃性。我们还可以探索潜在的相互作用，比如年龄和社会经济地位之间的相互作用。这两者都可以在一个步骤中完成，即更新模型。我们再次使用括号强制打印。

(m.ols3 <- update(m.ols2, . ~ . + AGE_W1 * SESCategory - Sex))

## Linear Regression Model
##
##  ols(formula = SWL_W1 ~ AGE_W1 + SESCategory + Employment_W1 +
##      AGE_W1:SESCategory, data = acl, x = TRUE)
##
##                  Model Likelihood    Discrimination
##                        Ratio Test           Indexes
##  Obs    3617    LR chi2    189.72    R2       0.051
##  sigma1.0266    d.f.           14    R2 adj   0.047
##  d.f.   3602    Pr(> chi2) 0.0000    g        0.256
##
##  Residuals
##
##       Min       1Q   Median       3Q      Max
##  -3.37389 -0.65254  0.04075  0.72383  2.60671
##
##
##                             Coef    S.E.   t     Pr(>|t|)
##  Intercept                  -1.2652 0.1568 -8.07 <0.0001
##  AGE_W1                      0.0116 0.0021  5.56 <0.0001
##  SESCategory=2              -0.0495 0.1566 -0.32 0.7518
##  SESCategory=3               0.6213 0.1678  3.70 0.0002
##  SESCategory=4               0.7440 0.2128  3.50 0.0005
##  Employment_W1=(2) 2500+HRS  0.5628 0.1092  5.15 <0.0001
##  Employment_W1=(3) 15002499  0.4643 0.0984  4.72 <0.0001
##  Employment_W1=(4) 500-1499  0.5592 0.1083  5.16 <0.0001
##  Employment_W1=(5) 1-499HRS  0.6284 0.1247  5.04 <0.0001
##  Employment_W1=(6) RETIRED   0.5280 0.0961  5.50 <0.0001
##  Employment_W1=(7) UNEMPLOY  0.2812 0.1232  2.28 0.0225
##  Employment_W1=(8) KEEP HS   0.6293 0.0978  6.43 <0.0001
##  AGE_W1 * SESCategory=2      0.0009 0.0026  0.35 0.7248
##  AGE_W1 * SESCategory=3     -0.0077 0.0029 -2.62 0.0088
##  AGE_W1 * SESCategory=4     -0.0107 0.0041 -2.61 0.0090
##

年龄和社会经济地位之间似乎确实存在相互作用。对于模型中的相互作用，年龄系数是最低四分位数 SES(伪编码 SES 的参照组，因此是 SES 与年龄相互作用的参照组)的年龄和对生活满意度的斜率。相互作用的系数表明，在较高的社会经济地位类别中，年龄与生活满意度之间的关系较低。为了更好地理解这种相互作用，我们可以通过显示生活满意度的期望值如何随年龄变化来绘制图表，结果按社会经济地位类别进行分类。

在R中绘制许多回归模型结果的快速方法是使用visreg包【15】中的visreg()函数。该功能允许快速预测和置信区间的生成和图形化。它有合理的缺省值，比如保存变量，你没有在它们的中间值或连续变量和分类变量的模式中绘图。visreg()函数最少只需要两个参数。第一，模型是什么，第二，你想在 x 轴上绘制哪个变量(xvar自变量)？我们将使用它来请求 x 轴上的年龄，并扩展它来请求不同的行by SESCategory，这将制作一个交互图。

默认情况下,visreg()包括部分残差，一个地毯图，并将交互图分成单独的面板。我们使用overlay = TRUE将所有面板合并成一个图，使用partial = FALSE关闭绘制部分残差，这样我们只有预测的线，使用rug = FALSE关闭显示数据点落在 x 轴上的地毯图，使用xlab和ylab添加一些更好的 x 轴和 y 轴标签，最后，更改线型，这样当图表不以彩色打印时，通过设置line = list(lty = 1:4)仍然可以阅读。所有定制结果如图 3-5 所示。

img/439480_1_En_3_Fig5_HTML.png

图 3-5

按社会经济地位类别对各年龄段生活满意度的估计。阴影区域表示回归估计的 95%置信区间。

plot(visreg(m.ols3, xvar = "AGE_W1", by = "SESCategory",
            plot = FALSE),
     overlay = TRUE, partial = FALSE, rug = FALSE,
     xlab = "Age (years)", ylab = "Predicted Life Satisfaction",
     line = list(lty = 1:4))

由于所有置信区间重叠，该图仍然有点混乱。就个人理解而言，这些可能是有帮助的。在演示文稿中使用时，可能会很难看到线条。我们可以使用另一个参数band = FALSE来关闭置信区间。我们可以通过将图表灰度化，传递四种颜色，用于四条线中的每一条线，来进一步修改它以便发布。进一步定制的结果如图 3-6 所示。

img/439480_1_En_3_Fig6_HTML.png

图 3-6

按社会经济地位类别对各年龄段生活满意度的估计。置信区间已删除。

plot(visreg(m.ols3, xvar = "AGE_W1", by = "SESCategory",
            plot = FALSE),
     overlay = TRUE, partial = FALSE, rug = FALSE,
     xlab = "Age (years)", ylab = "Predicted Life Satisfaction",
     line = list(
       lty = 1:4,
       col = c("black", "grey75", "grey50", "grey25")),
     band = FALSE)

使用visreg()功能既快速又相对容易，所以在大多数设置中它是一个很好的选择，当然也是帮助你自己理解结果的一种方式。为了更好的控制，我们可以手动制作同样的图表。为此，我们需要获得各种年龄和社会经济地位类别的预测值。我们还需要将模型中的其他变量，就业，保持在某个值上。在R中很容易从模型中获得预测值，但是我们首先需要创建一个小型数据集，其中包含我们希望用作预测输入的所有值。这可以通过使用expand.grid()功能轻松完成。重要的是，因子具有与模型中相同的水平，这最容易通过从真实数据中的因子提取levels()来实现。

newdata <- as.data.table(expand.grid(
  AGE_W1=quantile(acl$AGE_W1, .1):quantile(acl$AGE_W1, .9),
  SESCategory = factor(1:4, levels = levels(acl$SESCategory)),
  Employment_W1 = factor("(3) 15002499",
    levels = levels(acl$Employment_W1))))
newdata

##      AGE_W1 SESCategory Employment_W1
##   1:     30           1  (3) 15002499
##   2:     31           1  (3) 15002499
##   3:     32           1  (3) 15002499
##  ---
## 186:     74           4  (3) 15002499
## 187:     75           4  (3) 15002499
## 188:     76           4  (3) 15002499

现在我们可以使用predict()函数生成预测值。我们可以只提取预测值或者每个预测值的预测值和标准误差。标准误差有助于计算每个预测的置信区间，并显示估计的不确定性。为了得到标准误差，我们指定，se.fit = TRUE。结果是一个列表，其中第一个元素包含预测值的向量，第二个元素包含标准误差的向量，我们将其存储回数据表中。

newdata[, c("SWL_W1", "SE") :=
          predict(m.ols3, newdata = newdata, se.fit = TRUE)]
newdata

##      AGE_W1 SESCategory Employment_W1 SWL_W1    SE
##   1:     30           1  (3) 15002499 -0.453 0.076
##   2:     31           1  (3) 15002499 -0.441 0.075
##   3:     32           1  (3) 15002499 -0.430 0.073
##  ---
## 186:     74           4  (3) 15002499  0.014 0.121
## 187:     75           4  (3) 15002499  0.015 0.124
## 188:     76           4  (3) 15002499  0.016 0.128

置信区间的计算方法如下

$Estimate\pm SEx{z}_{\alpha /2}$

(3.17)

z 是指单位正态分布的分位数(通常称为 z 得分)。 z _{α/ 2} 是基于期望的α水平的分位数(例如，95%置信区间为 0.05)。更准确地说，可以使用具有适当自由度的 t 分布的分位数，尽管对于如此大的样本，t 分布实际上是正态的。这可以使用qnorm()功能在R中获得。

print(qnorm(.05/2), digits = 7)

## [1] -1.959964

print(qnorm(1 - (.05/2)), digits = 7)

## [1] 1.959964

接下来，我们使用ggplot2包和来自cowplot包的主题创建一个预测值的图表。代码有点复杂qnorm()功能，但产生了如图 3-7 所示的出版物质量图。虽然这比使用visreg()有些繁琐，但它让我们可以完全控制我们正在计算的预测值，并允许我们在绘制图表或演示之前对预测进行进一步的分析或工作。

img/439480_1_En_3_Fig7_HTML.png

图 3-7

按社会经济地位类别对各年龄段生活满意度的估计。阴影区域表示回归估计的 95%置信区间。

ggplot(newdata, aes(AGE_W1, SWL_W1, linetype=SESCategory)) +
  geom_ribbon(aes(ymin = SWL_W1 + SE * qnorm(.025),
                  ymax = SWL_W1 + SE * qnorm(.975)),
              alpha = .2) +
  geom_line(size = 1) +
  scale_x_continuous("Age (years)") +
  ylab("Satisfaction with Life") +
  theme_cowplot() +
  theme(
    legend.position = c(.8, .16),
    legend.key.width = unit(2, "cm"))

因为它们依赖于回归系数和标准误差(se)，置信区间与 p 值密切相关。然而，它们是显示真实总体回归系数估计的不确定性的一种有用方法。在R中，我们可以使用confint()函数计算每个回归系数的 95%置信区间。

confint(m.ols3)

##                              2.5 %  97.5 %
## Intercept                  -1.5726 -0.9579
## AGE_W1                      0.0075  0.0157
## SESCategory=2              -0.3566  0.2575
## SESCategory=3               0.2922  0.9504
## SESCategory=4               0.3267  1.1612
## Employment_W1=(2) 2500+HRS  0.3487  0.7768
## Employment_W1=(3) 15002499  0.2714  0.6572
## Employment_W1=(4) 500-1499  0.3468  0.7715
## Employment_W1=(5) 1-499HRS  0.3839  0.8729
## Employment_W1=(6) RETIRED   0.3396  0.7163
## Employment_W1=(7) UNEMPLOY  0.0396  0.5227
## Employment_W1=(8) KEEP HS   0.4376  0.8211
## AGE_W1 * SESCategory=2     -0.0041  0.0059
## AGE_W1 * SESCategory=3     -0.0134 -0.0019
## AGE_W1 * SESCategory=4     -0.0186 -0.0027

高性能线性回归

到目前为止，我们已经关注了具有全面输出的便利功能。这大概是大多数用户大部分时间需要的。线性回归在现代计算机上是如此之快，以至于计算时间在大多数时候都不是问题。然而，在某些情况下，计算速度是一个问题。Bootstrapping 是一个过程，我们将在稍后讨论机器学习时进行更深入的讨论，但简单来说，它需要从数据集重复采样，并估计重采样数据的一些参数，以生成经验参数分布。对于自举，我们可能只想提取回归系数，这可以使用coef()函数来完成。

人们经常采取数百或数千个 bootstrap 样本。出于时间和演示的原因，我们只取 500。首先，我们创建一个仅包含我们的变量的小型数据集，因为我们使用英特尔的 MKL 线性代数库，为了获得更纯粹的时间估计，我们将其设置为仅使用一个内核。这在单核机器或不使用 MKL 的机器上可以忽略。

tmpdat <- na.omit(acl[, .(SWL_W1, AGE_W1, SESCategory, Employment_W1)])
## use if using Microsoft R Open with Intel's MKL linear algebra library
setMKLthreads(1)

我们实际的代码相当简单。我们使用system.time()函数来跟踪需要多长时间，然后使用sapply()在 1 到 500 之间循环，创建要使用的行的索引，然后拟合我们的模型并提取系数。

set.seed(12345)
t1 <- system.time(ols.boot <- sapply(1:500, function(i) {
  index <- sample(nrow(tmpdat),
                  size = nrow(tmpdat), replace = TRUE)
  coef(ols(SWL_W1 ~ AGE_W1 * SESCategory + Employment_W1,
           data = tmpdat[index]))
}))

t1

##    user  system elapsed
##    4.27    0.06    4.33

使用ols()函数，花了 4.33 秒完成——没有长到不可能，但是长到足以显著减慢交互数据分析。对于任何实际的应用，我们可能至少需要几千次引导程序重采样。时间将随着重新采样的次数以线性方式增加，因此 10，000 次采样大约需要 86.6 秒。接下来，我们使用RcppEigen包【4】中的fastLm()函数。它使用C++来实现线性模型，以便更快更有效。

set.seed(12345)
t2 <- system.time(rcpp.boot1 <- sapply(1:500, function(i) {
  index <- sample(nrow(tmpdat), size = nrow(tmpdat), replace = TRUE)
  coef(fastLm(SWL_W1 ~ AGE_W1 * SESCategory + Employment_W1, data = tmpdat[index]))
}))

t2

##    user  system elapsed
##     2.5     0.0     2.5

现在，总时间减少到 2.52 秒，因此对于 10，000 次重新采样，大约需要 50.4 秒。最后，我们使用fastLmPure()函数，也来自RcppEigen包。fastLmPure()函数并不智能，它要求用户将结果作为向量和模型矩阵传递，而不是使用公式接口。我们将结果向量和模型矩阵的显式计算包含在我们的系统计时中，然后将自举重采样指数应用于这些预先计算的矩阵。

set.seed(12345)
t3 <- system.time({
  y <- tmpdat[, SWL_W1]
  X <- model.matrix(~ AGE_W1 * SESCategory + Employment_W1, data = tmpdat)
  N <- nrow(tmpdat)
  rcpp.boot2 <- sapply(1:500, function(i) {
    index <- sample.int(N, size = N, replace = TRUE)
    fastLmPure(X = X[index, ], y = y[index])$coefficients
  })
})

t3

##    user  system elapsed
##    0.48    0.02    0.50

使用这种方法，分析只需 0.5 秒。由于它如此之快，我们可以使用 10，000 个重采样轻松地重新运行它。

set.seed(12345)
t4 <- system.time({
  y <- tmpdat[, SWL_W1]
  X <- model.matrix(~ AGE_W1 * SESCategory + Employment_W1, data = tmpdat)
  N <- nrow(tmpdat)
  rcpp.boot3 <- sapply(1:10000, function(i) {
    index <- sample.int(N, size = N, replace = TRUE)
    fastLmPure(X = X[index, ], y = y[index])$coefficients
  })
})

t4

##    user  system elapsed
##    9.95    0.21   10.15

使用 10，000 个重采样需要 10.15 秒。有了并行处理，这个数字还可以进一步降低。我们不使用fastLmPure()进行交互式数据分析，但是对于计算量很大的任务，比如 bootstrapping，或者如果您正在尝试数百种不同的模型，速度的提高是有意义的。最后，我们可以检查使用all.equal()从所有模型中我们确实得到了相同的结果。设置check.attributes = FALSE忽略名称，因为ols()对虚拟系数的命名略有不同。

all.equal(ols.boot, rcpp.boot1, check.attributes = FALSE)

## [1] TRUE

all.equal(ols.boot, rcpp.boot2, check.attributes = FALSE)

## [1] TRUE

3.7 控制混杂因素

在科学领域，GLMs 通常用于研究一个变量对另一个变量的潜在影响。例如，自我效能感是指某人认为自己有能力改变或控制自己的生活的倾向。研究表明，自我效能高的人更容易改变行为(例如，开始锻炼计划、戒烟、报名并完成大学学位)。如果你想象一个自我效能感低的人，这是有道理的:他们倾向于相信他们不会成功地做出改变，并且倾向于认为他们的行为和环境不受他们的控制(例如，受环境、强大的他人等的控制)。).不管一个人是否真的能够或不能影响自己的生活，如果他们认为他们不能，他们可能会更快地放弃，更没有动力去尝试，因此更不可能开始或维持任何行为或追求自己的目标。现在假设自我效能高的人也不太可能经历抑郁症状。我们在引言中提到的 ACL 数据包括一个捕获自我效能的变量和另一个捕获两波数据收集中的抑郁症状的变量。一个自然的起点是测试第一波的自我效能是否能预测第二波的抑郁症状。下面的代码测试了这一点，我们可以看到，确实有一个统计上显著的负相关，在第一波自我效能较高的人往往在第二波有较低的抑郁症状。结果如表 3-7 所示。

表 3-7

统计模型

| |

模型 1

| | --- | --- | | 拦截 | Zero point zero two | | | (0.02) | | 自我效能 _W1 | -0.36 * | | | (0.02) | | 编号 obs。 | Two thousand eight hundred and sixty-seven | | R ² | Zero point one three | | 调整 R ² | Zero point one three | | L.R | Three hundred and ninety-nine point seven one |

<【0.001】【t】**

m0 <- ols(CESD11_W2 ˜ SelfEfficacy_W1, data = acl)

texreg(m0, label = "tglm1-olsunadj")

如果这纯粹是一个预测模型，我们可能会对目前的结果感到满意。然而，从科学的角度来看，这不足以表明两个变量是相关的。关联并不意味着一个变量导致另一个变量。这是一个重要的区别。如果自我效能感导致较低的抑郁症状，那么如果我们可以干预并增加某人的自我效能感，我们会期望他们有较少的抑郁症状。然而，如果自我效能不是原因，而只是一个预测因素或与抑郁症状有关，那么改变自我效能可能对抑郁症状没有影响。

形式上，这引入了一个通常称为混杂的概念。出于多种原因，两个变量可能相互关联。关联的一个原因是一个变量导致另一个变量。然而，在一个不准确的现实模型中，你也可以找到两个变量之间的关联，因为一些第三变量导致了这两个变量。例如，假设有慢性健康问题会导致较低的自我效能和较高的抑郁症状。如果慢性健康问题的存在不能以某种方式解释，那么自我效能和抑郁症状之间似乎存在关联。然而，一旦考虑到慢性健康问题的影响，自我效能和抑郁症状可能没有关联。

表示不同可能的因果配置的一种常见方式是通过因果图。这些可以理解为代表不同变量的圆圈和指示哪个变量导致哪个变量的有向箭头。对于这种模型的温和介绍和因果推理的更深入的概述，见[81]。因果图示例如图 3-8 所示。在图 3-8 中， Z 是 X 和 Y 的共同原因。从我们的概念示例来看， Z 将是慢性病， X 将是自我效能，而 Y 将是第二波的抑郁症状。如果不考虑 Z ，将会获得对 X 和 Y 的关联的不准确的、有偏差的估计。

img/439480_1_En_3_Fig8_HTML.png

图 3-8

示例图，其中变量 Z 是 X 和 y 的共同原因。如果不考虑 Z，则 X 和 y 之间似乎存在关联。

不准确的模型还会以其他方式产生有偏差的估计。图 3-9 显示了另一个图形模型，其中 Z 现在是 X 和 Y 的共同结果，称为碰撞器变量，因为来自 X 和 Y 的路径在 Z 上碰撞。如果事实如图 3-9 所示，那么直接测试 X 和 Y 之间的关联是合适的。然而，如果我们试图将 Z 作为混杂变量进行控制，通过将 Z 添加到预测 Y 的模型中，远未减少 Z 的混杂，该模型将导致 X 和 Y 之间的虚假关联。从这两个例子中得到的教训是，如果 Z 是一个共同的原因(例如，图 3-8 ，不恰当的排除路径会导致偏差。相反，如果 Z 是碰撞体(如图 3-9 )，包含路径的不恰当会导致偏差。

img/439480_1_En_3_Fig9_HTML.png

图 3-9

示例图，其中变量 Z 是 X 和 Y 的碰撞体(常见结果)。如果在检查 X 和 Y 之间的关联时忽略 Z，则该关联将被准确估计。但是，如果碰撞器 Z 被条件化，那么它将打开 X 和 y 之间的关联。

最后一个例子如图 3-10 所示。这里的 Z 是 X 冲击 Y 的机构。换句话说， Z 将 X 的效果传递给 Y 。在这种情况下，我们可以测试 X 通过 Z 对 Y 的间接影响，但是在模型中没有 Z 的情况下，我们应该会看到 X 与 Y 的关联，但是一旦 Z 被添加到模型中， X 应该不再与 Y 直接关联，因为

img/439480_1_En_3_Fig10_HTML.png

图 3-10

示例图，其中变量 Z 是 X 和 Y 的机制(中介)。因为 Z 将 X 的影响传递到 Y，所以对 Z 的调节将消除 X 和 Y 的关联。如果在检查 X 和 Y 之间的关联时忽略 Z，则该关联将被准确估计。

这些想法中的许多，无论是横截面数据还是纵向数据，都在为边际结构模型(MSMs) [80]开发的理论中有更深入的阐述。MSM 通常使用逆概率加权估值器进行估计[80]。在离散预测器的上下文中， X ，计算第个人的逆概率权重(IPW)的基本等式是

${w}_i=\frac{1}{P\left(X={x}_i|Z={z}_i\right)}$

(3.18)

有时，这些基本权重被扩展以计算所谓的稳定权重，这些权重只是针对 X 采用任何特定值的边际概率进行调整，由以下等式给出

$s{w}_i=\frac{P\left(X={x}_i\right)}{P\left(X={x}_i|Z={z}_i\right)}$

(3.19)

虽然逆概率加权(IPW)估计的使用来自于焦点预测因子是二分法处理的模型，但是 IPW 思想也很容易推广到连续预测因子。简而言之，这个想法是假设你正在研究一个焦点预测因子， X ，以及它与一些结果的关联， Y ，但是有一些已知的混杂变量， Z 。调整 Z 效果的一种方法是计算给定特定 Z 时 X 的概率。在连续的情况下，使用相同的基本思想，除了不是依赖于概率质量函数，我们必须使用概率密度函数用于一些假设的分布。此外，标准做法是使用稳定的砝码。假设我们相信 X 遵循正态分布， X ~ N ( μ，σ )，我们估计连续预测器的稳定权重为

$s{w}_i=\frac{f_X\left(X;{\mu}_1,{\sigma}_1\right)}{f_{X\mid Z}\left(X|Z={z}_i;{\mu}_2,{\sigma}_2\right)}$

(3.20)

其中 f _X (⋅)是正态分布的概率密度函数。本质上，通过使用基于 X 的无条件模型的概率密度函数估计，我们可以控制 X 的不同变化量。

使用 IPWs，我们可以获得我们感兴趣的模型，并使用 IPWs 对其进行估计，并且如果正确指定了用于生成权重的模型以及关联 X 和 Y 的模型，这被证明可以渐近地产生预测值 X 与结果 Y 的关联的无偏估计。权重模型(例如，未能包括所需的混杂变量，或未能指定正确的函数形式，如指定线性关联，而实际上它是二次的)或焦点模型的错误指定将导致有偏估计。在持续暴露的情况下，一个建议是缩小或调整底部和顶部 1%的权重，以减少与极端权重相关的噪声[68]。关于构建 IPWs 的更多信息可在参考文献[22]和[68]中获得。

为了在R中构造 IPWs，我们可以使用优秀的ipw包【96】。首先，我们必须决定我们希望在 IPW 调整模型中调整哪些变量。通常有几个步骤。也许作为我们的第一步，我们包括一些潜在的共同原因，我们认为(基于理论，以前的数据，希望比随机直觉稍微强一点)是自我效能和抑郁症状的潜在共同原因。性别、种族/民族和年龄可能是很好的选择，特别是因为这些都不可能是由自我效能或抑郁症状引起的。我们也可能包括一些慢性疾病。

为了计算单个时间点的 IPWs，我们使用ipw包中的ipwpoint()函数。该函数要求我们指定暴露的变量名；预期分布，在这种情况下我们假设正态分布；然后是分子和分母概率密度函数的模型。最后，当然我们必须告诉ipwpoint()使用哪个数据集。我们将结果保存为R中的变量w。可以从结果中访问 IPWs，如ipw.weights。在图 3-11 中显示了原始重量和将底部和顶部 1%百分点相加后的重量的快速分布图。

img/439480_1_En_3_Fig11_HTML.png

图 3-11

自我效能的原始和修整的反向概率权重

## weights
w <- ipwpoint(
  exposure = SelfEfficacy_W1,
  family = "gaussian",
  numerator = ~ 1,
  denominator = ~ 1 + Sex + RaceEthnicity + AGE_W1 + NChronic12_W1,
  data = acl)

plot_grid(
  testdistr(w$ipw.weights, plot = FALSE)$DensityPlot,
  testdistr(winsorizor(w$ipw.weights, .01),
            plot = FALSE)$DensityPlot,
  ncol = 1)

一旦我们有了权重，我们就可以通过将权重传递给ols()的weights参数来估计一个加权模型。使用 IPWs 将性别、种族/民族、年龄和慢性病作为潜在的混杂因素进行调整。为了进行比较，我们包括了未调整的模型，结果如表 3-8 所示。

表 3-8

统计模型

| |

模型 1

模型 2

| | --- | --- | --- | | 拦截 | Zero point zero two | Zero point zero two | | | (0.02) | (0.02) | | 自我效能 _W1 | 0.36 | 0.32 | | | (0.02) | (0.02) | | 编号 obs。 | Two thousand eight hundred and sixty-seven | Two thousand eight hundred and sixty-seven | | R ² | Zero point one three | Zero point one one | | 调整 R ² | Zero point one three | Zero point one one | | L.R | Three hundred and ninety-nine point seven one | Three hundred and twenty-five point two three |

^{** * *}<【0.001】、**<<【0.01】、

## unweighted, unadjusted
m0 <- ols(CESD11_W2 ~ SelfEfficacy_W1, data = acl)

## weighted, adjusted
m1 <- ols(CESD11_W2 ~ SelfEfficacy_W1, data = acl,
  weights = winsorizor(w$ipw.weights, .01))

texreg(list(m0, m1),
       label = "tglm1-weight1")

作为敏感性分析，我们也可以尝试进一步调整其他因素，我们认为这些因素可能是混杂因素，也可能是传递自我效能对抑郁症状影响的介质或机制。在这种情况下，我们添加了社会经济地位类别、就业、身体质量指数、吸烟状况和身体活动类别。我们估计权重，然后使用这些新的权重重新估计模型。表 3-9 中给出了未调整模型与部分和完全调整模型的对比。

表 3-9

统计模型

| |

模型 1

模型 2

模型 3

| | --- | --- | --- | --- | | 拦截 | Zero point zero two | Zero point zero two | Zero point zero two | | | (0.02) | (0.02) | (0.02) | | 自我效能 _W1 | 0.36 | 0.32 | 0.29 ^*** | | | (0.02) | (0.02) | (0.02) | | 编号 obs。 | Two thousand eight hundred and sixty-seven | Two thousand eight hundred and sixty-seven | Two thousand eight hundred and sixty-seven | | R ² | Zero point one three | Zero point one one | Zero point zero nine | | 调整 R ² | Zero point one three | Zero point one one | Zero point zero nine | | L.R | Three hundred and ninety-nine point seven one | Three hundred and twenty-five point two three | Two hundred and sixty-one point five two |

^{** * *}<【0.001】、**<<【0.01】、

# weighted, fully adjusted
w2 <- ipwpoint(
  exposure = SelfEfficacy_W1,
  family = "gaussian",
  numerator = ~ 1,
  denominator = ~ 1 + Sex + RaceEthnicity + AGE_W1 + NChronic12_W1 +
    SESCategory + Employment_W1 + BMI_W1 + Smoke_W1 + PhysActCat_W1,
  data = acl)

m2 <- ols(CESD11_W2 ~ SelfEfficacy_W1, data = acl,
  weights = winsorizor(w2$ipw.weights, .01))

texreg(list(m0, m1, m2),
       label = "tglm1-weight2")

调整潜在混杂因素的另一种方法是简单地将潜在混杂因素添加到模型中。以下代码显示了这样的例子，称为型号m1b和m2b,“b”表示它是 IPW“型号 1”和“型号 2”的替代产品:

m1b <- ols(CESD11_W2 ~ Sex + RaceEthnicity + AGE_W1 +
  NChronic12_W1 + SelfEfficacy_W1,
  data = acl)

m2b <- ols(CESD11_W2 ~ Sex + RaceEthnicity + AGE_W1 +
  NChronic12_W1 + SESCategory +
  Employment_W1 + BMI_W1 + Smoke_W1 + PhysActCat_W1 +
  SelfEfficacy_W1, data = acl)

最后，有些人建议使用所谓的双重稳健估计量。双重稳健估计简单地包括两个 IPW 权重，然后包括用于将权重显式构建到模型中的相同混杂。这方面的例子显示在下面的代码中，标记为m1c和m2c，因为它们是我们的两个调整模型的另一个变体。

m1c <- ols(CESD11_W2 ~ Sex + RaceEthnicity + AGE_W1 +
  NChronic12_W1 + SelfEfficacy_W1,
  data = acl,
  weights = winsorizor(w$ipw.weights, .01))

m2c <- ols(CESD11_W2 ~ Sex + RaceEthnicity + AGE_W1 +
  NChronic12_W1 + SESCategory +
  Employment_W1 + BMI_W1 + Smoke_W1 + PhysActCat_W1 +
  SelfEfficacy_W1, data = acl,
  weights = winsorizor(w2$ipw.weights, .01))

为了比较这些不同的方法，我们可以从每个模型中提取估计值和置信区间，然后绘制图表，这样我们就可以很容易地将差异可视化。这显示在下面的代码中，结果如图 3-12 所示。在这种情况下，所有的结果都非常相似。当变量只存在于一个时间点时，通常就是这种情况。然而，在特别推荐 IPWs 的边际结构模型的情况下，这些方法可能会有更大的差异。

img/439480_1_En_3_Fig12_HTML.png

图 3-12

不同模型中自我效能与抑郁症状关联的估计值和置信区间的比较。Covs =协变量调整模型。逆概率权重调整模型。Covs + IPW =在模型中再次明确包含反向概率权重和相同潜在混杂的模型。

## write an extract function
extractor <- function(obj, label) {
  b <- coef(obj)
  ci <- confint(obj)
  data.table(
    Type = label,
    B = b[["SelfEfficacy_W1"]],
    LL = ci["SelfEfficacy_W1", "2.5 %"],
    UL = ci["SelfEfficacy_W1", "97.5 %"])
}

allresults <- rbind(
  extractor(m0,  "M0: Unadjusted"),
  extractor(m1,  "M1: Partial IPW"),
  extractor(m1b, "M1: Partial Covs"),
  extractor(m1c, "M1: Partial Covs + IPW"),
  extractor(m2,  "M2: Full IPW"),
  extractor(m2b, "M2: Full Covs"),
  extractor(m2c, "M2: Full Covs + IPW"))
allresults[, Type := factor(Type, levels = Type)]

ggplot(allresults, aes(Type, y = B, ymin = LL, ymax = UL)) +
  geom_pointrange() +
  coord_flip() +
  xlab("") + ylab("Estimate + 95% CI")

3.8 案例研究:具有交互作用的多元线性回归

这个案例研究是模仿一篇期刊文章，研究人员对测试青少年睡眠和消极情绪的认知脆弱性模型感兴趣[5]。在这项研究中，大约 150 名青少年完成了负面情绪(情绪，抑郁和焦虑症状的复合物)、睡眠功能障碍信念(DBAS)、一般功能障碍态度(DAS)、学业压力(压力)和主观睡眠质量(SSQ)的测量，并佩戴了加速计来客观评估睡眠。具体来说，晚上入睡需要多少分钟(睡眠发作潜伏期；SOLacti)。测量是在学校和假期期间收集的，但这里我们关注的是假期。请注意，根据这些数据，主观睡眠质量得分越高，表明睡眠质量越差。

使用 Bei 及其同事[5]的表 2 中的标准化回归系数以及同一篇文章的表 1 中的平均值和标准差，我们可以模拟一个与该文章中使用的数据大致相似的数据集。

set.seed(12345)
adosleep <- data.table(
  SOLacti = rnorm(150, 4.4, 1.3)²,
  DBAS = rnorm(150, 72, 26),
  DAS = rnorm(150, 125, 32),
  Female = rbinom(150, 1, .53),
  Stress = rnorm(150, 32, 11))
adosleep[, SSQ := rnorm(150,
             (.36 * 3 / 12.5) * SOLacti +
             (.16 * 3 / 26) * DBAS +
             (.18 * 3 / .5) * Female +
             (.20 * 3 / 11) * Stress, 2.6)]
adosleep[, MOOD := rnorm(150,
             (-.07 / 12.5) * SOLacti +
             (.29  / 3) * SSQ +
             (.14  / 26) * DBAS +
             (.21  / 32) * DAS +
             (.12  / 32) * SSQ * (DAS-50) +
             (.44  / .5) * Female +
             (.28 / 11) * Stress, 2)]
adosleep[, Female := factor(Female, levels = 0:1,
                            labels = c("Males", "Females"))]

作为一个更大模型的一部分，研究人员假设主观睡眠质量与消极情绪有关，但这种关系会受到一般功能失调态度的调节。特别是，功能失调的态度被认为是一种弱点，因此，功能失调态度高的青少年也认为睡眠质量差，他们更容易受到负面情绪的影响。相反，对于那些功能失调态度水平较低的人，即使主观睡眠质量较差，他们也可能不太容易受到主观睡眠质量和消极情绪之间的影响(关系较弱)。通过活动描记法(一种客观测量)评估的睡眠开始潜伏期、压力、性别和对睡眠的功能障碍信念也包括在分析中。

首先，我们检查核心变量，以检查它们的分布并寻找任何异常值(图 3-13 )。

img/439480_1_En_3_Fig13_HTML.png

图 3-13

案例研究变量的分布

plot_grid(
  testdistr(adosleep$MOOD, extremevalues = "theoretical",
            plot=FALSE, varlab = "MOOD")$Density,
  testdistr(adosleep$SSQ, extremevalues = "theoretical",
            plot=FALSE, varlab = "SSQ")$Density,
  testdistr(adosleep$SOLacti, extremevalues = "theoretical",
            plot=FALSE, varlab = "SOLacti")$Density,
  testdistr(adosleep$DAS, extremevalues = "theoretical",
            plot=FALSE, varlab = "DAS")$Density,
  ncol = 2)

接下来，我们检查(连续)研究变量之间的二元相关性(图 3-14 )。

img/439480_1_En_3_Fig14_HTML.png

图 3-14

研究变量之间相关性的热图

plot(SEMSummary(
  ~ MOOD + SOLacti + DBAS + DAS + Stress + SSQ,
  data = adosleep), plot = "cor") +
  theme(axis.text.x = element_text(
          angle = 45, hjust = 1, vjust = 1))

接下来，我们为这项研究制作一个描述性统计表。尽管非常基本，但在大多数结果展示中，研究变量的描述性统计表格或图表是标准做法，以使读者更好地理解每个被测试变量的分布和范围。这里，我们利用JWileymisc包中的egltable()函数来计算和显示连续变量的平均值和标准偏差，以及离散变量的 N 和百分比(这里只是阴性)。

egltable(c("SOLacti", "SSQ", "MOOD", "Stress",
           "DBAS", "DAS", "Female"),
         data = as.data.frame(adosleep))

##                M (SD)/N (%)
## 1:   SOLacti  23.33 (13.60)
## 2:       SSQ    6.18 (3.00)
## 3:      MOOD    4.53 (2.49)
## 4:    Stress  32.84 (10.92)
## 5:      DBAS  72.10 (23.88)
## 6:       DAS 130.57 (30.45)
## 7:    Female
## 8:     Males      67 (44.7)
## 9:   Females      83 (55.3)

为了得到 Bei 及其同事论文[5]中的标准化估计值，我们可以标准化预测值。

adosleep[, zMOOD := as.vector(scale(MOOD))]
adosleep[, zDBAS := as.vector(scale(DBAS))]
adosleep[, zDAS := as.vector(scale(DAS))]
adosleep[, zSSQ := as.vector(scale(SSQ))]
adosleep[, zSOLacti := as.vector(scale(SOLacti))]
adosleep[, zStress := as.vector(scale(Stress))]

接下来，我们拟合三个不同的模型进行比较。首先我们从所有协变量开始。第二，我们添加没有交互作用的感兴趣的主要结构。第三，我们添加了主观睡眠质量和全球功能失调信念之间的假设相互作用。最后，我们使用screenreg()函数将所有结果放在一个表中。在这种情况下，screenreg()输出圆括号中的系数和标准误差，p 值阈值用星号表示。这种布局使得比较系数如何根据模型中的其他因素而变化变得非常容易。

m.adosleep1 <- ols(zMOOD ~ zSOLacti + zDBAS + Female + zStress,
                   data = adosleep)
m.adosleep2 <- update(m.adosleep1, . ~ . + zSSQ + zDAS)
m.adosleep3 <- update(m.adosleep2, . ~ . + zSSQ:zDAS)

screenreg(list(m.adosleep1, m.adosleep2, m.adosleep3))

##
## ==================================================
##                 Model 1     Model 2     Model 3
## --------------------------------------------------
## Intercept        -0.24 *     -0.28 **    -0.28 **
##                  (0.11)      (0.09)      (0.09)
## zSOLacti          0.17 *      0.04        0.03
##                  (0.08)      (0.07)      (0.07)
## zDBAS             0.14        0.07        0.08
##                  (0.08)      (0.06)      (0.06)
## Female=Females    0.44 **     0.50 ***    0.50 ***
##                  (0.15)      (0.13)      (0.13)
## zStress           0.26 ***    0.19 **     0.20 **
##                  (0.08)      (0.07)      (0.07)
## zSSQ                          0.34 ***    0.34 ***
##                              (0.07)      (0.07)
## zDAS                          0.41 ***    0.44 ***
##                              (0.06)      (0.06)
## zSSQ * zDAS                               0.14 *
##                                          (0.07)
## --------------------------------------------------
## Num. obs.       150         150         150
## R²               0.18        0.43        0.45
## Adj. R²          0.16        0.41        0.42
## L.R.             29.54       85.59       89.91
## ==================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05

回想一下，主观睡眠质量得分越高，表明睡眠质量越差，模型 2 显示，总体睡眠质量差和总体功能失调的态度与更多的负面情绪显著相关(p < .001)。转到模型 3，睡眠质量和功能失调态度之间的相互作用是积极而显著的，表明青少年的功能失调态度越多，主观睡眠质量差和消极情绪之间的关系就越强。

为了确保模型看起来合适，我们可以检查方差膨胀因子和残差的分布。合理的起点是最复杂的模型，因为它具有最多的变量，其中一些变量最有可能是共线的。

vif(m.adosleep3)

##       zSOLacti          zDBAS Female=Females        zStress
##            1.2            1.1            1.0            1.1
##           zSSQ           zDAS    zSSQ * zDAS
##            1.3            1.1            1.1

testdistr(resid(m.adosleep3), plot=FALSE, varlab = "Residuals")$QQPlot

没有一个方差膨胀因子特别高，残差呈现正态分布(图 3-15 )，所以我们继续。

img/439480_1_En_3_Fig15_HTML.png

图 3-15

模型残差的分布

为了解开交互，创建一个图是有帮助的。出于图示的目的，如果测量的尺度有意义(例如，以年为单位的年龄)，或者该测量常用于读者可能熟悉的原始尺度，则通常使用变量的原始尺度。我们还创建了一个新的数据集来生成图表预测。我们持有的所有协变量的均值或众数。作为“低”和“高”功能失调态度的示例值，我们使用平均值-一个标准偏差和平均值+一个标准偏差。

## refit model on raw data
m.adosleep.raw <- ols(MOOD ~ SOLacti + DBAS + Female +
                        Stress + SSQ * DAS,
                      data = adosleep)

## create a dataset
adosleep.newdat <- as.data.table(with(adosleep, expand.grid(
  SOLacti = mean(SOLacti),
  DBAS = mean(DBAS),
  Female = factor("Females", levels(Female)),
  Stress = mean(Stress),
  SSQ = seq(from = min(SSQ), to = max(SSQ), length.out = 100),
  DAS = mean(DAS) + c(1, -1) * sd(DAS))))

adosleep.newdat$MOOD <- predict(m.adosleep.raw,
                                newdata = adosleep.newdat,
                                se.fit = FALSE)

adosleep.newdat[, DAS := factor(round(DAS),
  levels = c(100, 161),
  labels = c("M - 1 SD", "M + 1 SD"))]

图 3-16 显示了主观睡眠质量和消极情绪之间的关系如何在低水平和高水平的功能失调态度中变化，对于具有较高水平功能失调态度的脆弱青少年，不良睡眠质量和消极情绪之间的关系被夸大了。

img/439480_1_En_3_Fig16_HTML.png

图 3-16

主观睡眠质量和预测消极情绪的整体功能失调信念之间的相互作用

ggplot(adosleep.newdat, aes(SSQ, MOOD, linetype=DAS)) +
  geom_line(size = 2) +
  scale_x_continuous("Subjective sleep quality\n(higher is worse)") +
  ylab("Negative Mood") +
  theme_cowplot() +
  theme(
    legend.position = c(.85, .15),
    legend.key.width = unit(2, "cm"))

总的来说，这为我们从这个案例研究开始的问题提供了一个相当彻底的检查。关于呈现结果的其他方式以及更详细的解释可能是什么样的想法，我们请读者参考在线免费提供的原始文章[5]:doi . org/10 . 5665/sleep . 4508。

3.9 摘要

本章介绍了广义线性模型(GLM)的概念背景，以及这个广泛的框架如何将许多常见的统计模型作为特例。还介绍了 GLM 的两个具体特例:方差分析和线性回归，当焦点结果或因变量是连续的且呈正态分布时，这两种方法是合适的。表 3-10 总结了本章介绍的一些主要建模功能。下一章将继续 GLM 主题，但重点是不遵循正态分布的结果的 GLMs，如二进制和计数数据。

表 3-10

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

功能

它的作用

| | --- | --- | | model.matrix() | 采用描述解释性预测变量的公式，并生成设计矩阵以适应 GLMs。 | | update() | 便于更新现有公式或模型对象，以添加或删除变量。 | | ezANOVA() | 来自ez包的函数，该包提供了一个框架，用于拟合多种类型的方差分析模型以及拟合指数和诊断。 | | anova() | 拟合方差分析模型的内置函数。也用于从 GLM 模型中获取方差分析汇总表。 | | ols() | 符合线性回归模型的rms包中的函数，以及全面的默认输出和诊断信息。 | | glm() | 用于拟合广义线性模型的内置R函数，包括线性回归等等。 | | summary() | 经常对拟合 GLM 的结果调用该函数，以产生额外的汇总信息。 | | coef() | 从线性回归或其他拟合 GLMs 中提取回归系数的函数。 | | vif() | 计算线性回归模型中每个预测值的方差膨胀因子的函数，用于确定由于模型中存在共线解释变量而导致参数协方差矩阵膨胀的程度。 | | predict() | 一种使用现有拟合模型对新数据生成预测值的功能。也可用于从模型生成交互图或预测图。 | | ipwpoint() | 计算离散或连续预测值/暴露变量在一个时间点的逆概率权重的函数。适用于计算边际结构模型和尝试说明潜在的混淆因素。 | | texreg() | 从各种模型制作表格的功能。也可以在一个表格中包含多个型号，以便于比较。 |****