R-统计编程和数据建模高级教程-二-

布客飞龙
85 阅读11分钟

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

原文：Advanced R Statistical Programming and Data Models

协议：CC BY-NC-SA 4.0

四、GLM 2

广义线性模型(GLMs)也可以适应非连续和正态分布的结果。事实上，GLMs 的一大优势是它们提供了一个统一的框架来理解应用于假设来自各种分布的变量的回归模型。对于这一章，我们将主要依靠一个优秀的R包VGAM,，它为向量广义线性模型(VGLMs)和向量广义加法模型(VGAMs)提供了实用程序【125】。VGLMs 和 VGAMs 是一类更加灵活的模型，其中可能有多个响应。然而，除了提供多参数的灵活性，VGAM包实现了超过 20 个链接函数，超过 50 个不同的模型/假设分布。在这一章中，我们将只触及VGAM包功能的表面，但是它的巨大灵活性意味着我们将不需要引入许多不同的包，也不需要许多不同的功能。如果您想更深入地了解 VGLMs 和 VGAMs，我们推荐一本由VGAM软件包[125]的作者撰写的优秀书籍。

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(ggthemes)
library(scales)
library(viridis)
library(VGAM)
library(ipw)
library(JWileymisc)
library(xtable)
library(texreg)

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

4.1 概念背景

本章涵盖了几种特定类型的 GLMs，它们共同构成了除最常见的线性模型之外的大多数其他常用 GLMs。

逻辑回归

离散数据的三种 GLMs 是二元、有序和多项式逻辑回归。它们共享规范链接函数 logit，其一般形式为:

{\log}_e\left(\frac{p}{1-p}\right)

(4.1)

其中 p 是某个概率。对于二元和多元逻辑回归，我们可以将概率写成

{\log}_e\left(\frac{P\left(Y=j|X\right)}{P\left(Y=M+1|X\right)}\right)

(4.2)

对于j1、...，M 这里有 M + 1 个离散层次的结果。在二元的情况下，有两个级别的结果， M = 1 和 M + 1 = 2，因为我们知道属于任何给定组(结果级别)的概率之和必须是 1:

\sum \limits_{j=1}^{M+1}P\left(Y=j|X\right)=1

(4.3)

那么二元结果的 logit 表达式简化为

{\log}_e\left(\frac{P\left(Y=1|X\right)}{1-P\left(Y=1|X\right)}\right)

(4.4)

对于有序逻辑回归情况，标准做法是使用累积逻辑，其形式为

{\log}_e\left(\frac{P\left(Y\le j|X\right)}{1-P\left(Y\le j|X\right)}\right)

(4.5)

概率的比率被称为几率:

Odds=\frac{P\left(Y=1|X\right)}{1-P\left(Y=1|X\right)}

(4.6)

例如，假设给定一些预测值和模型，P(Y= 1 |X)= 0.75，那么

Odds=\frac{.75}{1-.75}=\frac{.75}{.25}=3

(4.7)

如果赔率 = 3被解释为具有该组特定预测值的某人预期发生该事件( Y = 1)的可能性是不发生该事件( Y = 0)的三倍。logit 是概率的对数，确保至少在理论上有可能从∞到+∞。

回归系数也基于赔率。以具有单个预测器的模型的最简单情况为例， x 1 该模型将为

{\log}_e\left(\frac{P\left({Y}_i=1|X={x}_{1i}\right)}{1-P\left({Y}_i=1|X={x}_{1i}\right)}\right)={\beta}_0+{\beta}_1\ast {x}_{1i}

(4.8)

系数 β 1 则定义为

{\beta}_1={\log}_e\left(\frac{P\left({Y}_i=1|X={x}_{1i}+1\right)}{P\left({Y}_i=1|X={x}_{1i}\right)}\right)

(4.9)

这是给定的赔率比的自然对数xI1vsxI1+1。通常的做法是报告比值比而不是对数比值比，这是通过对系数求幂来实现的。

OddsRatio= OR={e}^{\beta_1}

(4.10)

假设β1= 0.5，那么e0*。* 5 = 1.65，我们将此解释为表明 x 1 中的一个单位变化与发生该事件的 1.65 倍几率相关。与具有相当自然的解释的优势相比，优势比在某种程度上更难解释，因为它们代表了优势的倍增变化，但没有告诉我们优势从何开始。为了说明这一点，2 的优势比同样可以从. 02/.01 和 2/1 中产生。在这两种情况下，赔率都是 2，但即使确实有人有两倍的赔率，在绝对基础上，赔率仍然很小，0.02，在另一种情况下，新的赔率表明某人发生该事件的可能性是不发生该事件的两倍。

计数回归

我们将在本章介绍的其他类型的 GLMs 是为计数数据设计的模型。计数数据是离散的，类似于逻辑回归中使用的结果，但与逻辑回归不同，计数结果可以采用许多值，并且是有序的。形式上，计数结果的范围是 0∞中的整数。计数结果出现在各种情况下。例如，在医疗环境中，特别是随着人口老龄化，越来越常见的是考虑共病的计数。在保险行业，建立一个人会遭遇多少次事故的模型是可取的。很多人可能是零事故，但有的会有一次，有的两次，更少的三次，四次等。在生产中，故障率和失败率很重要。如果一条生产线或一家工厂有大量不合格产品，这对于减少错误、降低成本和提高质量控制来说是非常有价值的信息。

计数数据的两种最常见的 GLMs 类型是泊松回归和负二项式回归。泊松分布只有一个参数，即比率或平均值，通常表示为 λ 。给定该参数，泊松的概率质量函数为

P\left(Y=y;\lambda \right)=\frac{e^{-\lambda }{\lambda}^y}{y!}

(4.11)

泊松分布的均值和方差都是 λ 。为了比较，这里有两个不同比率的泊松分布，如图 4-1 所示。注意，当我们有预先计算的值并且想要一个条形图时，我们使用geom_col()而不是geom_bar()。

img/439480_1_En_4_Fig1_HTML.png

图 4-1

λ= 2 和λ= 6 的泊松分布的密度

dpoisson <- data.table(X = 0:20)
dpoisson[, Lambda2 := dpois(X, lambda = 2)]
dpoisson[, Lambda6 := dpois(X, lambda = 6)]

ggplot(melt(dpoisson, id.vars = "X"),
       aes(X, value, fill = variable)) +
  geom_col(position = "dodge") +
  scale_fill_viridis(discrete = TRUE) +
  theme(legend.position = c(.7, .8)) +
  xlab("Y Score") + ylab("Poisson Density")

泊松回归和负二项式回归的标准关联函数都是自然对数。因此，泊松和负二项式回归的标准模型是

{\log}_e\left({Y}_i\right)={\beta}_0+{\beta}_1\ast {x}_1+\cdots +{\beta}_k\ast {x}_k

(4.12)

这给了他们的系数一个方便的解释。例如， β 1 将是 x 1 中的一个单位变化事件预计会发生多少次。

负二项式回归与泊松回归非常相似。唯一的变化是负二项分布包括一个额外的参数，允许方差不同于均值。负二项分布的概率质量函数比泊松分布的更复杂:

P\left(Y=y;\lambda, v\right)=\left(\begin{array}{c}y+v-1\\ {}y\end{array}\right){\left(\frac{\lambda }{\lambda +v}\right)}^y{\left(\frac{v}{v+\lambda}\right)}^v

(4.13)

负二项分布的均值仍然是 λ 。然而，方差现在由\lambda +\frac{\lambda²}{v}.{v}^{-1}给出，它被称为尺度或分散参数，有时也被称为辅助参数。随着 v 的增加，负二项分布变得越来越接近泊松，这样负二项 limv→∞就是泊松。

在许多应用案例中，负二项式是比泊松更好的选择，因为它只是稍微复杂一些，而且往往是更真实的数据匹配，因为泊松分布要求的均值和方差相同的假设经常被违反。

负二项式回归的解释或多或少与泊松回归相同，所以不需要额外的努力。

4.2 R示例

对于本章中的例子，我们将再次使用美国人的改变生活[45]研究数据。“简介”中标有“数据设置”的部分介绍了数据的读取和准备从技术上讲，数据具有采样权重，但为了简单起见，我们忽略这些权重。没有加权，分析仍然是正确的；它们只是不能反映抽样人口。

acl <- readRDS("advancedr_acl_data.RDS")

二元逻辑回归

要尝试二元逻辑回归模型，我们需要一个二元结果变量。我们通过将吸烟状态转换为当前吸烟者与非吸烟者(以前或从不吸烟)的二元关系来实现这一点。然后，我们使用带有family = binomialff()和 logit 链接的vglm()函数来运行我们的逻辑回归模型。与其他模型一样，summary()函数提供了模型和系数的摘要。

acl$CurSmoke <- as.integer(acl$Smoke_W1 == "(1) Cur Smok")

m.lr <- vglm(CurSmoke ~ Sex,
             family = binomialff(link = "logit"),
             data = acl, model = TRUE)
summary(m.lr)

##
## Call:
## vglm(formula = CurSmoke ~ Sex, family = binomialff(link = "logit"),
##     data = acl, model = TRUE)
##
##
## Pearson residuals:
##                Min   1Q Median  3Q Max
## logit(prob) -0.712 -0.603 -0.603 1.4 1.66
##
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)    -0.6788     0.0574  -11.82  < 2e-16 ***
## Sex(2) FEMALE  -0.3314     0.0746  -4.44   8.8e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1
##
## Number of linear predictors: 1
##
## Name of linear predictor: logit(prob)
##
## Residual deviance: 4356 on 3615 degrees of freedom
##
## Log-likelihood: -2178 on 3615 degrees of freedom
##
## Number of iterations: 4
##
## No Hauck-Donner effect found in any of the estimates

在这个简单的例子中，有一个二元结果和一个二元预测值，我们可以很容易地使用频率直接计算优势比。在 2×2 频率表中，比值比是比值比，逻辑回归报告的系数是比值比的自然对数。具体来说，使用表 4-1 中所示的表格，

优势比定义为

\frac{\frac{a}{c}}{\frac{b}{d}}

表 4-1

假设频率表

|

预言者

|

没有烟

|

烟

| | --- | --- | --- | | 男性的 | A | B | | 女性的 | C | D |

对于我们的数据，实际频率表在表 4-2 中:

表 4-2

观察频率表

|   |

Zero

|

one

| | --- | --- | --- | | (1)男性 | Nine hundred and one | Four hundred and fifty-seven | | (2)女性 | One thousand six hundred and fifty-six | Six hundred and three |

or.tab <- xtabs(~ Sex + CurSmoke, data = acl)
or.tab.res <- (or.tab[1,1]/or.tab[2,1])/(or.tab[1,2]/or.tab[2,2])
xtable(or.tab, caption = "Observed frequency table",
       label = "tglm2-obsfreq")

由此产生的优势比可以计算为

\frac{\frac{901}{1656}}{\frac{457}{603}}=0.72

几率的自然对数是-0.33，这与我们的逻辑回归模型的系数-0.33 相同。

理解比值比是如何计算的，有助于理解如何解释它们。2 x 2 频率表中比值比的简单公式

\frac{\frac{a}{c}}{\frac{b}{d}}

也解释了逻辑回归的要求，即不能有任何零单元格。如果任何单元格为零，如果 c 、 b 或 d 为零，则结果是未定义的(除以零)，或者如果 a 为零，则结果正好为零，并且零的对数是负无穷大，这意味着逻辑回归的系数将是负无穷大，这也是有问题的。

尽管能够阅读方程来计算比值比或对数比值比，但对大多数人来说，这不是一个直观的值来解释。许多人发现概率尺度更容易解释。例如，我们可以找到男性和女性吸烟的概率(比例或百分比)。回到 4-1，男性吸烟概率的等式是

\frac{b}{a+b}

女性吸烟概率的等式是

\frac{d}{c+d}

我们还可以报告两种概率的差异，以量化男性和女性吸烟的概率有多大差异。获得概率的最简单和最通用的方法是基于数据和我们的模型来预测它们。我们可以使用predict()函数。为了在概率尺度上而不是在对数概率尺度上获得预测，我们使用可选参数type = "response"，以便将结果转换回原始尺度，如图 4-2 所示。

img/439480_1_En_4_Fig2_HTML.png

图 4-2

按性别显示吸烟概率的图表

preddat <- data.table(Sex = levels(acl$Sex))
preddat$yhat <- predict(m.lr, newdata = preddat,
        type = "response")

ggplot(preddat, aes(Sex, yhat)) +
  geom_bar(stat = "identity") +
  scale_y_continuous("Smoking Probability", labels = percent) +
  theme_tufte()

我们可以使用xtable函数制作一个漂亮的结果表，结果如表 4-3 所示。

表 4-3

逻辑回归模型的总结，包括系数、标准误差和 p 值

|   |

估计

|

Std。错误

|

z 值

|

公关（>|z|）

| | --- | --- | --- | --- | --- | | (截取) | −0.68 | Zero point zero six | −11.82 | Zero | | 性别(2)女性 | −0.33 | Zero point zero seven | −4.44 | Zero |

xtable(coef(summary(m.lr)), digits = 2,
       caption = paste(
  "Summary of logistic regression model",
  "including coefficients, standard errors",
  "and p-values."), label = "tglm2-orsimple")

我们可以把这个结果解释为女性吸烟的几率比男性低 0.33 倍。我们也可以用比值比来解释这个结果。在这种情况下，我们可以说女性吸烟的几率是男性的 0.72 倍。

作为另一个例子，我们回到上一章所做的检查自我效能。一些潜在的混杂因素是性别、种族/民族和年龄。我们可以计算逆概率权重，并使用这些来调整我们的自我效能和吸烟模型。

## unadjusted model
m0.lr <- vglm(CurSmoke ~ SelfEfficacy_W1,
             family = binomialff(link = "logit"),
             data = acl, model = TRUE)

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

## adjusted logistic regression model
m1.lr <- vglm(CurSmoke ~ SelfEfficacy_W1,
             family = binomialff(link = "logit"),
             data = acl, model = TRUE,
             weights = winsorizor(w$ipw.weights, .01))

然后我们可以使用xtable()函数制作一个表格，比较原始模型和调整后模型的估计值，如表 4-4 所示。在这种情况下，在对性别、种族和年龄进行调整后，结果实际上稍微强一些。

表 4-4

未调整(原始)和调整回归模型的比较

|   |

类型

|

估计

|

Std。错误

|

z 值

|

公关（>|z|）

| | --- | --- | --- | --- | --- | --- | | one | 生的 | −0.88 | Zero point zero four | −24.13 | Zero | | Two | 生的 | −0.06 | Zero point zero three | −1.71 | Zero point zero nine | | three | 形容词 | −0.88 | Zero point zero four | −24.12 | Zero | | four | 形容词 | −0.08 | Zero point zero three | −2.36 | Zero point zero two |

xtable(rbind(
  data.table(Type = "Raw", coef(summary(m0.lr))),
  data.table(Type = "Adj", coef(summary(m1.lr)))),
  digits = 2,
  caption = paste("Comparison of unadjusted (raw)",
    "and adjusted regression models"),
  label = "tglm2-lrcompare")

我们可以将结果解释为表明自我效能增加一个单位与吸烟的对数几率降低-0.08 相关。我们也可以用比值比来解释这个结果。在这种情况下，我们会说，自我效能感每增加一个单位，吸烟的几率就会增加 0.92 倍。

通常，因为比值比没有绝对的解释，所以用绝对概率来表示结果是有帮助的。我们可以再次依赖于predict()函数来生成特定的预测概率并绘制它们。为了做到这一点，我们可能会产生一系列自我效能值的预测。图 4-3 中的结果显示，随着自我效能感的提高，当前吸烟者的概率下降。在这种情况下，即使在概率尺度上，它也是近似线性的，因为没有概率接近 0 或 1。

img/439480_1_En_4_Fig3_HTML.png

图 4-3

根据自我效能显示吸烟概率的图表

preddat2 <- data.table(SelfEfficacy_W1 =
  seq(from = min(acl$SelfEfficacy_W1, na.rm = TRUE),
      to = max(acl$SelfEfficacy_W1, na.rm = TRUE),
      length.out = 1000))
preddat2$yhat <- predict(m1.lr, newdata = preddat2,
                        type = "response")

ggplot(preddat2, aes(SelfEfficacy_W1, yhat)) +
  geom_line() +
  scale_x_continuous("Self-Efficacy") +
  scale_y_continuous("Smoking Probability", label = percent) +
  theme_tufte() + coord_cartesian(ylim = c(.25, .40))

最后，有时人们会根据数据集计算概率的平均变化。因为在概率尺度上，结果不是线性的(尽管在图 4-3 中它们是近似线性的)，与自我效能感变化相关的概率变化取决于一个人的初始自我效能感水平。此外，如果模型中有其他变量，变化也将取决于这些其他变量。处理这个问题的一种方法是使用实际数据集生成预测概率，然后使用实际数据集，但稍微增加每个人的自我效能。在这种情况下，我们不必对人们的自我效能或其他预测分数做出任何不切实际的假设，我们使用他们的实际分数。然后我们可以找到每个人吸烟概率的预测变化，最后对所有这些进行平均，得到概率的平均边际变化。

## delta value for change in self efficacy
delta <- .01

## create a copy of the dataset
## where we increase everyone's self-efficacy by delta
aclalt <- copy(acl)
aclalt$SelfEfficacy_W1 <- aclalt$SelfEfficacy_W1 + delta

## calculate predicted probabilities
p1 <- predict(m1.lr, newdata = acl, type = "response")
p2 <- predict(m1.lr, newdata = aclalt, type = "response")

## calculate the average, marginal change in probabilities
## per unit change in self efficacy
## in percents and rounded
round(mean((p2 - p1) / delta) * 100, 1)

## [1] -1.7

代码显示，在这个样本中，平均边际效应是这样的，自我效能增加一个单位，预计会导致当前吸烟的机会降低 1.7%。

有序逻辑回归

当一个结果是离散的和分类的时，有序逻辑回归是有用的，但在分类有一个自然的顺序时也是有用的。在 ACL 数据集中，身体活动分为五类，从最不活跃到最活跃。

acl$PhysActCat_W2 <- factor(acl$PhysActCat_W2, ordered = TRUE)

## adjusted ordered logistic regression model
m0.or <- vglm(PhysActCat_W2 ~ SelfEfficacy_W1,
              family = propodds(),
              data = acl)

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

## adjusted ordered logistic regression model
m1.or <- vglm(PhysActCat_W2 ~ SelfEfficacy_W1,
             family = propodds(),
             data = acl, model = TRUE,
             weights = winsorizor(w$ipw.weights, .01))

在有序逻辑回归模型中，有多个截距，比结果中唯一级别的数量少一个。因为我们使用了比例优势模型，假设自我效能与结果的关联在所有水平上都是成比例的，所以自我效能只估计了一个系数。表 4-5 显示了原始模型、未调整模型和调整模型的比较。

表 4-5

未调整(原始)和调整有序逻辑回归模型的比较

|   |

类型

|

标签

|

估计

|

Std。错误

|

z 值

|

公关（>|z|）

| | --- | --- | --- | --- | --- | --- | --- | | one | 生的 | (截距):1 | Zero point eight | Zero point zero four | Nineteen point seven | Zero | | Two | 生的 | (截距):2 | Zero point zero eight | Zero point zero four | Two point zero three | Zero point zero four | | three | 生的 | (截距):3 | −1.14 | Zero point zero four | −26.11 | Zero | | four | 生的 | (截距):4 | −1.88 | Zero point zero six | −34.20 | Zero | | five | 生的 | 自我效能 _W1 | Zero point two two | Zero point zero three | Six point six | Zero | | six | 形容词 | (截距):1 | Zero point seven nine | Zero point zero four | Nineteen point five six | Zero | | seven | 形容词 | (截距):2 | Zero point zero seven | Zero point zero four | One point nine three | Zero point zero five | | eight | 形容词 | (截距):3 | −1.14 | Zero point zero four | −26.15 | Zero | | nine | 形容词 | (截距):4 | −1.89 | Zero point zero six | −34.22 | Zero | | Ten | 形容词 | 自我效能 _W1 | Zero point one nine | Zero point zero three | Five point eight seven | Zero |

xtable(rbind(
  data.table(Type = "Raw",
             Labels = rownames(coef(summary(m0.or))),
             coef(summary(m0.or))),
  data.table(Type = "Adj",
             Labels = rownames(coef(summary(m1.or))),
             coef(summary(m1.or)))),
  digits = 2,
  caption = paste("Comparison of unadjusted (raw) and",
   "adjusted ordered logistic regression models"),
  label = "tglm2-orcompare")

与二元逻辑回归一样，绘制预测概率有助于提供更直接的数据解释视图。因为有多个类别，所以我们最终得到多个概率，然后使用melt()函数将它们合并成一个长数据集。结果图如图 4-4 所示。我们可以看到，随着自我效能的提高，成为前四类成员的概率有适度的增加，但被成为最低类成员的概率的急剧下降所抵消。

img/439480_1_En_4_Fig4_HTML.png

图 4-4

图表按自我效能显示不同身体活动类别的概率

preddat3 <- data.table(SelfEfficacy_W1 =
  seq(from = min(acl$SelfEfficacy_W1, na.rm = TRUE),
      to = max(acl$SelfEfficacy_W1, na.rm = TRUE),
      length.out = 1000))
preddat3 <- cbind(preddat3,
  predict(m1.or, newdata = preddat3,
          type = "response"))
preddat3 <- melt(preddat3, id.vars = "SelfEfficacy_W1")

ggplot(preddat3, aes(SelfEfficacy_W1, value,
                     colour = variable, linetype = variable)) +
  geom_line(size = 2) +
  scale_color_viridis(discrete = TRUE) +
  scale_x_continuous("Self-Efficacy") +
  scale_y_continuous("Activity Probability", label = percent) +
  coord_cartesian(ylim = c(0, .6), expand = FALSE) +
  theme_tufte() +
  theme(legend.position = c(.7, .8),
        legend.key.width = unit(2, "cm"))

与二元逻辑回归一样，我们可以计算自我效能单位变化的预测概率的平均边际变化。对于多个类别，我们得到每个类别的平均边际变化。

## delta value for change in self efficacy
delta <- .01

## create a copy of the dataset
## where we increase everyone's self-efficacy by delta
aclalt <- copy(acl)
aclalt$SelfEfficacy_W1 <- aclalt$SelfEfficacy_W1 + delta

## calculate predicted probabilities
p1 <- predict(m1.or, newdata = acl, type = "response")
p2 <- predict(m1.or, newdata = aclalt, type = "response")

## average marginal change in probability of
## membership in each category
round(colMeans((p2 - p1) / delta) * 100, 1)

## (1) Low_5th (2) 2Low_5th (3) 3Low_5th (4) 4Low_5th (5) Hi_5th
##        -4.2         -0.6          1.3          1.3        2.2

多项式逻辑回归

多项逻辑回归类似于有序逻辑回归，因为它适用于结果变量有两个以上的水平。然而，与假设比例优势的有序逻辑回归不同，多项逻辑回归模型不假设任何比例效应或水平排序。然而，这种灵活性是以大量参数和解释结果的复杂性增加为代价的。

为此，我们将查看 ACL 数据中的雇佣信息。ACL 根据员工的工作时数对他们进行编码。为简单起见，我们将把它归为一个单独的就业类别。

acl[, EmployG_W2 := as.character(Employment_W2)]
acl[EmployG_W2 %in% c(
  "(2) 2500+HRS", "(3) 15002499",
  "(4) 500-1499", "(5) 1-499HRS"),
  EmployG_W2 := "(2) EMPLOYED"]
acl[, EmployG_W2 := factor(EmployG_W2)]

重新编码后，得到的频率表如表 4-6 所示。

表 4-6

就业频率表

|   |

Var1

|

频率

| | --- | --- | --- | | one | (1)残疾人 | One hundred and twenty-two | | Two | (2)就业 | One thousand four hundred and seventy-six | | three | (6)退休 | Seven hundred and twenty-four | | four | (7)失业 | Eighty-six | | five | (8)保持 HS | Four hundred and fifty-nine |

xtable(as.data.frame(table(acl$EmployG_W2)),
       caption = "Frequency table of employment",
       label = "tglm2-freqtab")

接下来，我们可以像以前一样使用vglm()函数来估计多项式逻辑回归模型。唯一的变化是我们指定了参数family = multinomial()。正如我们前面的例子一样，我们可以估计未调整和调整后的模型，其中调整后的模型使用 IPWs 来解释性别、种族/民族和年龄的混杂。

## unadjusted multinomial logistic regression model

m0.mr <- vglm(EmployG_W2 ~ SelfEfficacy_W1,
              family = multinomial(),
              data = acl, model = TRUE)

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

## adjusted multinomial logistic regression model

m1.mr <- vglm(EmployG_W2 ~ SelfEfficacy_W1,

              family = multinomial(),
             data = acl, model = TRUE,
             weights = winsorizor(w$ipw.weights, .01))

接下来，我们可以制作一个表格，比较未调整(原始)模型和调整模型的估计值和系数。在多项逻辑回归中，没有假设预测因子的效果在不同的结果水平上是相等的。取而代之的是，k–1 个独立的参数被估计用于每个预测器，其中 k 是结果的唯一级别的数量。考虑多项逻辑回归的另一种方式是，如果选择一个水平的结果作为参考组，那么实际上运行一系列的k–1 二元逻辑回归。唯一真正的变化是，还有一个约束，即属于任何一个群体的概率总和必须为 1，这反映了人们只能属于一个群体，而每个人都必须属于某个群体的现实。来自vglm()的系数用数字标注，这些系数基于因子等级的顺序。结果如表 4-7 所示。

表 4-7

未调整(原始)和调整的多项式逻辑回归模型的比较

|   |

类型

|

标签

|

估计

|

Std。错误

|

z 值

|

Pr( > z)

| | --- | --- | --- | --- | --- | --- | --- | | one | 生的 | (截距):1 | −1.44 | Zero point one one | −12.98 | Zero | | Two | 生的 | (截距):2 | One point one eight | Zero point zero five | Twenty-one point eight three | Zero | | three | 生的 | (截距):3 | Zero point four six | Zero point zero six | Seven point seven three | Zero | | four | 生的 | (截距):4 | −1.72 | Zero point one two | −14.00 | Zero | | five | 生的 | 自我效能 _W1:1 | −0.33 | Zero point zero nine | −3.70 | Zero | | six | 生的 | 自我效能 _W1:2 | Zero point two two | Zero point zero five | Four point two eight | Zero | | seven | 生的 | 自我效能 _W1:3 | Zero point two three | Zero point zero six | Three point nine one | Zero | | eight | 生的 | 自我效能 _W1:4 | −0.17 | Zero point one one | −1.66 | Zero point one | | nine | 形容词 | (截距):1 | −1.44 | Zero point one one | −13.01 | Zero | | Ten | 形容词 | (截距):2 | One point one seven | Zero point zero five | Twenty-one point eight | Zero | | Eleven | 形容词 | (截距):3 | Zero point four six | Zero point zero six | Seven point six five | Zero | | Twelve | 形容词 | (截距):4 | −1.73 | Zero point one two | −14.07 | Zero | | Thirteen | 形容词 | 自我效能 _W1:1 | −0.40 | Zero point zero nine | −4.50 | Zero | | Fourteen | 形容词 | 自我效能 _W1:2 | Zero point one five | Zero point zero five | Two point nine six | Zero | | Fifteen | 形容词 | 自我效能 _W1:3 | Zero point one seven | Zero point zero six | Two point eight two | Zero | | Sixteen | 形容词 | 自我效能 _W1:4 | −0.23 | Zero point one | −2.21 | Zero point zero three |

xtable(rbind(
  data.table(Type = "Raw",
             Labels = rownames(coef(summary(m0.mr))),
             coef(summary(m0.mr))),
  data.table(Type = "Adj",
             Labels = rownames(coef(summary(m1.mr))),
             coef(summary(m1.mr)))),
  digits = 2,
  caption = paste("Comparison of unadjusted (raw) and",
   "adjusted multinomial logistic regression models"),
  label = "tglm2-mrcompare")

我们可以画一个图来显示，在任何特定的就业类别中，预测的概率如何随着自我效能的函数而变化。结果如图 4-5 所示。这个数字告诉我们，随着自我效能的提高，人们不太可能残疾，更有可能就业或退休。它还强调了在这些模型中，随着时间的推移并不总是线性变化的。失能的可能性从自我效能的-4 到-2 迅速下降，然后在自我效能较高的时候下降得更慢。

img/439480_1_En_4_Fig5_HTML.png

图 4-5

图表按自我效能显示不同就业类别的概率

preddat4 <- data.table(SelfEfficacy_W1 =
  seq(from = min(acl$SelfEfficacy_W1, na.rm = TRUE),
      to = max(acl$SelfEfficacy_W1, na.rm = TRUE),
      length.out = 1000))
preddat4 <- cbind(preddat4,
  predict(m1.mr, newdata = preddat4,
          type = "response"))
preddat4 <- melt(preddat4, id.vars = "SelfEfficacy_W1")

ggplot(preddat4, aes(
  SelfEfficacy_W1, value,
  colour = variable, linetype = variable)) +
  geom_line(size = 2) +
  scale_color_viridis(discrete = TRUE) +
  scale_x_continuous("Self-Efficacy") +
  scale_y_continuous("Probability", label = percent) +
  coord_cartesian(ylim = c(0, .65), expand = FALSE) +
  theme_tufte() +
  theme(legend.position = c(.18, .82),
        legend.key.width = unit(2, "cm"))

最后，我们可以计算自我效能的单位变化的预测概率的平均边际变化。对于多个类别，我们得到每个类别的平均边际变化。这些结果表明，平均而言，自我效能感每增加一个单位，就业变化最大(平均增加 2.9%)，其次是残疾(平均减少 2.1%)，其他类别的变化较小。

## delta value for change in self efficacy
delta <- .01

## create a copy of the dataset
## where we increase everyone's self-efficacy by delta
aclalt <- copy(acl)

aclalt$SelfEfficacy_W1 <- aclalt$SelfEfficacy_W1 + delta

## calculate predicted probabilities
p1 <- predict(m1.mr, newdata = acl, type = "response")
p2 <- predict(m1.mr, newdata = aclalt, type = "response")

## average marginal change in probability of
## membership in each category
round(colMeans((p2 - p1) / delta) * 100, 1)

## (1) DISABLED (2) EMPLOYED (6) RETIRED (7) UNEMPLOY (8) KEEP HS
##         -2.1          2.9         1.7         -1.0        -1.5

泊松和负二项式回归

对于计数结果，我们可以使用泊松回归。在 ACL 数据中，一个变量是过去 12 个月中经历的慢性疾病的计数。这种变量可能很适合泊松回归。

我们可能做的第一件事是查看分布，并获得一些基本的描述性统计数据。对于计数结果，平均值和标准偏差可能没有意义，因此中值和四分位间距通常是更好的总结。我们可以使用egltable()函数快速汇总中位数和四分位间距，如下所示。

egltable(c("NChronic12_W1", "NChronic12_W2"),
         data = acl, parametric = FALSE)

##                    Mdn (IQR)
## 1: NChronic12_W1 1.00 (2.00)
## 2: NChronic12_W2 1.00 (2.00)

我们可以使用ggplot()功能来绘制每种慢性疾病的频率柱状图，以更广泛地了解每种波的分布情况。结果如图 4-6 所示。

img/439480_1_En_4_Fig6_HTML.png

图 4-6

图表显示了 ACL 数据中每个波的各种慢性疾病的频率

plot_grid(
  ggplot(acl, aes(NChronic12_W1)) +
  geom_bar() + theme_tufte(),
  ggplot(acl, aes(NChronic12_W2)) +
  geom_bar() + theme_tufte(),
  ncol = 1,
  labels = c("Wave 1", "Wave 2"),
  label_x = .8)

## Warning: Removed 750 rows containing non-finite values (stat_count).

接下来，我们可以使用vglm()函数来估计泊松回归模型。对于泊松回归，我们将族参数指定为family = poissonff()。对summary()的调用提供了模型结果和估计的快速总结。

## unadjusted poisson regression model

m0.pr <- vglm(NChronic12_W2 ~ SelfEfficacy_W1,
              family = poissonff(),
              data = acl, model = TRUE)

summary(m0.pr)

##
## Call:
## vglm(formula = NChronic12_W2 ~ SelfEfficacy_W1, family = poissonff(),
##     data = acl, model = TRUE)
##
##
## Pearson residuals:
##                Min    1Q Median    3Q  Max
## loge(lambda) -1.34 -1.01 -0.117 0.811 5.67
##
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)
## (Intercept)       0.0954     0.0179    5.33  9.8e-08 ***
## SelfEfficacy_W1  -0.1347     0.0165   -8.16  3.3e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1
##
## Number of linear predictors: 1
##
## Name of linear predictor: loge(lambda)
##
## Residual deviance: 4075 on 2865 degrees of freedom
##
## Log-likelihood: -4126 on 2865 degrees of freedom
##
## Number of iterations: 5
##
## No Hauck-Donner effect found in any of the estimates

但是，在继续之前，最好检查一下泊松回归模型的假设是否合理。也就是说，经常会出现过度离差，方差与均值相同的假设被违反。要检验这一假设，最简单的方法是也拟合负二项式回归模型，然后比较两个模型的相对拟合度，以确定负二项式是否提高了拟合度。

为了比较泊松和负二项式结果，我们首先需要拟合负二项式回归模型。唯一需要的改变是将family = poissonff()改为family = negbinomial()。然后，我们可以使用AIC()和BIC返回基于模型似然性的 Akaike 信息标准(AIC)和贝叶斯信息标准(BIC ),这些标准因参数数量而受到惩罚。较低的 AIC 和 BIC 分数表明更好的拟合，即使在考虑模型的复杂性之后。使用 AIC 和 BIC 比使用简单的模型拟合度量更可取，因为通常更复杂的模型提供更好的拟合。我们想知道的是，合身性的提高是否值得增加复杂性；因此，一些参数的惩罚项是需要的，AIC 和 BIC 都包括。

比较 AIC 和 BIC 发现，负二项式回归模型具有较低的 AIC 和较低的 BIC，表明对于这些数据，负二项式模型优于泊松模型。

## unadjusted negative binomial regression model
m0.nbr <- vglm(NChronic12_W2 ~ SelfEfficacy_W1,
              family = negbinomial(),
              data = acl, model = TRUE)

AIC(m0.nbr) - AIC(m0.pr)

## [1] -97

BIC(m0.nbr) - BIC(m0.pr)

## [1] -91

另一个有用的健全性检查是检查模型的模拟值是否与真实的观察数据一致。我们可以使用内置在VGAM包中的simulate()函数轻松实现这一点。它只需要一个模型，但我们也可以指定要生成的模拟的数量，我们只需要一个，并设置随机种子，这样结果是可重复的。接下来，我们使用真实结果分数、泊松模拟和负二项式模型模拟构建一个数据集。最后，我们将所有这些绘制在图 4-7 中。该图告诉我们，我们的两个模型都不能完美地再现真实的分布。然而，我们也可以看到，负二项式模型的模拟比泊松模型的模拟更接近事实。如图 4-7 所示的曲线图对于比较模型和评估模型是否是观测数据的合理近似非常有用。有时候“最好”的模型可能仍然是一个糟糕的模型，我们想要提前知道。

img/439480_1_En_4_Fig7_HTML.png

图 4-7

图表显示了基于真实数据、负二项式模型模拟和泊松回归模型模拟的各种慢性疾病的频率。

test.pr <- simulate(m0.pr, nsim = 1, seed = 1234)$sim_1
test.nbr <- simulate(m0.nbr, nsim = 1, seed = 1234)$sim_1
test.all <- data.table(
  Type = rep(c("Truth", "Poisson", "Negative\nBinomial"),
             times = c(
               nrow(model.frame(m0.pr)),
               length(test.pr),
               length(test.nbr))),
  Score = c(
    model.frame(m0.pr)$NChronic12_W2,
    test.pr,
    test.nbr))

ggplot(test.all, aes(Score, fill = Type)) +
  geom_bar(position = "dodge") +
  scale_fill_viridis(discrete = TRUE) +
  theme_tufte() +
  theme(legend.position = c(.8, .8))

在这一点上，通过比较 AIC 和 BIC 的得分，以及通过可视化泊松和负二项式回归模型的模拟值，很明显我们应该继续使用负二项式模型。如果我们想比较未调整和调整后的结果，我们可以计算 IPWs 并使用这些来估计调整后的模型，考虑性别、种族/民族和年龄的影响。

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

## adjusted negative binomial regression model
m1.nbr <- vglm(NChronic12_W2 ~ SelfEfficacy_W1,
              family = negbinomial(),
             data = acl, model = TRUE,
             weights = winsorizor(w$ipw.weights, .01))

接下来，我们可以制作一个表格，比较未调整(原始)模型和调整模型的估计值和系数。负二项式回归模型包括两个截距，一个是位置截距，称为 μ ，另一个是过度分散参数截距，称为大小参数。对于慢性病平均数的模型，只有第一个截距是相关的。结果如表 4-8 所示。

表 4-8

未经调整(原始)和经过调整的负二项式回归模型的比较

|   |

类型

|

标签

|

估计

|

Std。错误

|

z 值

|

公关（>|z|）

| | --- | --- | --- | --- | --- | --- | --- | | one | 生的 | (截距):1 | Zero point one | Zero point zero two | Four point six five | Zero | | Two | 生的 | (截距):2 | One point two three | Zero point one two | Ten point four one | Zero | | three | 生的 | 自我效能 _W1 | −0.13 | Zero point zero two | −6.98 | Zero | | four | 形容词 | (截距):1 | Zero point one | Zero point zero two | Four point seven three | Zero | | five | 形容词 | (截距):2 | One point two three | Zero point one two | Ten point four one | Zero | | six | 形容词 | 自我效能 _W1 | −0.13 | Zero point zero two | −6.69 | Zero |

xtable(rbind(
  data.table(Type = "Raw",
             Labels = rownames(coef(summary(m0.nbr))),
             coef(summary(m0.nbr))),
  data.table(Type = "Adj",
             Labels = rownames(coef(summary(m1.nbr))),
             coef(summary(m1.nbr)))),
  digits = 2,
  caption = paste("Comparison of unadjusted (raw) and",
   "adjusted negative binomial regression models"),
  label = "tglm2-nbrcompare")

因为泊松和负二项式回归模型都使用自然对数关联函数，所以系数是对数尺度的。如果我们查看表 4-8 中未调整的结果，自我效能系数表明，自我效能增加一个单位与慢性疾病数量变化-0.13 个对数单位相关。或者，我们可以指数化该系数，在这种情况下，解释是自我效能增加一个单位与 0.87 倍的慢性疾病相关。

img/439480_1_En_4_Fig8_HTML.png

图 4-8

图表显示慢性疾病的预测数量作为自我效能的函数

如果我们愿意，我们也可以将预测的平均条件数绘制成自我效能的函数。当生成预测时，与逻辑回归一样，我们指定type = "response"来表示我们希望在原始尺度上进行预测，而不是在链接尺度上，这里是对数变换尺度。结果如图 4-8 所示。这个数字向我们表明，随着自我效能的提高，平均而言，人们预期会有更少的慢性疾病。

preddat5 <- data.table(SelfEfficacy_W1 =
  seq(from = min(acl$SelfEfficacy_W1, na.rm = TRUE),
      to = max(acl$SelfEfficacy_W1, na.rm = TRUE),
      length.out = 1000))
preddat5$yhat <- predict(m1.nbr, newdata = preddat5,
          type = "response")

ggplot(preddat5, aes(SelfEfficacy_W1, yhat)) +
  geom_line() +
  scale_x_continuous("Self-Efficacy") +
  scale_y_continuous("Expected Number Conditions") +
  theme_tufte()

4.3 案例研究:多项逻辑回归

默认情况下，在多项逻辑回归中，参数(例如，优势比)是相对于参照组计算的。虽然这足以说明模型，但在实践中应用多项逻辑回归来考虑所有(或至少关键的)组间成对比较是常见的[14，87]。例如，仅仅知道 B 组和 C 组与 A 组显著不同并不能告诉你 B 组和 C 组是否彼此不同。

同时评估几个预测因子的影响也很常见(例如[14，87])，这需要与评估单个预测因子稍有不同的处理。在本案例研究中，我们将构建一个从提问到最终呈现结果和解释的完整示例。

ACL 数据集包括第一波和第二波的吸烟状态。除了观察一波的吸烟状况，一个有趣的问题是谁随着时间的推移改变了(开始或停止)吸烟，以及哪些因素可能预测这种变化。首先，我们需要为吸烟创造一个新的变量，表明随着时间的变化或稳定。我们如下进行重新编码。结果频率表如表 4-9 所示。

表 4-9

一段时间内吸烟频率表

|   |

Var1

|

频率

| | --- | --- | --- | | one | 稳定从不吸烟 | One thousand two hundred and ninety-two | | Two | 稳定的前吸烟者 | Seven hundred and five | | three | 稳定电流吸烟者 | Six hundred and forty-one | | four | 最近戒烟 | One hundred and sixty-seven | | five | 新烟民 | Sixty-two |

acl[, Smoke_W2W1 := NA_character_]
acl[Smoke_W1 == "(3) Nevr Smo" &
    Smoke_W2 == "(3) W2 Never Smoker",
    Smoke_W2W1 := "Stable Never Smoker"]
acl[Smoke_W1 == "(2) Past Smo" &
    Smoke_W2 == "(2) W2 Former Smoker",
    Smoke_W2W1 := "Stable Former Smoker"]
acl[Smoke_W1 == "(1) Cur Smok" &
    Smoke_W2 == "(1) W2 Current Smoker",
    Smoke_W2W1 := "Stable Current Smoker"]
acl[Smoke_W1 %in% c("(2) Past Smo", "(3) Nevr Smo") &
    Smoke_W2 == "(1) W2 Current Smoker",
    Smoke_W2W1 := "New Smoker"]
acl[Smoke_W1 == "(1) Cur Smok" &
    Smoke_W2 == "(2) W2 Former Smoker",
    Smoke_W2W1 := "Recently Quit Smoker"]

acl[, Smoke_W2W1 := factor(Smoke_W2W1,
  levels = c("Stable Never Smoker", "Stable Former Smoker",
             "Stable Current Smoker", "Recently Quit Smoker",
             "New Smoker"))]

xtable(as.data.frame(table(acl$Smoke_W2W1)),
       caption = "Frequency table of smoking over time",
       label = "tglm2-freqtab-smoke")

在本章的前面，我们只关注了一个预测因子。在现实环境中，我们可能对几个潜在的预测因素感兴趣。一个有趣的问题是，随着时间的推移，社会人口统计学、心理社会或健康类型变量是否是吸烟的更好预测因素。我们将像之前一样使用vglm()函数和family = multinomial()来估计模型，以获得多项式结果。

acl[, SES := as.numeric(SESCategory)]

mr.ses <- vglm(Smoke_W2W1 ~ Sex + SES + AGE_W1,
  family = multinomial(),
  data = acl, model = TRUE)

mr.psych <- vglm(Smoke_W2W1 ~ SWL_W1 + InformalSI_W1 +
  FormalSI_W1 + SelfEfficacy_W1 + CESD11_W1,
  family = multinomial(),
  data = acl, model = TRUE)

mr.health <- vglm(Smoke_W2W1 ~ PhysActCat_W1 +
  BMI_W1 + NChronic12_W1,
  family = multinomial(),
  data = acl, model = TRUE)

我们可以比较每个模型的相对性能，使用 AIC 和 BIC 来惩罚复杂性，如表 4-10 所示。结果表明，社会人口因素是吸烟状况和随时间变化的最佳预测因素。

表 4-10

模型比较

|   |

模型

|

美国化学师学会(American Institute of Chemists)

|

比克

| | --- | --- | --- | --- | | one | 社会人口统计 | Seven thousand and fifty-six point three four | Seven thousand one hundred and fifty-one point seven two | | Two | 社会心理的 | Seven thousand two hundred and three point seven three | Seven thousand three hundred and forty-six point seven nine | | three | 健康 | Seven thousand three hundred and forty point five four | Seven thousand five hundred and seven point four five |

xtable(
  data.table(
  Model = c("Sociodemographics", "Psychosocial", "Health"),
  AIC = c(AIC(mr.ses), AIC(mr.psych), AIC(mr.health)),
  BIC = c(BIC(mr.ses), BIC(mr.psych), BIC(mr.health))),
  caption = "Model Comparisons",
  label = "tglm2-modelcomparisons")

我们可以使用summary()函数检查社会人口统计模型中的各个系数。然而，默认情况下，这些只是与参考水平的比较，默认情况下，VGAM包中的参考水平是最后一个水平，对我们来说是“新烟民”

summary(mr.ses)

##
## Call:
## vglm(formula = Smoke_W2W1 ~ Sex + SES + AGE_W1, family = multinomial(), ##     data = acl, model = TRUE)
##
##
## Pearson residuals:
##                      Min     1Q Median    3Q  Max
## log(mu[,1]/mu[,5]) -7.39 -0.744 -0.412 0.811 1.97
## log(mu[,2]/mu[,5]) -6.95 -0.441 -0.306 -0.186 2.98
## log(mu[,3]/mu[,5]) -6.33 -0.420 -0.289 -0.184 3.58
## log(mu[,4]/mu[,5]) -5.88 -0.202 -0.155 -0.118 4.82
##
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)
## (Intercept):1   -0.79722    0.61458   -1.30  0.19457
## (Intercept):2   -1.16638    0.62965   -1.85  0.06397 .
## (Intercept):3    1.41356    0.61771    2.29  0.02212 *
## (Intercept):4   -1.18639    0.70648   -1.68  0.09309 .
## Sex(2) FEMALE:1  0.76073    0.27200    2.80  0.00516 **
## Sex(2) FEMALE:2 -0.46076    0.27545   -1.67  0.09437 .
## Sex(2) FEMALE:3 -0.04184    0.27459   -0.15  0.87888
## Sex(2) FEMALE:4  0.02589    0.30782    0.08  0.93297
## SES:1            0.51292    0.14821    3.46  0.00054 ***
## SES:2            0.50412    0.15079    3.34  0.00083 ***
## SES:3            0.20550    0.15032    1.37  0.17159
## SES:4            0.39726    0.16691    2.38  0.01731 *
## AGE_W1:1         0.04439    0.00858    5.18  2.3e-07 ***
## AGE_W1:2         0.05430    0.00877    6.19  5.9e-10 ***
## AGE_W1:3         0.01111    0.00871    1.28  0.20181
## AGE_W1:4         0.02732    0.00971    2.82  0.00487 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1
##
## Number of linear predictors:  4
##
## Names of linear predictors:
## log(mu[,1]/mu[,5]), log(mu[,2]/mu[,5]), log(mu[,3]/mu[,5]), log(mu[,4]/mu[,5])
##
## Residual deviance: 7024 on 11452 degrees of freedom
##
## Log-likelihood: -3512 on 11452 degrees of freedom
##
## Number of iterations: 6
##
## No Hauck-Donner effect found in any of the estimates
##
## Reference group is level 5 of the response

在观察这些结果时，我们可以注意到的另一件事是，尽管它在统计上是显著的，但年龄的系数相当小，因为年龄的 1 岁变化是相对较小的变化。我们可以考虑将年龄转换成十岁，这样一个单位的差异更有意义。

acl[, AGE_W1 := AGE_W1 / 10]

如果我们想要明确地导出其他类别之间的对比，我们可以更改参考级别。这可以作为multinomial()函数的选项参数来完成。例如，设置refLevel = 1将使第一个类别成为参考，在我们的例子中是“稳定的从不吸烟者”我们也可以重新运行“稳定的前吸烟者”级别 2 和“稳定的当前吸烟者”级别 3。请注意，从数学上讲，所有这些模型都是相同的，不同的是结果中的默认比较。当重新运行时，我们不需要重新指定整个模型，我们可以使用update()功能简单地更新一个现有的模型。

mr.ses1 <- vglm(Smoke_W2W1 ~ Sex + SES + AGE_W1,
              family = multinomial(refLevel = 1),
              data = acl, model = TRUE)
mr.ses2 <- update(mr.ses1,
                  family = multinomial(refLevel = 2))
mr.ses3 <- update(mr.ses1,
                  family = multinomial(refLevel = 3))

接下来，通常报告比值比而不是对数比值，这是默认输出。报告置信区间也很常见，我们可以使用confint()函数来计算置信区间。我们可以通过组合系数和置信区间来创建一个结果表，在对它们求幂之后，我们就有了比值比和比值比的置信区间。

例如，如果我们观察参照组为“稳定的从不吸烟者”时AGE_W1:1的比值比，它告诉我们，年龄增加 10 岁，成为“稳定的前吸烟者”的几率是“稳定的从不吸烟者”的 1.1 倍

如果我们观察参照组是“稳定的从不吸烟者”时AGE_W1:2的比值比，它告诉我们，年龄增加 10 岁，成为“稳定的当前吸烟者”的几率是“稳定的从不吸烟者”的 0.72 倍

相比之下，如果我们观察参照组是“稳定的当前吸烟者”时AGE_W1:2的比值比，它告诉我们，年龄增加 10 岁，是“稳定的前吸烟者”成为“稳定的当前吸烟者”的 1.54 倍

最后，如果我们观察参照组是“稳定的当前吸烟者”时AGE_W1:3的比值比，它告诉我们年龄增加 10 岁与“最近戒烟者”是“稳定的当前吸烟者”的 1.18 倍相关

data.table(
  Ref = "Stable Never Smoker",
  Term = names(coef(mr.ses1)),
  OR = exp(coef(mr.ses1)),
  exp(confint(mr.ses1)))

##                    Ref            Term   OR 2.5 % 97.5 %
## 1: Stable Never Smoker   (Intercept):1 0.69  0.42   1.15
## 2: Stable Never Smoker   (Intercept):2 9.12  5.55  14.98
## 3: Stable Never Smoker   (Intercept):3 0.68  0.30   1.54
## 4: Stable Never Smoker   (Intercept):4 2.22  0.67   7.40
## 5: Stable Never Smoker Sex(2) FEMALE:1 0.29  0.24   0.36
## 6: Stable Never Smoker Sex(2) FEMALE:2 0.45  0.36   0.55
## 7: Stable Never Smoker Sex(2) FEMALE:3 0.48  0.34   0.67
## 8: Stable Never Smoker Sex(2) FEMALE:4 0.47  0.27   0.80
## 9: Stable Never Smoker           SES:1 0.99  0.90   1.10
## 10: Stable Never Smoker          SES:2 0.74  0.66   0.82
## 11: Stable Never Smoker          SES:3 0.89  0.75   1.06
## 12: Stable Never Smoker          SES:4 0.60  0.45   0.80
## 13: Stable Never Smoker       AGE_W1:1 1.10  1.04   1.17
## 14: Stable Never Smoker       AGE_W1:2 0.72  0.67   0.76
## 15: Stable Never Smoker       AGE_W1:3 0.84  0.76   0.94
## 16: Stable Never Smoker       AGE_W1:4 0.64  0.54   0.76

data.table(
  Ref = "Stable Current Smoker",
  Term = names(coef(mr.ses3)),
  OR = exp(coef(mr.ses3)),
  exp(confint(mr.ses3)))

##                      Ref            Term    OR  2.5 % 97.5 %
## 1: Stable Current Smoker   (Intercept):1 0.110  0.067  0.18
## 2: Stable Current Smoker   (Intercept):2 0.076  0.043  0.13
## 3: Stable Current Smoker   (Intercept):3 0.074  0.032  0.17
## 4: Stable Current Smoker   (Intercept):4 0.243  0.072  0.82
## 5: Stable Current Smoker Sex(2) FEMALE:1 2.231  1.811  2.75
## 6: Stable Current Smoker Sex(2) FEMALE:2 0.658  0.525  0.82
## 7: Stable Current Smoker Sex(2) FEMALE:3 1.070  0.752  1.52
## 8: Stable Current Smoker Sex(2) FEMALE:4 1.043  0.609  1.79
## 9: Stable Current Smoker           SES:1 1.360  1.221  1.51
## 10: Stable Current Smoker          SES:2 1.348  1.195  1.52
## 11: Stable Current Smoker          SES:3 1.211  1.006  1.46
## 12: Stable Current Smoker          SES:4 0.814  0.606  1.09
## 13: Stable Current Smoker       AGE_W1:1 1.395  1.309  1.49
## 14: Stable Current Smoker       AGE_W1:2 1.540  1.432  1.66
## 15: Stable Current Smoker       AGE_W1:3 1.176  1.054  1.31
## 16: Stable Current Smoker       AGE_W1:4 0.895  0.754  1.06

呈现结果的另一种方式是计算预测概率。然而，当有多个预测因子时，这就变得更复杂了，因为我们把其他预测因子放在哪里会影响结果。相反，对于多个预测因子，计算平均边际概率可能是最明智的，它将预测因子保持在它们的观察值，并且一次改变一个预测因子。

## delta value for change in age and SES
delta <- .01

## create a copy of the dataset
## where we increase everyone's age by delta
aclage <- copy(acl)
aclage[, AGE_W1 := AGE_W1 + delta]

## create a copy of the dataset
## where we increase everyone's SES by delta
aclses <- copy(acl)
aclses[, SES := SES + delta]

## create two copies of the data
## one where we set everyone to "female" and another to "male"
aclfemale <- copy(acl)
aclfemale[, Sex := factor("(2) FEMALE",
                          levels = levels(acl$Sex))]

aclmale <- copy(acl)
aclmale[, Sex := factor("(1) MALE",
                        levels = levels(acl$Sex))]

## calculate predicted probabilities
p.ref <- predict(mr.ses1, newdata = acl,
                 type = "response")
p.age <- predict(mr.ses1, newdata = aclage,
                 type = "response")
p.ses <- predict(mr.ses1, newdata = aclses,
                 type = "response")
p.female <- predict(mr.ses1, newdata = aclfemale,
                    type = "response")
p.male <- predict(mr.ses1, newdata = aclmale,
                    type = "response")

最后，我们可以计算预测概率的所有平均边际变化，并将它们放在一个表格中，以便于展示。最终结果如表 4-11 所示。这突出了性的强大影响，女性更有可能是稳定的从不吸烟者。我们还可以看到年龄较大和社会经济地位较高如何与成为稳定的当前吸烟者的概率降低约 5%相关联。

表 4-11

预测概率的平均边际变化

|   |

水平

|

年龄

|

(美)工程科学学会(Society of Engineering Science)

|

女性的

| | --- | --- | --- | --- | --- | | one | 稳定从不吸烟 | Two point eight three | Three point six six | Twenty-three point three four | | Two | 稳定的前吸烟者 | Three point nine four | One point eight five | −16.96 | | three | 稳定电流吸烟者 | −5.46 | −4.52 | −5.06 | | four | 最近戒烟 | −0.55 | −0.12 | −0.93 | | five | 新烟民 | −0.76 | −0.87 | −0.39 |

xtable(
data.table(
  Level = colnames(p.ref),
  Age = colMeans((p.age - p.ref) / delta) * 100,
  SES = colMeans((p.ses - p.ref) / delta) * 100,
  Female = colMeans(p.female - p.male) * 100),
  digits = 2,
  caption = "Average marginal change in predicted probability",
  label = "tglm2-margprobs")

虽然这需要一些额外的工作，但创建这样一个表来显示关键预测者的预测概率的平均边际变化，是一种非常有用的方式，可以用比优势比更直观的格式呈现结果。结合优势比和置信区间来估计不确定性，这提供了一个相对全面的结果介绍。

4.4 总结

本章展示了如何使用广义线性模型(GLMs)来构建离散结果的回归模型，包括二分结果、有序和无序分类结果，以及事件计数或数量结果。虽然这些结果是 GLMs 的一些最常见的用途，但 GLMs 可以适应许多其他类型的结果变量和许多其他分布，而不是这里介绍的分布。优秀的VGAM包支持常见和不常见的分布类型，因此如果您的数据看起来不是正态分布，也不是本章介绍的任何分布，很可能仍然可以使用不同分布的vglm()进行建模。关于VGAM软件包特性的更多细节和报道，请参阅其作者的书。

本章还介绍了一些工具和函数，通过生成预测或在更直接可解释的尺度(如概率)上获得效果来帮助我们解释 GLMs。虽然这样的代码严格来说不是 GLMs 的一部分，但它通常可以使结果更清晰，并帮助分析师和读者理解模型的含义。表 4-12 总结了引入的主要功能及其作用。

表 4-12

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

|

功能

|

它的作用

| | --- | --- | | vglm() | 灵活的函数通过“系列”函数拟合多种类型的特定广义线性模型，提供数百种分布和链接函数组合。 | | binomialff() | VGAM 家族函数将 GLMs 拟合到二项式、离散/分类结果(如吸烟或不吸烟)。 | | propodds() | VGAM 家族函数将 GLMs 拟合到有序、离散/分类结果(例如，身体活动水平)。 | | multinomial() | VGAM 家族函数将 GLMs 拟合到离散/分类结果，没有任何内在顺序的假设(例如，随着时间的推移，不同的预测因子可能对特定结果有或多或少的影响)。 | | poissonff() | VGAM 家族函数将 GLMs 拟合到泊松结果(如慢性疾病)。 | | negbinomial() | VGAM 家族函数拟合 GLMs 来计算方差超过均值的数据(广义泊松)。 | | summary() | 打印对象概要的通用功能，包括vglm型号。 | | coef() | 从模型中提取系数的通用函数，包括vglm模型。 | | confint() | 计算模型置信区间的通用函数，包括vglm模型。 | | update() | 更新现有模型，无需重写未更改的部分(例如，更改吸烟水平的默认参考)。 | | ipwpoint() | 估计逆概率权重。 | | xtable() | 将表格很好地导出到 LATEXor HTML。 | | winsorizor() | 在指定的百分位裁剪异常值。 | | predict() | 采用新的数据预测值，对其应用模型以估计最可能的响应结果，当与type = "response",一起使用时，转换回原始数据的规模。 | | simulate() | 从模型中生成模拟数据，可用于比较vglm模型的模型和原始数据分布。 | | AIC() | 返回基于 Akaike 信息标准的模型似然性(由参数计数决定)。 | | BIC() | 返回基于贝叶斯信息标准的模型似然性(由参数计数决定)。 |

五、通用代数建模系统

广义可加模型(gam)是我们在前面章节中讨论的广义线性模型(glm)的扩展。像 GLMs 一样，gam 适应连续和离散的结果。然而，与完全参数模型的 glm 不同，gam 是半参数模型。GAMs 允许在结果和预测之间混合参数和非参数的关联。对于这一章，我们将主要依靠一个优秀的R包VGAM,，它为向量广义线性模型(VGLMs)和向量广义加法模型(VGAMs)提供了实用程序【125】。VGAMs 是比 gam 更灵活的一类模型，gam 可能有多个响应。然而，除了提供多参数的灵活性，VGAM包实现了超过 20 个链接函数，超过 50 个不同的模型/假设分布。在这一章中，我们将只触及VGAM包功能的表面，但是它的巨大灵活性意味着我们将不需要引入许多不同的包，也不需要许多不同的功能。如果你想更深入地了解 VGAMs，我们推荐一本由VGAM软件包[125]的作者写的优秀书籍。

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(ggthemes)
library(scales)
library(viridis)
library(car)
library(mgcv)
library(VGAM)
library(ipw)
library(JWileymisc)
library(xtable)

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

5.1 概念概述

广义加性模型(gam)是半参数加性模型，它使用非参数平滑函数放松了广义线性模型(GLMs)的线性假设[40]。非参数平滑函数的使用为 gam 提供了极大的灵活性，允许他们对预测值和结果之间的关联进行建模，即使函数形式未知，也不需要函数形式的(正确)规范。这种灵活性使得 gam 被多个学科采用，包括心理学[53]和医学[75]。

下一节介绍平滑样条的概念，这是 gam 的基础。如果你对R中关于现代游戏的优秀、全面的介绍感兴趣，请看【122】。要阅读关于 gam 的估计和推断的基础细节，请参见[57，121，123]。在本章的范围之外，gam 还被扩展到不仅模拟位置，还模拟分布的规模和形状。有关更多信息，请参见描述理论方面的优秀论文[78]和[88]，有关R中的实践方面，请参见GAMLSS包[88]。

平滑样条

gam 背后的一个关键概念是平滑样条，它允许 gam 模拟未知的函数形式。平滑样条部分依赖于多项式。可以证明，一个足够高阶的多项式可以逼近任何函数。然而，通常适当的近似可以采用非常高阶的多项式。虽然理论上是正确的，但实际上通常很难建立一个足够高阶的多项式来近似观察到的数据。特别是，通常情况下，多项式在任一极端都提供非常差的近似。图 5-1 显示了 ACL 数据中年龄与抑郁症状之间关系的单截距、线性、二次、三次、四次和十次多项式示例。与高阶多项式相比，低阶多项式在其极值处产生更适度的预测。即使是二次多项式也会对最年轻和最年长的参与者产生相对极端的预测(图 5-1 )。

img/439480_1_En_5_Fig1_HTML.png

图 5-1

显示仅截距(扁平线)和渐进高阶多项式的图形

acl <- readRDS("advancedr_acl_data.RDS")

ggplot(acl, aes(AGE_W1, CESD11_W1)) +
  stat_smooth(method = "lm", formula = y ˜ 1,
    colour = viridis(6)[1], linetype = 1, se = FALSE) +
  stat_smooth(method = "lm", formula = y ˜ x,
    colour = viridis(6)[2], linetype = 4, se = FALSE) +
  stat_smooth(method = "lm", formula = y ˜ poly(x, 2),
    colour = viridis(6)[3], linetype = 2, se = FALSE) +
  stat_smooth(method = "lm", formula = y ˜ poly(x, 3),
    colour = viridis(6)[4], linetype = 3, se = FALSE) +
  stat_smooth(method = "lm", formula = y ˜ poly(x, 4),
    colour = viridis(6)[5], linetype = 1, se = FALSE) +
  stat_smooth(method = "lm", formula = y ˜ poly(x, 10),
    colour = viridis(6)[6], linetype = 5, se = FALSE)

回归样条试图解决需要非常高阶的多项式来适应不同的函数形式的问题，以及在分布的下端和上端的极端预测的问题。样条最初并不用于统计学。样条最初指的是薄木片，弯曲后在节之间形成平滑的曲线。回归样条使用这种思想，因为它们本质上是分段模型，其中每一段都是多项式模型，并且它们在每一段的边缘(节点)处平滑连接。最简单的样条模型是阶跃函数。为了创建这个步骤，我们用边界点定义逻辑语句，在本例中:x>42和 x* ≤ 65 和 x > 65 和 x ≤ 96。为了在R中达到这些界限，我们使用了一个扩展的逻辑操作符:%gle%。左边应该是一些向量。右边应该是一个长度为 2 的向量，如果左边的数字或向量大于最小值并且小于或等于右边的最大值，它将返回TRUE或FALSE。下面的简单示例演示了这些功能。*

## > and <
1:5 %gl% c(2, 4)

## [1] FALSE FALSE  TRUE FALSE FALSE

## > and <=
1:5 %gle% c(2, 4)

## [1] FALSE FALSE  TRUE  TRUE FALSE

## >= and <
1:5 %gel% c(2, 4)

## [1] FALSE  TRUE  TRUE FALSE FALSE

## >= and <=
1:5 %gele% c(2, 4)

## [1] FALSE  TRUE  TRUE  TRUE FALSE

随着多项式次数的增加，每个结之间会出现线性和二次趋势。图 5-2 显示了两个内部节点和阶跃、线性、二次和三次多项式的结果。

img/439480_1_En_5_Fig2_HTML.png

图 5-2

显示阶跃函数样条、线性样条和二次样条的图形，所有样条都有两个内部节点

ggplot(acl, aes(AGE_W1, CESD11_W1)) +
  stat_smooth(method = "lm",
    formula = y ∼ 1 +
      ifelse(x %gle% c(42, 65), 1, 0) +
      ifelse(x %gle% c(65, 96), 1, 0),
    colour = viridis(6)[1], linetype = 1, se = FALSE) +
  stat_smooth(method = "lm",
    formula = y ∼ bs(x, df = 3, degree = 1L),
    colour = viridis(6)[2], linetype = 2, se = FALSE) +
  stat_smooth(method = "lm",
    formula = y ∼ bs(x, df = 4, degree = 2L),
    colour = viridis(6)[3], linetype = 3, se = FALSE) +
  stat_smooth(method = "lm",
    formula = y ∼ bs(x, df = 5, degree = 3L),
    colour = viridis(6)[4], linetype = 4, se = FALSE)

用于基础的 b 样条或基础样条试图通过在给定区域的基础函数中具有相对小的重叠来工作，这有助于它们在计算上更加稳定，并且有助于它们成为非常受欢迎的样条类型。以下代码通过显示具有固定结的 B 样条曲线的基函数来帮助可视化重叠。因为这些图是相同的，只是它们基于不同的数据集，并且具有不同的标题，我们没有重复代码来绘制这些图(这相当长)，而是将结果存储在一个R对象中，p1。然后我们使用%+%操作符用一个新的数据集替换原始数据。结果如图 5-3 所示。

img/439480_1_En_5_Fig3_HTML.png

图 5-3

显示 B 样条(基本样条)的图形

knots <- c(33, 42, 57, 65, 72)
x <- seq(from = min(acl$AGE_W1),
         to = max(acl$AGE_W1), by = .01)

p1 <- ggplot(melt(bs(x, degree = 1,
          knots = knots, intercept = TRUE)),
          aes(Var1, value, colour = factor(Var2))) +
  geom_line() +
  scale_color_viridis("Basis", discrete = TRUE) +
  theme_tufte()

plot_grid(
  p1 +
    ggtitle("5 Knots, Degree = 1"),
  p1 %+% melt(bs(x, degree = 2,
          knots = knots, intercept = TRUE)) +
    ggtitle("5 Knots, Degree = 2"),
  p1 %+% melt(bs(x, degree = 3,
          knots = knots, intercept = TRUE)) +
    ggtitle("5 Knots, Degree = 3"),
  p1 %+% melt(bs(x, degree = 4,
          knots = knots, intercept = TRUE)) +
    ggtitle("5 Knots, Degree = 4"),
  ncol = 2)

样条的延伸称为光滑样条。平滑样条的基本思想是，我们可以自动学习适当的平滑度，而不是直接指定节点和多项式次数，这需要事先知道需要什么。

自动学习适当平滑度的过程通常通过允许许多结和高度灵活性，并基于某些标准使用惩罚来降低灵活性(施加平滑度)来进行。平滑样条的一种常用方法是广义交叉验证(GCV)标准，或者，如果尺度是已知的(通常是未知的)，本质上是 Akaike 信息标准(AIC)的一种变体，称为无偏风险估计(UBRE)。关于 GCV 的详细情况，见[104]。最后一个选择是使用受限最大似然(REML ),其中平滑分量被视为随机效应，每个平滑有一个“方差分量”。

关于光滑样条的最后一点说明是，通常希望有一些方法来量化它们有多光滑或灵活。这样的值在描述上是有用的，并且在计算值(如模型 AIC)时也起作用，在这种情况下，可能性会受到复杂性的影响。一般的解决方案是使用“有效自由度”或 EDF。根据 EDF 的计算是否包括截距(常数)项，1 或 2 的 EDF 可能对应于线性函数。有时报告 EDF 计算截距，在这种情况下，EDF = 2 对应于线性趋势。然而，如果不计算截距，截距有时被称为有效非线性自由度(ENDF)，则 ENDF = 1 对应于线性趋势。随着 EDF/ENDF 的增加，拟合的灵活性同样增加。关于回归中样条和平滑的更多细节，见[39]。

出于本章的目的，对样条的粗略理解为广义可加模型(gam)提供了基础。gam 通过允许参数(例如，假设一个预测器的线性关联)和非参数(例如，对另一个使用平滑样条)的混合来扩展 GLMs。随着每一项的增加，gam 仍然是可加的，它们的一般形式如下:

g(y)=\eta ={b}_0+{f}_1\left({x}_1\right)+{f}_2\left({x}_2\right)+\cdots +{f}_k\left({x}_k\right)

(5.1)

在这个参数化中，我们有熟悉的截距；然而，代替每个预测值的回归系数，我们有函数， f 1 等。这些函数可以预先指定，例如，iff1(x1)=b1∑x1，或者这些函数可以是平滑样条函数，其中平滑度是基于某种标准从数据中学习的，例如 GCV 或 UBRE。要成为 GAM，通常必须至少有一个平滑项。然而，像 GLMs 一样，gam 在理论上适应许多平滑项和许多常规参数项。混合光滑项和非光滑项的能力使 GAMs 成为一类高度灵活的模型，可以应用于许多环境。例如，平滑项可能具有实质性的意义，如按年龄划分的儿童体重增长图，其中增长通常是非线性的，但函数形式未知。然而，GAMs 也可以应用于其他情况。例如，一个实质性问题可能涉及一个已知或假定参数形式的变量的影响。然而，可能存在关于混杂的担忧，并且混杂的影响具有未知的函数形式。在这种情况下，混杂变量的任何参数本质上都是有害参数，因为唯一的目标是充分模拟它们的(未知)函数形式。没有兴趣真正理解它们的功能形式。在这种情况下,“假设检验”可能是针对具有预先指定的参数形式的变量，从而避免过度拟合的风险，但混杂变量可以使用平滑项稳健地捕获，因为关于它们的统计推断没有什么意义。

除了使用平滑样条，gam 本质上与其他 glm 功能相同。因此，分布假设和可以使用的不同分布族是可以比较的。诚然，由于确定样条适当平滑度的数据驱动性质，自由度和标准误差的计算往往是不同的，可能更接近 gam。下一节将介绍R中的 gam，包括如何估计它们，绘制结果，以及呈现或使用估计的模型。

5.2 游戏在R

高斯结果

基本游戏

gam 可以使用vgam()函数来拟合高斯结果，该函数的设置与我们在前面章节中熟悉的vglm()函数几乎相同。主要区别是平滑样条的使用，使用另一个函数s()将平滑样条添加到模型中。s()函数采用参数df,，该参数控制变量的平滑样条的最大灵活性。下面的例子用两个预测值拟合 GAM:性别和年龄的平滑样条。summary()函数主要提供平滑项的汇总，p 值测试参数是否与线性趋势有显著差异。请注意，uninormal()系列可以模拟正态分布的位置和比例，尽管默认情况下所有预测值都只针对比例。因此有两个截距，一个是位置截距，GLMs 和 gam 共有的“通常”截距，另一个是分布的标度，它基于方差的自然对数。其原因是VGAM软件包正准备允许预测分布的位置和规模，尽管“经典”的焦点只是预测分布的位置。

mgam <- vgam(CESD11_W1 ∼ Sex + s(AGE_W1, df = 3), data = acl,
        family = uninormal(), model = TRUE)

summary(mgam)

##
## Call:
## vgam(formula = CESD11_W1 ∼ Sex + s(AGE_W1, df = 3), family = uninormal(),
##    data = acl, model = TRUE)
##
##
## Number of linear predictors:     2
##
## Names of linear predictors: mean, loge(sd)
##
## Dispersion Parameter for uninormal family:   1
##
## Log-likelihood: -5290 on 7228 degrees of freedom
##
## Number of iterations:  4
##
## DF for Terms and Approximate Chi-squares for Nonparametric Effects
##
##                   Df Npar Df Npar Chisq P(Chi)
## (Intercept):1      1
## (Intercept):2      1
## Sex                1
## s(AGE_W1, df = 3)  1       2         20      0

使用coef()功能可以查看 GAM 的系数。然而，平滑项的系数不容易解释。但是参数项的系数，在这个模型中只有性别，就像常规的 GLMs 一样可以解释，除了性别的影响控制着年龄的平滑样条函数，而不仅仅是年龄的线性趋势。

coef(mgam)

##     (Intercept):1    (Intercept):2    Sex(2) FEMALE
##            0.2158           0.0435           0.2393
## s(AGE_W1, df = 3)
##           -0.0047

我们可以使用car包中的linearHypothesis()函数得到参数系数的假设检验。

## test parametric coefficient for sex
linearHypothesis(mgam, "Sex(2) FEMALE",
  coef. = coef(mgam), vcov = vcov(mgam))

## Linear hypothesis test
##
## Hypothesis:
## Sex(2) FEMALE = 0
##
## Model 1: restricted model
## Model 2: CESD11_W1 ∼ Sex + s(AGE_W1, df = 3)
##
## Note: Coefficient covariance matrix supplied.
##
##   Res.Df Df Chisq Pr(>Chisq)
## 1   7229
## 2   7228  1  44.1    3.2e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1

可以使用linearHypothesis()来测试参数系数的更复杂的线性假设。例如，我们可以测试截距和性别系数是否同时等于零。

## test parametric coefficient for
## intercept and sex simultaneously
linearHypothesis(mgam,
  c("(Intercept):1", "Sex(2) FEMALE"),
  coef. = coef(mgam), vcov = vcov(mgam))

## Linear hypothesis test
##
## Hypothesis:
## (Intercept):1 = 0
## Sex(2) FEMALE = 0
##
## Model 1: restricted model
## Model 2: CESD11_W1 ∼ Sex + s(AGE_W1, df = 3)
##
## Note: Coefficient covariance matrix supplied.
##
##   Res.Df Df Chisq Pr(>Chisq)
## 1   7230
## 2   7228  2  79.1     <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1

我们可以基于 2 se 以 95%的置信区间来可视化结果。这些应被视为大约 95%的置信区间，因为自由度不是直接估计的，因此依赖于中心极限定理，并且由于平滑样条，标准误差(se)的估计是近似的。使用包含VGAM对象方法的plot()函数可以很容易地绘制图形。结果如图 5-4 所示，显示了性别和年龄平滑项的参数效应。默认情况下会添加地毯图，这对于性别来说不是很有用，但有助于显示年龄数据的分布。颜色来自viridis()调色板。

img/439480_1_En_5_Fig4_HTML.png

图 5-4

将性别作为参数项，年龄作为平滑样条的广义加性模型的模型结果图

par(mfrow = c(1, 2))
plot(mgam, se = TRUE,
     lcol = viridis(4)[1], scol = viridis(4)[2])

将我们 GAM 的结果与更熟悉的替代方案进行比较会有所帮助。为此，我们将安装两个常规 glm。第一个包括年龄的线性项，第二个包括年龄的二次项，使用poly()函数生成二次多项式。

mlin <- vglm(CESD11_W1 ∼ Sex + AGE_W1, data = acl,
        family = uninormal(), model = TRUE)
mquad <- vglm(CESD11_W1 ∼ Sex + poly(AGE_W1, 2), data = acl,
        family = uninormal(), model = TRUE)

为了将这些 glm 与我们的 GAM 进行比较，我们可以制作一个由两个图组成的面板，如图 5-5 所示。深紫色是 GAM 拟合，左边是叠加的线性拟合，右边是施加的二次多项式拟合。显然，线性拟合与 GAM 拟合有很大不同。二次拟合(在图 5-5 的右边板上)相对接近 GAM，除了在最右边的尾部。GAM 在尾部开始变平，而二次趋势继续快速增长。

img/439480_1_En_5_Fig5_HTML.png

图 5-5

双面板图显示了根据广义相加模型预测的抑郁症症状水平，左侧为线性拟合，右侧为二次拟合。

par(mfrow = c(1, 2))
plot(mgam, se = TRUE, which.term = 2,
     lcol = viridis(4)[1], scol = viridis(4)[1])
plot(as(mlin, "vgam"), se = TRUE, which.term = 2,
     lcol = viridis(4)[2], scol = viridis(4)[2],
     overlay = TRUE, add = TRUE)

plot(mgam, se = TRUE, which.term = 2,
     lcol = viridis(4)[1], scol = viridis(4)[1])
plot(as(mquad, "vgam"), se = TRUE, which.term = 2,
     lcol = viridis(4)[3], scol = viridis(4)[3],
     overlay = TRUE, add = TRUE)

有时，gam 可用于选择更简单的趋势函数。例如，我们可以根据图表决定年龄的二次多项式就足够了，并切换到多项式模型。然而，如果结果不能被任何简单的多项式很好地近似，我们可能希望保持 GAM 作为我们的最终模型。在这种情况下，从游戏中得到一些推论是有帮助的。

首先，让我们看看另一个游戏玩家的例子。在这里，我们从性别预测第二波的抑郁症状，第一波的抑郁症状的平滑样条和第一波的年龄的平滑样条。在这个新模型中，我们看到summary()表明年龄的非线性在统计上不显著。

mgam2 <- vgam(CESD11_W2 ∼ Sex +
               s(CESD11_W1, df = 3) +
               s(AGE_W1, df = 3), data = acl,
       family = uninormal(), model = TRUE)

summary(mgam2)

##
## Call:
## vgam(formula = CESD11_W2 ∼ Sex + s(CESD11_W1, df = 3) + s(AGE_W1,
##     df = 3), family = uninormal(), data = acl, model = TRUE)
##
##
## Number of linear predictors:    2
##
## Names of linear predictors: mean, loge(sd)
##
## Dispersion Parameter for uninormal family:  1
##
## Log-likelihood: -3657 on 5725 degrees of freedom
##
## Number of iterations:  5
##
## DF for Terms and Approximate Chi-squares for Nonparametric Effects
##
##                      Df Npar Df Npar Chisq P(Chi)
## (Intercept):1         1
## (Intercept):2         1
## Sex                   1
## s(CESD11_W1, df = 3)  1       2         31    0.0
## s(AGE_W1, df = 3)     1       2          4    0.1

如果其中一个平滑样条与一个线性项没有明显的不同，我们可以考虑回退到该项的线性拟合，如下面的代码所示。

mgam3 <- vgam(CESD11_W2 ∼ Sex +
               s(CESD11_W1, df = 3) +
               AGE_W1, data = acl,
        family = uninormal(), model = TRUE)

summary(mgam3)

##
## Call:
## vgam(formula = CESD11_W2 ∼ Sex + s(CESD11_W1, df = 3) + AGE_W1,
##     family = uninormal(), data = acl, model = TRUE)
##
##
## Number of linear predictors:    2
##
## Names of linear predictors: mean, loge(sd)
##
## Dispersion Parameter for uninormal family:  1
##
## Log-likelihood: -3659 on 5727 degrees of freedom
##
## Number of iterations:  5
##
## DF for Terms and Approximate Chi-squares for Nonparametric Effects
##
##                      Df Npar Df Npar Chisq P(Chi)
## (Intercept):1         1
## (Intercept):2         1
## Sex                   1
## s(CESD11_W1, df = 3)  1       2         31  2e-07
## AGE_W1                1

我们已经看到了如何使用linearHypothesis()函数测试参数项的统计显著性。现在我们可以用它来测试年龄和性别。我们可以使用names()函数和coef()函数来获取每个参数的名称，以便通过测试。

names(coef(mgam3))

## [1] "(Intercept):1"      "(Intercept):2"
## [3] "Sex(2) FEMALE"      "s(CESD11_W1, df = 3)"
## [5] "AGE_W1"

linearHypothesis(mgam3,
  "Sex(2) FEMALE",
  coef. = coef(mgam3), vcov = vcov(mgam3))

## Linear hypothesis test
##
## Hypothesis:
## Sex(2) FEMALE = 0
##
## Model 1: restricted model
## Model 2: CESD11_W2 ∼ Sex + s(CESD11_W1, df = 3) + AGE_W1
##
## Note: Coefficient covariance matrix supplied.
##
##   Res.Df Df Chisq Pr(>Chisq)
## 1   5728
## 2   5727  1  4.09      0.043 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1

linearHypothesis(mgam3,
  "AGE_W1",
  coef. = coef(mgam3), vcov = vcov(mgam3))

## Linear hypothesis test
##
## Hypothesis:
## AGE_W1 = 0
##
## Model 1: restricted model
## Model 2: CESD11_W2 ∼ Sex + s(CESD11_W1, df = 3) + AGE_W1
##
## Note: Coefficient covariance matrix supplied.
##
##   Res.Df Df Chisq Pr(>Chisq)
## 1   5728
## 2   5727  1  3.56      0.059 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1

我们在图 5-6 中展示了最终模型的结果。这里我们可以看到，对于在第一波低于约 2 的抑郁症状水平，较高的抑郁症状预示着在第二波较高的抑郁症状。然而，对于高于约 2 分的抑郁症状，没有太多的关联，这由预测值的变平来证明。

par(mfrow = c(2, 2))
plot(mgam3, se = TRUE,
     lcol = viridis(4)[1],
     scol = viridis(4)[2])

在其他情况下，我们可以为游戏生成预测。因为不可能用语言来描述平滑样条，所以如果对平滑样条的结果感兴趣，那么用图形来表示 gam 是标准的。首先，我们为预测建立一个新的数据集。我们得到了因子变量的所有级别，Sex，第一波的抑郁症状序列，从最低分到最高分，具有 1000 个点的均匀间隔网格，以及第一波的年龄的五个数字摘要。在R中，预测和往常一样使用predict()函数，该函数为VGAM包中的模型准备了方法。

img/439480_1_En_5_Fig6_HTML.png

图 5-6

将性别和年龄作为参数项，将第一波抑郁症状作为平滑样条的广义加性模型的模型结果图

## generate new data for prediction
## use the whole range of sex and depression symptoms
## and a five number summary of age
## (min, 25th 50th 75th percentiles and max)
newdat <- as.data.table(expand.grid(
  Sex = levels(acl$Sex),
  CESD11_W1 = seq(
    from = min(acl$CESD11_W1, na.rm=TRUE),
    to = max(acl$CESD11_W1, na.rm=TRUE),
    length.out = 1000),
  AGE_W1 = fivenum(acl$AGE_W1)))

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

## Warning in `<-.data.table`(x, j = name, value = value): 2 column matrix RHS of := will be treated

## Warning in `[<-.data.table`(x, j = name, value = value): Supplied 20000 items to be assigned to 10000

一旦我们有了一个预测数据集，我们就可以使用ggplot()来制作最终的图表。结果如图 [5-7 所示。该图清楚地强调了后续抑郁症状的最强预测因素是先前的抑郁症状。虽然不同年龄和性别之间有一些微小的差异，但相比之下这些都是小巫见大巫。

img/439480_1_En_5_Fig7_HTML.png

图 5-7

预测不同年龄和性别的第一波抑郁症状水平的第二波抑郁症状

ggplot(newdat,
       aes(CESD11_W1, yhat,
           colour = factor(AGE_W1),
           linetype = factor(AGE_W1))) +
  geom_line() +
  scale_color_viridis("Age", discrete = TRUE) +
  scale_linetype_discrete("Age") +
  facet_wrap(∼ Sex) +
  theme(legend.position = c(.75, .2),
        legend.key.width = unit(1.5, "cm")) +
  xlab("Depression Symptoms (Wave 1)") +
  ylab("Depression Symptoms (Wave 2)")

最后，尽管VGAM包具有广泛的功能，并且是矢量游戏的唯一选择之一，但对于我们在这里展示的单结果游戏，还有一些游戏的特性尚未实现。GAMs 理论和实现的领导者之一 Simon Wood 编写的mgcv包具有一些有用的附加特性。我们将在这里简要概述一下mgcv包的使用，但是更多的细节可以在西蒙·伍德的关于 GAMs 的书中找到。

在我们使用mgcv包之前，我们必须处理一些冲突。具体来说，VGAM和mgcv包都实现了一个名为s()的平滑样条函数，但它们不是同一个函数。加载了两个包后，最后加载的包“屏蔽”了前面包的功能，这意味着当我们在R控制台中键入s()时，我们会从最后加载的包中得到结果，而不一定是我们想要的包。在我们的例子中，解决这个问题最简单的方法是分离不想要的包，并确保加载了想要的包。我们使用下面的代码来实现。请注意，运行后，您将无法使用vgam()功能，直到您重新加载VGAM包。

detach("package:VGAM")
library(mgcv)

现在我们可以使用mgcv包中的gam()函数来拟合一个 GAM。我们再次使用s()函数表示平滑样条，但是注意对于mgcv，控制最大灵活性的参数是k而不是df。适当的家庭功能还有gaussian()。尽管代码的其余部分可能看起来相似，但在默认情况下使用的平滑样条和估计类型方面还是有一些差异。具体来说，mgcv软件包默认使用薄板回归样条和 GCV 准则来学习适当的灵活性程度。

mgam4 <- gam(CESD11_W2 ∼ Sex +
               s(CESD11_W1, k = 3) +
               s(AGE_W1, k = 3), data = acl,
        family = gaussian())

使用gam()功能的一个好处是，默认摘要包括各种有用的信息。具体来说，它自动计算参数项的统计推断，并为平滑项提供近似的显著性测试。请注意，与测试平滑项是否显著不同于线性趋势的vgam()不同，gam()测试平滑项的总体显著性(即，它包括线性趋势以及任何非线性)。

summary(mgam4)

##
## Family: gaussian
## Link function: identity
##
## Formula:
## CESD11_W2 ∼ Sex + s(CESD11_W1, k = 3) + s(AGE_W1, k = 3)
##
## Parametric coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)    -0.0202     0.0272   -0.74    0.457
## Sex(2) FEMALE   0.0681     0.0342    1.99    0.046 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1
##
## Approximate significance of smooth terms:
##               edf Ref.df      F p-value
## s(CESD11_W1) 1.95   2.00 514.03  <2e-16 ***
## s(AGE_W1)    1.63   1.86   2.04   0.085 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1
##
## R-sq.(adj) =  0.271   Deviance explained = 27.2%
## GCV = 0.75609  Scale est. = 0.75461   n = 2867

还因为mgcv提供了对估计自由度的快速访问，更容易检查我们是否应该允许更大的灵活性。例如，如果我们将k = 3增加到k = 4并改装模型，我们可以看到，估计的自由度随年龄变化很小，但随抑郁症状增加。

mgam5 <- gam(CESD11_W2 ∼ Sex +
               s(CESD11_W1, k = 4) +
               s(AGE_W1, k = 4), data = acl,
        family = gaussian())

summary(mgam5)

##
## Family: gaussian
## Link function: identity
##
## Formula:
## CESD11_W2 ∼ Sex + s(CESD11_W1, k = 4) + s(AGE_W1, k = 4)
##
## Parametric coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)    -0.0207     0.0272   -0.76    0.447
## Sex(2) FEMALE   0.0688     0.0342    2.01    0.044 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1
##
## Approximate significance of smooth terms:
##               edf Ref.df      F p-value
## s(CESD11_W1) 2.86   2.99 344.7   <2e-16 ***
## s(AGE_W1)    1.78   2.16   2.4    0.076 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1
##
## R-sq.(adj) =  0.272   Deviance explained = 27.3%
## GCV = 0.75508  Scale est. = 0.75333   n = 2867

我们在图 5-8 中并排绘制了这两个模型。结果显示，在第一波中，抑郁症状的趋势有细微的差异，在我们设定的 GAM 中有更明显的平台期k = 4.

img/439480_1_En_5_Fig8_HTML.png

图 5-8

改变光滑样条最大柔度的两个广义加法模型的模型结果图

par(mfrow = c(2, 2))
plot(mgam4, se = TRUE, scale = 0, main = "k = 3")
plot(mgam5, se = TRUE, scale = 0, main = "k = 4")
互动游戏

另一个在mgcv包中可用但在VGAM包中还不可用的特性是包含交互平滑样条的能力。例如，假设我们认为抑郁症状或年龄的影响可能因性别而异。我们可以相对容易地实现这一点，方法是将参数by = Sex添加到我们认为可能因性别而异的s()函数中。这一总结并没有显示出很大的差异，尽管在性别上有一些差异，但看起来相似的曲线证明了这一点(图 5-9 )。

img/439480_1_En_5_Fig9_HTML.png

图 5-9

允许样条随性别变化的广义可加模型的模型结果图

mgam6 <- gam(CESD11_W2 ∼ Sex +
              s(CESD11_W1, k = 4, by = Sex) +
              s(AGE_W1, k = 4, by = Sex),
             data = acl,
        family = gaussian())

summary(mgam6)

##
## Family: gaussian
## Link function: identity
##
## Formula:
## CESD11_W2 ∼ Sex + s(CESD11_W1, k = 4, by = Sex) + s(AGE_W1, k = 4,
##     by = Sex)
##
## Parametric coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)    -0.0218     0.0276   -0.79    0.428
## Sex(2) FEMALE   0.0693     0.0343    2.02    0.044 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1
##
## Approximate significance of smooth terms:
##                             edf Ref.df      F p-value
## s(CESD11_W1):Sex(1) MALE   2.52   2.83 119.68  <2e-16 ***
## s(CESD11_W1):Sex(2) FEMALE 2.76   2.96 234.18  <2e-16 ***
## s(AGE_W1):Sex(1) MALE      1.80   2.17   1.58     0.2
## s(AGE_W1):Sex(2) FEMALE    1.00   1.00   2.68     0.1
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1
##
## R-sq.(adj) =  0.272   Deviance explained = 27.4%
## GCV = 0.75643  Scale est. = 0.75377   n = 2867

par(mfrow = c(2, 2))
plot(mgam6, ask = FALSE, scale = 0)

我们可以使用 Akaike 信息标准或贝叶斯信息标准作为比较两个模型之间结果的快速方法。在这种情况下，两个指数都指向没有性别交互的模型，作为平衡适合度和节俭度的高级模型。

AIC(mgam5, mgam6)

##         df   AIC
## mgam5  7.6  7333
## mgam6 11.1  7338

BIC(mgam5, mgam6)

##         df   BIC
## mgam5  7.6  7378
## mgam6 11.1  7404

两个连续变量相互作用的光滑样条变得更加复杂。然而，mgcv包允许通过张量积平滑。张量积平滑背后的细节不容易理解，但大致可以认为是取每个变量，通过纽结将其分离并拟合多项式(就像常规样条一样)。这些被假定为乘积项，尽管在整个数据范围内这可能是不现实的，但希望它提供一个合理的近似值，因为可能的空间被两个变量上的结所分解。从实际角度来看，想象一个未知的三维表面，其中深度和宽度由两个预测变量定义，高度是结果的级别。然后想象在顶部覆盖一层厚布。厚重的材料会提供一定程度的“光滑度”，但形状会随着你向任何方向移动而灵活变化，如果不指明另一个变量的水平，你就不能谈论一个变量的“效果”。另一个实际注意事项是张量积平滑对计算要求更高，因此拟合速度较慢。

在下面的代码中，我们拟合了一个 GAM，其中性别作为参数项，张量积在第一波的抑郁症状和自尊之间平滑，预测第二波的抑郁症状。

mgam7 <- gam(CESD11_W2 ∼ Sex +
               te(CESD11_W1, SelfEsteem_W1, k = 4ˆ2),
             data = acl,
        family = gaussian())

summary(mgam7)

##
## Family: gaussian
## Link function: identity
##
## Formula:
## CESD11_W2 ∼ Sex + te(CESD11_W1, SelfEsteem_W1, k = 4ˆ2)
##
## Parametric coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)    -0.0226     0.0268   -0.84    0.400
## Sex(2) FEMALE   0.0718     0.0337    2.13    0.033 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1
##
## Approximate significance of smooth terms:
##                              edf Ref.df    F p-value
## te(CESD11_W1,SelfEsteem_W1) 12.1   14.5 77.4  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1
##
## R-sq.(adj) =  0.286   Deviance explained = 28.9%
## GCV = 0.74276  Scale est. = 0.7391    n = 2867

总结表明，总体而言，张量积平滑在统计上是显著的，尽管它没有说明哪个变量的贡献最大。这种交互的另一个挑战是它们更加难以可视化。现在我们使用vis.gam()功能，它可以制作 3D 透视图或等高线图。首先，我们将在一组图表中从几个不同的角度绘制 3D 透视图。我们还缩小了默认的利润率。结果如图 5-10 所示。

img/439480_1_En_5_Fig10_HTML.png

图 5-10

3D 透视图显示了第 1 波抑郁症状和自尊之间的张量积平滑结果，预测了第 2 波的抑郁

par(mfrow = c(2, 2), mar = c(.1, .1, .1, .1))
vis.gam(mgam7,
  view = c("CESD11_W1", "SelfEsteem_W1"),
  theta = 210, phi = 40,
  color = "topo",
  plot.type = "persp")
vis.gam(mgam7,
  view = c("CESD11_W1", "SelfEsteem_W1"),
  theta = 150, phi = 40,
  color = "topo",
  plot.type = "persp")
vis.gam(mgam7,
  view = c("CESD11_W1", "SelfEsteem_W1"),
  theta = 60, phi = 40,
  color = "topo",
  plot.type = "persp")
vis.gam(mgam7,
  view = c("CESD11_W1", "SelfEsteem_W1"),
  theta = 10, phi = 40,
  color = "topo",
  plot.type = "persp")

在二维空间中更容易可视化的图是等高线图。等值线图在 x 轴和 y 轴上显示预测值，但使用线条和颜色显示第三维。每条线或等高线代表相同的预测值，曲线展示了如何通过改变两个预测值的组合来实现相同的预测值。图 5-11 中显示了一个等高线图示例。

img/439480_1_En_5_Fig11_HTML.png

图 5-11

显示在第 1 波抑郁症状和自尊之间的张量积平滑结果的等高线图预测在第 2 波的抑郁。

par(mfrow = c(1, 1), mar = c(5.1, 4.1, 4.1, 2.1))
vis.gam(mgam7,
  view = c("CESD11_W1", "SelfEsteem_W1"),
  color = "topo",
  plot.type = "contour")

如果我们想要尝试方差的粗略分解，我们可以使用张量积相互作用，使用ti()函数。结果显示在下面的代码中，它们表明抑郁症状和自尊之间的相互作用并没有在抑郁症状和自尊的平滑项之间提供多少附加值。

mgam8 <- gam(CESD11_W2 ∼ Sex +
               ti(CESD11_W1, k = 4) +
               ti(SelfEsteem_W1, k = 4) +
               ti(CESD11_W1, SelfEsteem_W1, k = 4ˆ2),
             data = acl,
        family = gaussian())

summary(mgam8)

##
## Family: gaussian
## Link function: identity
##
## Formula:
## CESD11_W2 ∼ Sex + ti(CESD11_W1, k = 4) + ti(SelfEsteem_W1, k = 4) +
##     ti(CESD11_W1, SelfEsteem_W1, k = 4ˆ2)
##
## Parametric coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)    -0.0156     0.0281   -0.55    0.579
## Sex(2) FEMALE   0.0681     0.0338    2.02    0.044 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1
##
## Approximate significance of smooth terms:
##                               edf Ref.df      F p-value
## ti(CESD11_W1)                1.72   2.08 260.86 < 2e-16 ***
## ti(SelfEsteem_W1)            1.00   1.00  22.77 1.9e-06 ***
## ti(CESD11_W1,SelfEsteem_W1) 21.77  31.08   1.03    0.43
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1
##
## R-sq.(adj) =  0.284   Deviance explained =   29%
## GCV = 0.74778  Scale est. = 0.74087   n = 2867

我们将检查的mgcv包的最后一个附加特性是运行一些快速检查最大允许平滑度是否受限的能力。虽然平滑度是已知的，但参数 k 控制允许的最大值。通常，如果估计的自由度比*k–1 低得多，增加 k 不太可能有任何好处，因为模型已经确定简单的结构就足够了。然而，这并不总是正确的，特别是如果估计的自由度接近于k–1，这可能表明施加的人为限制导致了过度约束的模型，如果我们增加 k 我们可能会得到一组不同的结果。为了检验这一点，我们将回到我们早期的模型，在没有任何交互作用的情况下，检验抑郁症状和预测随后抑郁的年龄。为了方便起见，我们在这里复制了这个模型。请注意，抑郁症症状的估计自由度接近于k–*1 = 4–1 = 3，这表明可能存在一些问题。

mgam5 <- gam(CESD11_W2 ∼ Sex +
               s(CESD11_W1, k = 4) +
               s(AGE_W1, k = 4), data = acl,
        family = gaussian())

summary(mgam5)

##
## Family: gaussian
## Link function: identity
##
## Formula:
## CESD11_W2 ∼ Sex + s(CESD11_W1, k = 4) + s(AGE_W1, k = 4)
##
## Parametric coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)    -0.0207     0.0272   -0.76    0.447
## Sex(2) FEMALE   0.0688     0.0342    2.01    0.044 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1
##
## Approximate significance of smooth terms:
##               edf Ref.df     F p-value
## s(CESD11_W1) 2.86   2.99 344.7  <2e-16 ***
## s(AGE_W1)    1.78   2.16   2.4   0.076 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1
##
## R-sq.(adj) =  0.272   Deviance explained = 27.3%
## GCV = 0.75508  Scale est. = 0.75333   n = 2867

为了检查我们是否想要增加 k ，我们可以使用gam.check()函数。它需要的只是一个合适的游戏。因为gam.check()依赖于一些模拟，它可以根据随机种子而变化。为了确保再现性，设置随机种子，就像我们接下来使用set.seed()所做的那样。打印结果，并绘制一些图，如图 5-12 所示。结果表明，对于第一波的抑郁症状来说，k可能还不够高。一般来说，由于平滑，不需要猜测 k 完全正确，但它需要足够大，以使函数形式不会受到不适当的限制。

img/439480_1_En_5_Fig12_HTML.png

图 5-12

广义加性模型的诊断图

par(mfrow = c(2, 2))
set.seed(12345)
gam.check(mgam5)

##
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 9 iterations.
## The RMS GCV score gradient at convergence was 6.7e-07 .
## The Hessian was positive definite.
## Model rank = 8 / 8
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
##                k'  edf k-index p-value
## s(CESD11_W1) 3.00 2.86    0.97   0.025 *
## s(AGE_W1)    3.00 1.78    0.99   0.230
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1

根据来自gam.check()的信息，我们可能会重新调整我们的模型，增加抑郁症状的 k 。我们选择一个新的值， k = 20。现在结果显示估计的自由度远低于k–1 = 20–1 = 19。虽然我们不需要这样做，但是我们也增加了年龄的 k 来演示改变参数的影响。

mgam5b <- gam(CESD11_W2 ˜ Sex +
               s(CESD11_W1, k = 20) +
               s(AGE_W1, k = 20), data = acl,
        family = gaussian())

summary(mgam5b)

##
## Family: gaussian
## Link function: identity
##
## Formula:
## CESD11_W2 ˜ Sex + s(CESD11_W1, k = 20) + s(AGE_W1, k = 20)
##
## Parametric coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)    -0.0202     0.0271   -0.75    0.456
## Sex(2) FEMALE   0.0680     0.0341    2.00    0.046 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1
##
## Approximate significance of smooth terms:
##                edf Ref.df     F p-value
## s(CESD11_W1) 11.68  13.99 76.43    <2e-16 ***
## s(AGE_W1)     1.69   2.13  1.87      0.14
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1
##
## R-sq.(adj) =  0.278   Deviance explained = 28.2%
## GCV = 0.7508  Scale est. = 0.74678   n = 2867

我们将结果绘制在图 5-13 中，图中显示了抑郁症状趋势比之前包含的更大的灵活性。还要注意，尽管增加了 k ，年龄的趋势并没有发生有意义的变化。这发生在更简单的拟合已经足够的时候；因此，增加 k 对结果几乎没有影响，因为额外的灵活性已经被抵消了。

par(mfrow = c(1, 2))
plot(mgam5b, se = TRUE, scale = 0)

尽管很简短，但本节涵盖了 gam 对于连续的、正态分布的数据的许多基本用途和特性。接下来的部分利用相同的想法，但是应用于其他类型的结果数据，并假设不同的分布。然而，模型拟合、模型比较和可视化的基本步骤往往是相似的。

img/439480_1_En_5_Fig13_HTML.png

图 5-13

增加抑郁症状的 k 值后，广义加法模型的模型结果图

二元结果

二元结果已经在讨论 GLMs 的章节中介绍过了。二元结果的 gam 依赖于与 GLMs 相同的理论，并使用相同的家族(伯努利/二项式)。二元结果 GAMs 的新颖之处在于使用了平滑样条，其功能与连续、正态分布的结果相同，只是平滑样条是在链接尺度上，在二元结果的情况下通常是 logit(对数概率)。

为了检验 gam 的二元结果，我们将使用吸烟作为结果，将当前吸烟者(1)与以前或从不吸烟者(0)进行比较。我们从年龄开始预测。该总结表明，年龄与吸烟状况之间存在某种非线性关系。

library(VGAM)

##
## Attaching package: 'VGAM'

## The following object is masked from 'package:car':
##
##     logit

## The following objects are masked from 'package:rms':
##
##     calibrate, lrtest

## The following object is masked from 'package:mgcv':
##
##     s

## The following objects are masked from 'package:boot':
##
##     logit, simplex

acl$CurSmoke <- as.integer(acl$Smoke_W1 == "(1) Cur Smok")

mgam.lr1 <- vgam(CurSmoke ˜ s(AGE_W1, df = 3),
             family = binomialff(link = "logit"),
             data = acl, model = TRUE)

summary(mgam.lr1)

##
## Call:
## vgam(formula = CurSmoke ˜ s(AGE_W1, df = 3), family = binomialff(link = "logit"),
##     data = acl, model = TRUE)
##
##
## Number of linear predictors:    1
##
## Name of linear predictor: logit(prob)
##
## (Default) Dispersion Parameter for binomialff family:    1
##
## Residual deviance:  4173 on 3613 degrees of freedom
##
## Log-likelihood: -2087 on 3613 degrees of freedom
##
## Number of iterations:  5
##
## DF for Terms and Approximate Chi-squares for Nonparametric Effects
##
##                   Df Npar Df Npar Chisq P(Chi)
## (Intercept)        1
## s(AGE_W1, df = 3)  1        2          44  2e-10

为了更好地理解结果，我们可以绘制平滑样条曲线。

img/439480_1_En_5_Fig14_HTML.png

图 5-14

年龄和当前吸烟状况的广义相加模型

par(mfrow = c(1, 1))
plot(mgam.lr1, se = TRUE,
     lcol = viridis(4)[1],
     scol = viridis(4)[2])

这表明在较年轻的年龄阶段变化不大，但在 60 岁左右急剧下降。然而，目前，该图处于链接比例。我们可以得到概率尺度上的预测，以便用更直观的度量来绘图。

## generate new data for prediction
## use the whole range of age
newdat <- as.data.table(expand.grid(
  AGE_W1 = seq(
    from = min(acl$AGE_W1, na.rm=TRUE),
    to = max(acl$AGE_W1, na.rm=TRUE),
    length.out = 1000)))

newdat$yhat <- predict(mgam.lr1,
                       newdata = newdat,
                       type = "response")

一旦我们有了一个预测数据集，我们就可以使用ggplot()来制作最终的图表。结果如图 5-15 所示。该图显示了在这种情况下与链路规模上的结果类似的图片。然而，转换后的结果更容易解释。

img/439480_1_En_5_Fig15_HTML.png

图 5-15

各年龄段吸烟的预测概率

ggplot(newdat, aes(AGE_W1, yhat)) +
  geom_line() +
  scale_y_continuous(labels = percent) +
  xlab("Age (years)") +
  ylab("Probability of Smoking") +
  coord_cartesian(xlim = range(acl$AGE_W1),
                  ylim = c(0, .4),
                  expand = FALSE)

一个限制是，没有一种内置的方法来获得来自VGAM包的新数据的 GAM 预测的置信区间。生成置信区间的一种方法是使用自举。有了这样一个简单的模型，这并不需要太多时间。我们还可以利用并行处理来加速它，特别是对于更复杂的模型，或者如果我们要进行大量的 bootstrap 重采样。

nboot <- 500

out <- matrix(NA_real_, ncol = nboot, nrow = nrow(newdat))

start.time <- proc.time()
set.seed(12345)
for (i in 1:500) {
  tmp <- vgam(CurSmoke ˜ s(AGE_W1, df = 3),
             family = binomialff(link = "logit"),
             data = acl[sample(nrow(acl), replace = TRUE)], model = TRUE)
  out[, i] <- predict(tmp,
                      newdata = newdat,
                      type = "response")
}
stop.time <- proc.time()

## time to bootstrap 500 times
stop.time - start.time
##    user  system elapsed
##   19.12    0.03   19.21

现在，我们可以从自举预测中生成一些摘要。首先，有时人们会将 bootstrap 预测的平均值与实际模型进行比较，以查看是否存在系统偏差。下面的代码只是一个计算平均绝对差的快速检查。在这种情况下，它非常小。

mean(abs(newdat$yhat - rowMeans(out)))

## [1] 0.00031

接下来，我们可以计算置信区间，但是要取 bootstrap 样本的百分位数。

newdat$LL <- apply(out, 1, quantile,
  probs = .025, na.rm = TRUE)

newdat$UL <- apply(out, 1, quantile,
  probs = .975, na.rm = TRUE)

最后，我们可以重新制作我们的预测概率图，但现在增加了置信区间。结果如图 5-16 所示。由于抽取的 bootstrap 样本数量相对较少，置信区间略有起伏。基于分位数的置信区间往往会随着样本数量的增加而变得平滑，但即使如此，它们也提供了一个相对快速且非常有用的补充，即通过年龄的平滑样条来指示吸烟预测概率的不确定性程度。

img/439480_1_En_5_Fig16_HTML.png

图 5-16

各年龄段吸烟的预测概率

ggplot(newdat, aes(AGE_W1, yhat)) +
  geom_ribbon(aes(ymin = LL, ymax = UL), fill = "grey80") +
  geom_line(size = 2) +
  scale_y_continuous(labels = percent) +
  xlab("Age (years)") +
  ylab("Probability of Smoking") +
  theme_tufte() +
  coord_cartesian(xlim = range(acl$AGE_W1),
                  ylim = c(0, .5),
                  expand = FALSE)

在关于 GLMs 的一章中，我们引入了一个测量方法，即预测概率的平均边际变化，作为当结果是分类的时，对预测值的更直观的总结。虽然这可以处理由于链接函数引起的非线性，但是在平滑样条的情况下，产生这样的测量通常是没有意义的，因为它不仅包含由于开始概率引起的差异，而且包含预测器的效果是非线性的事实。因此，我们在图 5-16 中产生的具有置信区间的数字通常是二元结果博弈的最终结果。

无序的结果

在关于 GLMs 的章节中，我们看到了超过两个水平的无序分类结果的多项逻辑回归模型。对于具有无序分类结果的 gam，过程是相似的。我们首先做一些数据管理，以生成第二波的崩溃就业变量，然后拟合一个 GAM，从年龄的平滑样条预测这一点。该总结揭示了几种对比的显著非线性的证据。一个有趣的特性是，现在不是对平滑项中是否存在非线性进行单一测试，而是进行*k–*1 次测试，其中 k 是结果中类别的数量。在我们的例子中，我们有五个类别；一个用作参考，因此有四个非线性测试。

acl[, EmployG_W2 := as.character(Employment_W2)]
acl[EmployG_W2 %in% c(
  "(2) 2500+HRS", "(3) 15002499",
  "(4) 500-1499", "(5) 1-499HRS"),
  EmployG_W2 := "(2) EMPLOYED"]
acl[, EmployG_W2 := factor(EmployG_W2)]

mgam.mr1 <- vgam(EmployG_W2 ˜ s(AGE_W1, k = 5),
               family = multinomial(),
               data = acl, model = TRUE)

summary(mgam.mr1)

##
## Call:
## vgam(formula = EmployG_W2 ˜ s(AGE_W1, k = 5), family = multinomial(),
##     data = acl, model = TRUE)
##
##
## Number of linear predictors:    4
##
## Names of linear predictors:
## log(mu[,1]/mu[,5]), log(mu[,2]/mu[,5]), log(mu[,3]/mu[,5]), log(mu[,4]/mu[,5])
##
## Dispersion Parameter for multinomial family:    1
##
## Residual deviance:  5261 on 11450 degrees of freedom
##
## Log-likelihood: -2631 on 11450 degrees of freedom
##
## Number of iterations:  8
##
## DF for Terms and Approximate Chi-squares for Nonparametric Effects
##
##                    Df  Npar  Df  Npar  Chisq  P(Chi)
## (Intercept):1       1
## (Intercept):2       1
## (Intercept):3       1
## (Intercept):4       1
## s(AGE_W1, k = 5):1  1         3           16    0.0
## s(AGE_W1, k = 5):2  1         2           83    0.0
## s(AGE_W1, k = 5):3  1         2           71    0.0
## s(AGE_W1, k = 5):4  1         2            5    0.1

我们可以再次绘制结果，默认情况下是在链接范围内，我们现在再次得到四个图，而不是一个，因为我们对分类结果的每个级别都有一个图。结果如图 5-17 所示。一些比较显示出比其他比较更大程度的非线性，并且趋势明显不同。多项式逻辑回归的 gam 允许平滑样条的形状和灵活性在结果的所有级别上变化。

img/439480_1_En_5_Fig17_HTML.png

图 5-17

年龄和就业状况的广义相加模型是一个五级无序分类结果，导致四种不同的年龄效应

par(mfrow = c(2, 2))
plot(mgam.mr1, se = TRUE,
     lcol = viridis(4)[1],
     scol = viridis(4)[2])

同样，我们可能更喜欢生成预测概率，而不是逻辑图。我们照常生成预测。这里，我们使用cbind()函数将预测概率与我们的数据集相结合，因为不是预测概率的一个向量，而是返回一个矩阵，因为在结果的每个级别中都有一个隶属概率。然后，我们将数据融合成一个长数据集，用于绘图。新的数据集有三个变量，一个是年龄，另一个是结果的水平，第三个是实际预测的概率。

## generate new data for prediction
## use the whole range of age
newdat <- as.data.table(expand.grid(
  AGE_W1 = seq(
    from = min(acl$AGE_W1, na.rm=TRUE),
    to = max(acl$AGE_W1, na.rm=TRUE),
    length.out = 1000)))

newdat <- cbind(newdat, predict(mgam.mr1,
                newdata = newdat,
                type = "response"))

newdatlong <- melt(newdat, id.vars = "AGE_W1")

summary(newdatlong)

##      AGE_W1           variable       value
## Min.    :24  (1)  DISABLED:1000  Min.   :0.00
## 1st Qu. :42  (2)  EMPLOYED:1000  1st Qu.:0.03
## Median  :60  (6)  RETIRED :1000  Median :0.08
## Mean    :60  (7)  UNEMPLOY:1000  Mean   :0.20
## 3rd Qu. :78  (8)  KEEP HS :1000  3rd Qu.:0.28
## Max.    :96                      Max.   :0.84

最后，我们可以使用ggplot()绘制结果图。结果如图 5-18 所示。这些发现强调了一个众所周知但不一定被线性模型很好地捕捉到的东西:人们倾向于在 60 岁后退休。因为当许多人(不是所有人，而是许多人)退休时，有一个相对狭窄的年龄窗口，这些模型显然不是线性的。相反，它们相对平坦，短暂地大幅变化，然后又回到相对平坦的状态。鉴于我们对退休年龄的了解，我们可能会考虑围绕该年龄增加更多的灵活性或有目的的分段模型，因为平滑样条仍将试图平滑事实上可能是相对离散的过程。尽管如此，即使没有这样的额外努力，我们也可以看到 GAM 在相对较好地捕捉这种快速转变方面的价值，尽管我们没有为这种转变可能发生的模型提供指导。

img/439480_1_En_5_Fig18_HTML.png

图 5-18

各年龄段就业状况的预测概率

ggplot(newdatlong, aes(
  AGE_W1, value,
  colour = variable, linetype = variable)) +
  geom_line(size = 2) +
  scale_color_viridis(discrete = TRUE) +
  scale_x_continuous("Age (years)") +
  scale_y_continuous("Probability", label = percent) +
  coord_cartesian(ylim = c(0, 1), expand = FALSE) +
  theme_tufte() +
  theme(legend.position = c(.2, .5),
        legend.key.width = unit(2, "cm"))

统计结果

在讨论 GLMs 的章节中已经介绍了计数结果。计数结果的 gam 依赖于与 GLMs 相同的理论，并使用相同的家族(泊松、负二项式)。gam 用于计数结果的新颖之处在于使用了平滑样条，其功能与用于连续、正态分布结果的功能相同，只是平滑样条是在链接尺度上，在计数结果的情况下通常是自然对数。在 GLM 一章中，我们看到了泊松分布的局限性，过度分散可能经常发生并且不能被充分捕捉。因此，这里我们将直接使用负二项式分布来处理 gam，这种分布允许过度分散。尽管我们只研究了一个单一的预测因子，但我们并不局限于此。在下面的例子中，我们检查了慢性疾病的数量，并根据性别和年龄的平滑样条来预测。该总结显示了年龄非线性的有力证据。

## negative binomial regression model
mgam.nbr1 <- vgam(NChronic12_W2 ˜ Sex + s(AGE_W1, k = 5),
              family = negbinomial(),
              data = acl, model = TRUE)

summary(mgam.nbr1)

##
## Call:
## vgam(formula = NChronic12_W2 ˜ Sex + s(AGE_W1, k = 5), family = negbinomial(),
##     data = acl, model = TRUE)
##
##
## Number of linear predictors:  2
##
## Names of linear predictors: loge(mu), loge(size)
##
## Dispersion Parameter for negbinomial family:  1
##
## Log-likelihood: -3636 on 5727 degrees of freedom
##
## Number of iterations:    8
##
## DF for Terms and Approximate Chi-squares for Nonparametric Effects
##
##                  Df Npar Df Npar Chisq P(Chi)
## (Intercept):1    1
## (Intercept):2    1
## Sex              1
## s(AGE_W1, k = 5) 1        3        112      0

为了更好地理解结果，我们可以绘制平滑样条曲线，如图 5-19 所示。我们可以看到，女性比男性报告更多的慢性疾病，报告的疾病数量在年轻时随着年龄的增长而迅速增加，但在老年时增加速度会放缓(图 5-19 )。

img/439480_1_En_5_Fig19_HTML.png

图 5-19

性别、年龄和慢性病数量的广义相加模型

par(mfrow = c(1, 2))
plot(mgam.nbr1, se = TRUE,
     lcol = viridis(4)[1],
     scol = viridis(4)[2])

与二元和多项逻辑回归一样，目前，该图处于链接级别。我们可以在原始计数尺度上获得预测，以便以更直观的度量标准绘图。

## generate new data for prediction
## use the whole range of age and sex
newdat <- as.data.table(expand.grid(
  Sex = levels(acl$Sex),
  AGE_W1 = seq(
    from = min(acl$AGE_W1, na.rm=TRUE),
    to = max(acl$AGE_W1, na.rm=TRUE),
    length.out = 1000)))

newdat$yhat <- predict(mgam.nbr1,
                       newdata = newdat,
                       type = "response")

一旦我们有了一个预测数据集，我们就可以使用ggplot()来绘制图表。结果如图 5-20 所示。虽然性别的影响在链接量表上是不变的，但在最初的反应量表上，随着预测分数的增加，它在绝对量级上增加；因此，预计老年男女之间的差距会更大。但是，注意，目前，年龄和性别之间没有交互作用。因此，在图 5-20 中，年龄对男女的影响是相同的。

img/439480_1_En_5_Fig20_HTML.png

图 5-20

按性别分列的各年龄段慢性疾病的预测数量

ggplot(newdat, aes(AGE_W1, yhat, colour = Sex)) +
  geom_line(size = 2) +
  scale_color_viridis(discrete = TRUE) +
  xlab("Age (years)") +
  ylab("Number Chronic Conditions") +
  theme_tufte() +
  coord_cartesian(xlim = range(acl$AGE_W1),
                  ylim = c(0, 2.5),
                  expand = FALSE) +
  theme(legend.position = c(.2, .8),
        legend.key.width = unit(1, "cm"))

我们可能想知道女性是否真的在年轻时报告了更多的慢性疾病，她们是否比男性更快达到稳定状态。这暗示了一种互动。由于VGAM中目前不支持平滑样条的交互，我们分离包并切换到mgcv。如下图所示，可使用gam()功能拟合模型。

detach("package:VGAM")
library(mgcv)

mgam.nbr2 <- gam(NChronic12_W2 ˜ Sex + s(AGE_W1, k = 10, by = Sex),
              family = nb(), data = acl)

summary(mgam.nbr2)

##
## Family: Negative Binomial(20719.179)
## Link function: log
##
## Formula:
## NChronic12_W2 ˜ Sex + s(AGE_W1, k = 10, by = Sex)
##
## Parametric coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)    -0.2674     0.0395   -6.77  1.3e-11 ***
## Sex(2) FEMALE   0.2661     0.0477    5.58  2.4e-08 ***
## ---
## Signif. codes:    0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1
##
## Approximate significance of smooth terms:
##                          edf Ref.df Chi.sq p-value
## s(AGE_W1):Sex(1) MALE   4.36   5.36  295  <2e-16 ***
## s(AGE_W1):Sex(2) FEMALE 4.11   5.10  447  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '␣' 1
##
## R-sq.(adj) =  0.248   Deviance explained = 25.6%
## -REML = 3649.7  Scale est. = 1         n = 2867

现在，我们再次生成预测，使用的代码与我们给模型安装vgam()时使用的代码相同。

## generate new data for prediction
## use the whole range of age and sex
newdat <- as.data.table(expand.grid(
  Sex = levels(acl$Sex),
  AGE_W1 = seq(
    from = min(acl$AGE_W1, na.rm=TRUE),
    to = max(acl$AGE_W1, na.rm=TRUE),
    length.out = 1000)))

newdat$yhat <- predict(mgam.nbr2,
                       newdata = newdat,
                       type = "response")

最后，我们再次绘制结果，如图 5-21 所示。然而，研究结果显示，如果有什么不同的话(尽管结果可能并不可靠)，女性的平台期远远没有女性的平台期长，男性平台期更早，甚至在老年女性预计会报告更多的慢性疾病。

img/439480_1_En_5_Fig21_HTML.png

图 5-21

根据一个交互模型，按性别预测不同年龄的慢性疾病数量

ggplot(newdat, aes(AGE_W1, yhat, colour = Sex)) +
  geom_line(size = 2) +
  scale_color_viridis(discrete = TRUE) +
  xlab("Age (years)") +
  ylab("Number Chronic Conditions") +
  theme_tufte() +
  coord_cartesian(xlim = range(acl$AGE_W1),
                  ylim = c(0, 2.7),
                  expand = FALSE) +
  theme(legend.position = c(.2, .8),
        legend.key.width = unit(1, "cm"))

与二元逻辑博弈一样，如果我们愿意，我们可以使用 bootstrapping 来生成预测的置信区间。请注意，这比二元逻辑博弈的例子花费的时间要长一些，可能是由于软件的差异，但也因为这是一个更复杂的模型，有两个预测因子和作为性别函数的年龄平滑样条。由于运行时间较长，在这个示例中，如果运行并行处理进行实际分析，可能会更有优势，因为您可能会使用至少几千个引导样本。

nboot <- 500

out <- matrix(NA_real_, ncol = nboot, nrow = nrow(newdat))

start.time <- proc.time()
set.seed(12345)
for (i in 1:500) {
  tmp <- gam(NChronic12_W2 ˜ Sex + s(AGE_W1, k = 10, by = Sex),
              family = nb(),
             data = acl[sample(nrow(acl), replace = TRUE)])
  out[, i] <- predict(tmp,
                      newdata = newdat,
                      type = "response")
}
stop.time <- proc.time()

## time to bootstrap 500 times
stop.time - start.time

##    user  system elapsed
##  167.18    0.73  168.08

现在，我们可以从自举预测中生成一些摘要。首先，有时人们会将 bootstrap 预测的平均值与实际模型进行比较，以查看是否存在系统偏差。下面的代码只是一个计算平均绝对差的快速检查。在这种情况下，它非常小。

mean(abs(newdat$yhat - rowMeans(out)))

## [1] 0.0094

接下来，我们可以计算置信区间，但是要取 bootstrap 样本的百分位数。

newdat$LL <- apply(out, 1, quantile,
  probs = .025, na.rm = TRUE)

newdat$UL <- apply(out, 1, quantile,
  probs = .975, na.rm = TRUE)

最后，我们可以重新制作预测计数图，但现在增加了置信区间。结果如图 5-22 所示。

img/439480_1_En_5_Fig22_HTML.png

图 5-22

使用 bootstrapped 置信区间按性别对不同年龄的慢性病进行预测计数

ggplot(newdat, aes(AGE_W1, yhat)) +
  geom_ribbon(aes(ymin = LL, ymax = UL, fill = Sex), alpha = .2) +
  geom_line(aes(colour = Sex), size = 2) +
  scale_color_viridis(discrete = TRUE) +
  scale_fill_viridis(discrete = TRUE) +
  xlab("Age (years)") +
  ylab("Number Chronic Conditions") +
  theme_tufte() +
  coord_cartesian(xlim = range(acl$AGE_W1),
                  ylim = c(0, 4),
                  expand = FALSE) +
  theme(legend.position = c(.2, .8),
       legend.key.width = unit(2, "cm"))

总的来说，在不同年龄的趋势形态上，女性和男性之间存在一些差异，但在预测中也存在相当大的不确定性，特别是在最年长的年龄，这表明女性超过 80 岁的明显增加可能是也可能不是一个非常可靠的趋势。同样，尽管男性预计在 80 岁后相对稳定，但他们有一个很宽的置信区间，可能包括一个明显的增长。在这种情况下，根本没有足够的数据来确定。事实上，如果我们看看第一波中 80 岁或以上的人在第二波中完成慢性病报告的人数，我们看到只有 27 名男性和 73 名女性，没有太多数据，特别是男性，可以据此对过去 80 年的趋势做出强有力的推断，因此有很大的置信区间。

xtabs(˜Sex + I(AGE_W1 > 80), data = acl[!is.na(NChronic12_W2)])

##             I(AGE_W1 > 80)
## Sex          FALSE TRUE
##   (1) MALE    1010   27
##   (2) FEMALE  1757   73

5.3 总结

本章简要介绍了回归的多项式、样条和光滑样条，然后介绍了一类灵活的模型，广义可加模型(gam)。它展示了 GAMs 如何被用来扩展参数广义线性模型(GLMs)以捕捉预测器的未知函数形式。特别是，当有连续的预测值时，gam 通常是有用的，并且担心与结果的关联可能不能被线性或多项式趋势充分捕获，或者没有足够的信息来推测多项式趋势的程度。在这些情况下，当有足够大的样本量时，GAMs 会表现出色，允许人们捕捉和模拟这些未知的趋势。本章还展示了如何检查结果，以及可视化和呈现 gam 结果的基础知识，包括如何使用 bootstrapping 获得不确定性估计值，以及如何绘制和可视化交互。表 5-1 简要概述了所使用的一些关键功能。

表 5-1

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

|

功能

|

它的作用

| | --- | --- | | vgam() | 来自VGAM包的矢量广义加性模型，用于拟合半参数模型，该模型既包括广义线性模型等标准参数项，也包括一个或多个项的光滑样条。 | | s() | 函数来指示应该对哪个预测变量应用平滑样条。注意，这个函数在VGAM包和mgcv包中以相同的名称出现，但是在VGAM包中，控制灵活性的参数是df，而在mgcv包中，相同的参数是k。 | | gam() | mgcv包中的广义加性模型；见前面的vgam()，因为它们是相似的。 | | plot() | 当应用于广义加性模型时，通常会使用参数化或平滑样条结果生成每个预测变量的图。 | | linearHypothesis() | 测试来自car包的线性假设的函数。允许我们测试关于来自一个vgam()模型的参数项的假设。 | | predict() | 在R中为广义加性模型准备了方法的通用函数，根据模型从原始或新数据生成预测得分。由于平滑样条难以用语言概括，所以在表示广义可加模型时经常使用预测值图。 | | vis.gam() | 通过mgcv包中的gam()功能拟合模型的透视图(3D)和等高线图。对于可视化交互尤其有用。 | | gam.check() | 用于检查所允许的最大灵活性是否(可能)足够，或者灵活性参数 k 是否应该增加。 |

avatar
文章
阅读
粉丝