visualization-practice.qmd

# 探索实践 {#sec-visualization-practice}

```{r}
#| echo: false

source("_common.R")
```

## 分析老忠实间歇泉喷发规律 {#sec-faithful}

@fig-faithful-bkde2d 展示美国怀俄明州黄石国家公园[老忠实间歇泉](https://en.wikipedia.org/wiki/Old_Faithful)喷发规律，横轴表示喷发持续时间（以分钟计），纵轴表示等待时间（以分钟计），点的亮暗程度（白到黑）代表附近点密度的高低，亮度值通过二维核密度估计方法得到，具体实现借助了 **KernSmooth** [@KernSmooth1995] 包提供的 `bkde2D()` 函数，设置了喷发时间的窗宽为 0.7 分钟，等待时间的窗宽为 7分钟。不难看出，每等待 55 分钟左右间歇泉喷发约 2 分钟，或者每等待 80 分钟左右间歇泉喷发 4.5 约分钟，非常守时，表现得很老实，故而得名。说实话，二维核密度估计在这里有点大材小用了，因为数据点比较少，肉眼也能分辨出来哪里聚集的点多，哪里聚集的点少。

```{r}
#| label: fig-faithful-bkde2d
#| echo: !expr knitr::is_html_output()
#| code-fold: true
#| par: true
#| fig-cap: "二维核密度估计"
#| fig-subcap: 
#| - "faithful 数据集的散点图"
#| - "点的亮暗表示核密度估计值的大小"
#| - "等高线表示核密度估计值"
#| - "等高线表示核密度估计值"
#| fig-width: 4.5
#| fig-height: 4.5
#| fig-showtext: true
#| message: false
#| layout-ncol: 2

# faithful 添加二维核密度估计 density 列
library(KernSmooth)
den <- bkde2D(x = faithful, bandwidth = c(0.7, 7), gridsize = c(51L, 51L))
faithful2d <- expand.grid(eruptions = den$x1, waiting = den$x2) |>
  transform(density = as.vector(den$fhat))

plot(faithful,
  pch = 20, panel.first = grid(), cex = 1, ann = FALSE,
  xlim = c(0.5, 6.5),
  ylim = c(35, 100)
)
title(xlab = "喷发时间", ylab = "等待时间", family = "Noto Serif CJK SC")

plot(faithful,
  pch = 20, panel.first = grid(), cex = 1, ann = FALSE,
  xlim = c(0.5, 6.5),
  ylim = c(35, 100),
  col = densCols(faithful,
    bandwidth = c(0.7, 7),
    nbin = c(51L, 51L), colramp = hcl.colors
  )
)
title(xlab = "喷发时间", ylab = "等待时间", family = "Noto Serif CJK SC")

plot(faithful,
  pch = 20, panel.first = grid(), cex = 1, ann = FALSE,
  xlim = c(0.5, 6.5),
  ylim = c(35, 100),
  col = densCols(faithful,
    bandwidth = c(0.7, 7),
    nbin = c(51L, 51L), colramp = hcl.colors
  )
)
contour(den$x1, den$x2, den$fhat, nlevels = 10, add = TRUE, family = "sans")
title(xlab = "喷发时间", ylab = "等待时间", family = "Noto Serif CJK SC")

# 散点添加颜色
mkBreaks <- function(u) u - diff(range(u)) / (length(u) - 1) / 2
# faithful 划入网格内
xbin <- cut(faithful[, 1], mkBreaks(den$x1), labels = FALSE)
ybin <- cut(faithful[, 2], mkBreaks(den$x2), labels = FALSE)
# 网格对应的核密度估计值即为 faithful 对应的核密度估计值
faithful$dens <- den$fhat[cbind(xbin, ybin)]
# 若是 faithful 数据点没有划分，则置为 0 
faithful$dens[is.na(faithful$dens)] <- 0

library(ggplot2)
library(ggnewscale)
ggplot() +
  geom_point(
    data = faithful, aes(x = eruptions, y = waiting, color = dens),
    shape = 20, size = 2, show.legend = FALSE
  ) +
  scale_colour_viridis_c(option = "D") +
  new_scale_color() +
  geom_contour(data = faithful2d, aes(
    x = eruptions, y = waiting,
    z = density, colour = after_stat(level)
  ), bins = 14, linewidth = 0.45, show.legend = FALSE) +
  scale_colour_viridis_c(option = "C", direction = -1, begin = 0.2, end = 0.8) +
  # colorspace::scale_color_continuous_sequential(palette = "Grays") +
  scale_x_continuous(breaks = 1:6) +
  scale_y_continuous(breaks = 10 * 4:10) +
  coord_cartesian(xlim = c(0.5, 6.5), ylim = c(35, 100)) +
  labs(x = "喷发时间", y = "等待时间", colour = "密度") +
  theme_bw(base_size = 13) +
  theme(
    legend.title = element_text(family = "Noto Serif CJK SC"),
    axis.title = element_text(family = "Noto Serif CJK SC"),
    axis.title.x = element_text(
      margin = margin(b = 0, l = 0, t = 20, r = 0)
    ),
    axis.title.y = element_text(
      margin = margin(b = 0, l = 0, t = 0, r = 20)
    ),
    panel.border = element_rect(color = "black"),
    panel.grid = element_blank(),
    panel.grid.major = element_line(
      color = "lightgray",
      linetype = 3, linewidth = 0.5
    ),
    axis.ticks.length = unit(0.25, "cm"),
    axis.text.x = element_text(
      family = "sans", color = "black",
      vjust = -1.5, size = rel(1.25)
    ),
    axis.text.y = element_text(
      family = "sans", color = "black",
      angle = 90, vjust = 1.5, hjust = 0.5,
      size = rel(1.25)
    )
  )
```

::: callout-tip
函数 `bkde2D()` 实现二维带窗宽的核密度估计（2D Binned Kernel Density Estimate），R 语言存在多个版本，**grDevices** 包的函数 `densCols()` 直接调用 **KernSmooth** 包的函数 `bkde2D()`，**graphics** 包的函数 `smoothScatter()` 与函数 `densCols()` 一样，内部也是调用 `bkde2D()` 函数，**ggplot2** 包的图层 `geom_density_2d()` 采用 **MASS** 包的函数 `kde2d()`，在算法实现上，`MASS::kde2d()` 与 `KernSmooth::bkde2D()` 不同，前者是二维核密度估计（Two-Dimensional Kernel Density Estimation）。Tarn Duong 的著作 [《Multivariate Kernel Smoothing and Its Applications 》](https://www.mvstat.net/mvksa/) [@Chacon2018] 对多元核平滑方法及其应用作了专门的论述，相关实现见书籍配套的 **ks** 包。
:::

```{r}
#| eval: false
#| echo: false

# 三维交互式图形
dens <- with(faithful, MASS::kde2d(x = waiting, y = eruptions))
library(plotly)
with(dens, plot_ly(x = x, y = y, z = z, type = "surface"))
```

## 中国地区级男女性别比分布 {#sec-china-household-sex}

@fig-china-household-sex-1 展示 2020 年中国各省、自治区和直辖市分城市、镇和乡村的性别比数据。数据来自中国国家统计局发布的 2021 年统计年鉴，在数据量不太大的情况下，尽可能展示原始数据，可以捕捉到更加细致的差异。社会经济相关的数据常常显示有马太效应，对原始数据适当做一些变换有利于比较和展示数据，@fig-china-household-sex-2 展示对数变换后的性别比数据的分布。大部分地区的性别比数据都在 100:100 左右，流动人口的性别比波动大很多。

```{r}
#| label: fig-china-household-sex
#| echo: !expr knitr::is_html_output()
#| code-fold: true
#| fig-cap: "2020 年中国地区级男女性别比分布"
#| fig-subcap: 
#|  - 原始性别比数据
#|  - 对数变换后的性别比数据
#| fig-width: 7
#| fig-height: 4.5
#| fig-showtext: true

china_household_sex <- readRDS(file = "data/china-household-sex-2020.rds")
ggplot(data = china_household_sex, aes(x = `户口登记状况`, y = `男性` / `女性`)) +
  geom_jitter(aes(color = `区域`), width = 0.3) +
  theme_classic()

ggplot(data = china_household_sex, aes(x = `户口登记状况`, y = `男性` / `女性`)) +
  geom_jitter(aes(color = `区域`), width = 0.3) +
  scale_y_continuous(trans = "log10") +
  theme_classic()
```

-   「住本乡、镇、街道，户口在本乡、镇、街道」土著和已获得当地户口的。性别比分布有明显的层次差异，性别比均值从大到小依次是城市、乡村、镇。城市里，男性略多于女性，镇里，男性明显少于女性，乡村里，男性略低于女性。
-   「住本乡、镇、街道，户口待定」黑户或其它。性别比分布有明显的层次差异。同上。
-   「住本乡、镇、街道，户口在外乡、镇、街道，离开户口登记地半年以上」流出人口，流出乡、镇、街道。城市、镇、乡村的性别比数据充分混合了，性别比分布没有明显的层次差异。
-   「居住在港澳台或国外，户口在本乡、镇、街道」流出人口，流出国。同上。

## 美国历年各年龄死亡率变化 {#sec-usa-mortality}

<!-- 函数型数据探索、分析和可视化，趋势分析，时间 -->

@fig-usa-mortality 展示美国 1933-2020 年男性分年龄的死亡率数据[^visualization-practice-1]。图分上下两部分，上半部分展示死亡率原值随年龄的变化情况，以 ggplot2 默认的调色板给各个年份配色，下半部分展示死亡率对数变换后随年龄的变化情况，并以红、橙、黄、绿、青、蓝、紫构造彩虹式的调色板给各个年份配色。作图过程中，使用对数变换和调用彩虹式的调色板，帮助我们观察到更多的细节、层次。对数变换后，更加清晰地展示死亡率的变化，尤其是 0-20 岁之间的死亡率起伏变化。调用彩虹式的调色板后，约 20 年为一个阶段，每个阶段内呈现梯度变化，多个阶段体现层次性，更加清晰地展示死亡率曲线的变动趋势。透过层次看到 80 多年来，美国在医疗和公共卫生方面取得的显著改善。

[^visualization-practice-1]: 数据来自德国马克斯普朗克人口研究所、美国加州大学伯克利分校、法国人口研究所共同建立的人类死亡率数据库 (<https://www.mortality.org/>)。

```{r}
#| label: fig-usa-mortality
#| fig-cap: "1933-2020 年美国男性死亡率曲线"
#| fig-width: 6
#| fig-height: 6
#| echo: !expr knitr::is_html_output()
#| fig-showtext: true
#| code-fold: true

usa_mortality <- readRDS(file = "data/usa-mortality-2020.rds")
library(patchwork)
p1 <- ggplot(data = usa_mortality, aes(x = Age, y = Male, group = Year)) +
  geom_vline(xintercept = "100", colour = "gray", lty = 2) +
  geom_line(aes(color = Year), linewidth = 0.25) +
  scale_x_discrete(
    breaks = as.character(20 * 0:5),
    labels = as.character(20 * 0:5)
  ) +
  theme_classic() 
p2 <- p1 +
  labs(x = "年龄", y = "死亡率", color = "年份")
p3 <- p1 +
  scale_y_log10(labels = scales::label_log()) +
  scale_colour_gradientn(colors = RColorBrewer::brewer.pal(name = "Spectral", n = 11)) +
  labs(x = "年龄", y = "死亡率（对数尺度）", color = "年份")
p2 / p3
```

@fig-usa-mortality 也展示了很多基础信息，出生时有很高的死亡率，出生后死亡率迅速下降，一直到10岁，死亡率才又开始回升，直到 20 岁，死亡率才回到出生时的水平。之后，在青年阶段死亡率缓慢增加，直至老年阶段达到很高的死亡率水平。相比于老年阶段，医疗水平的改善作用主要体现在婴儿、儿童、青少年阶段。

@fig-usa-mortality 还展示了一个潜在的数据质量问题，在 100 岁之后，死亡率波动程度明显在变大，这是因为高龄人数变得很少，即死亡率的分母变得很小，分子的细小波动会被放大，也因为同样的原因，100 岁以上的死亡率主要依赖模型估计，甚至出现死亡率大于 1 的罕见情况。因此，就对比医疗和公共卫生水平的变化而言，从数据的实际情况出发，100 岁以上的情况可以不参与比较。

@fig-usa-mortality-heatmap 以年份为横轴，以年龄为纵轴绘制网格，网格内部根据男性死亡率数据填充颜色制作热力图，死亡率数据是对数尺度，颜色的变化和死亡率的变化关系同前，采用了相同的调色板。更加深入的分析和建模，详见 @Marron2022 的第一章。

```{r}
#| label: fig-usa-mortality-heatmap
#| fig-cap: "1933-2020 年美国男性死亡率热力图"
#| fig-width: 6.5
#| fig-height: 4.5
#| echo: !expr knitr::is_html_output()
#| fig-showtext: true
#| code-fold: true

ggplot(data = usa_mortality, aes(x = Year, y = Age, fill = Male)) +
  scale_fill_gradientn(
    colors = RColorBrewer::brewer.pal(name = "Spectral", n = 11),
    trans = "log10", labels = scales::percent_format()
  ) +
  geom_tile(linewidth = 0.4) +
  scale_y_discrete(
    breaks = as.character(10 * 0:10),
    labels = as.character(10 * 0:10),
    expand = c(0, 0)
  ) +
  scale_x_continuous(
    breaks = 1940 + 10 * 0:8,
    labels = 1940 + 10 * 0:8,
    expand = c(0, 0)
  ) + 
  theme_classic() +
  labs(x = "年份", y = "年龄", fill = "死亡率")
```

## 美国弗吉尼亚州城乡死亡率 {#sec-virginia-deaths}

VADeaths 数据来自 Base R 内置的 datasets 包，记录 1940 年美国弗吉尼亚州死亡率，如下表。

```{r}
#| label: tbl-virginia-deaths
#| tbl-cap: 1940 年美国弗吉尼亚州死亡率
#| echo: false

knitr::kable(VADeaths, row.names = TRUE, 
             col.names = c("农村男", "农村女", "城市男", "城市女"))
```

死亡率数据是按年龄段、城乡、性别分组统计的，这是一个三因素交叉统计表，表格中第1行第1列的数据表示 50-54 岁乡村男性的死亡率为 11.7 ‰ ，即在 50-54 岁乡村男性群体中，1000 个人中死亡 11.7 个，这是抽样调查出来的数字。下图分年龄段、城乡、性别展示弗吉尼亚州死亡率数据，从图例来看，突出的是各年龄段的对比，图主要传递的信息是各年龄段的死亡率差异。无论城市还是乡村，也无论男性还是女性，年龄越大死亡率越高，这完全是一个符合生物规律的客观事实，这是众人皆知的，算不上洞见。

```{r}
#| label: fig-virginia-age-group
#| fig-cap: 弗吉尼亚州各年龄段死亡率的对比
#| fig-width: 6
#| fig-height: 5
#| fig-showtext: true
#| echo: !expr knitr::is_html_output()
#| code-fold: true

dat <- transform(expand.grid(
  site = c("乡村", "城镇"), sex = c("男", "女"), 
  age = ordered(c("50-54", "55-59", "60-64", "65-69", "70-74"))
), deaths = as.vector(t(VADeaths)) / 1000)

library(ggplot2)
# (\u2030) 表示千分号
ggplot(data = dat, aes(x = sex, y = deaths, fill = age)) +
  geom_col(position = "dodge2") +
  scale_y_continuous(labels = scales::label_percent(scale = 1000, suffix = "\u2030")) +
  scale_fill_brewer(palette = "Spectral") +
  facet_wrap(~site, ncol = 1) +
  theme_bw(base_size = 13) +
  labs(x = "性别", y = "死亡率", fill = "年龄")
```

将对比对象从年龄段转变为城乡，描述城乡差距在死亡率上的体现，是不是一下子更深刻了呢？城市降低各个年龄段的死亡率，暗示着城市的居住条件、医疗水平比乡村好，提高城市化率增加全民的寿命。严格来说，就这个粗糙的数据集不能如此快地下这个结论，但是，它暗示这个信息，同样也符合人们的常识。

```{r}
#| label: fig-virginia-population-group
#| fig-cap: 弗吉尼亚州城乡死亡率的对比
#| fig-width: 6
#| fig-height: 5
#| fig-showtext: true
#| echo: !expr knitr::is_html_output()
#| code-fold: true

ggplot(data = dat, aes(x = age, y = deaths, fill = site)) +
  geom_col(position = "dodge2") +
  scale_y_continuous(labels = scales::label_percent(scale = 1000, suffix = "\u2030")) +
  scale_fill_brewer(palette = "Spectral") +
  facet_wrap(~sex, ncol = 1) +
  theme_bw(base_size = 13) +
  labs(x = "年龄", y = "死亡率", fill = "城乡")
```

## 如何用图表示累积概率分布 {#sec-diamonds-distr}

不失一般性，考虑连续随机变量的条件分布和累积分布。不妨设随机变量 $x$ 的概率分布函数和概率密度函数分别是 $F(x)$ 和 $f(x)$ 。已知概率分布函数和概率密度函数之间有如下关系。

$$
F(x) = \int_{-\infty}^{x} f(t) \mathrm{dt}
$$

diamonds 数据来自 ggplot2 包，记录了约 54000 颗钻石的价格、重量、切工、颜色、纯净度、尺寸等属性信息。@fig-diamonds-distr 展示这批不同切工的钻石随价格的分布，在这个示例中，如何表达累积分布？概率分布的密度曲线是根据直方图估计得来的，根据不同价格区间内钻石的数量占总钻石的比例估计概率。固定窗宽，即在同一价格区间内累积不同切工的钻石数量，得到累积概率，最后获得累积概率密度曲线，更多理论细节见数据可视化陷阱 [@pu2020] 。

```{r}
#| label: fig-diamonds-distr
#| fig-cap: 不同切工的钻石随价格的分布
#| fig-width: 6
#| fig-height: 5
#| echo: !expr knitr::is_html_output()
#| fig-showtext: true
#| code-fold: true

library(ggplot2)
library(patchwork)
p1 <- ggplot(diamonds, aes(x = price, y = after_stat(density), fill = cut)) +
  geom_density(position = "stack", colour = "white") +
  scale_fill_brewer(palette = "Spectral") +
  scale_y_continuous(
    labels = expression(0, 5~"·"~10^-4, 10 ~ "·" ~ 10^-4, 15 ~ "·" ~ 10^-4),
    breaks = c(0, 5, 10, 15) * 10^(-4)
  ) +
  theme_bw(base_family = "Noto Serif CJK SC") +
  labs(x = "价格", y = "概率密度", fill = "切工", tag = "坏") +
  theme(
    axis.text.x = element_text(family = "sans", color = "black"),
    axis.text.y = element_text(
      family = "sans", color = "black",
      angle = 90, vjust = 1.5, hjust = 0.5
    ),
    legend.text = element_text(family = "sans"),
    plot.tag = element_text(family = "Noto Serif CJK SC", color = "red"),
    plot.tag.position = "topright"
  )

p2 <- ggplot(diamonds, aes(x = price, y = after_stat(density * n), fill = cut)) +
  geom_density(position = "stack", colour = "white") +
  scale_fill_brewer(palette = "Spectral") +
  theme_bw(base_family = "Noto Serif CJK SC") +
  labs(x = "价格", y = "概率质量", fill = "切工", tag = "好") +
  theme(
    axis.text.x = element_text(family = "sans", color = "black"),
    axis.text.y = element_text(
      family = "sans", color = "black",
      angle = 90, vjust = 1.5, hjust = 0.5
    ),
    legend.text = element_text(family = "sans"),
    plot.tag = element_text(family = "Noto Serif CJK SC", color = "black"),
    plot.tag.position = "topright"
  )

p1 / p2
```

## 解释二项总体参数的置信带 {#sec-confidence-belt}

0 和 1 是世界的原点，蕴含大道真意，从 0-1 分布也叫伯努利分布，独立同 0-1分布之和可以衍生出二项分布。在一定条件下，可以用泊松分布近似二项分布。根据中心极限定理，独立同二项分布的极限和与正态分布可以发生关系。在二项分布、正态分布的基础上，可以衍生出超几何分布、贝塔分布等等。本书多个地方以二项分布为例介绍基本统计概念和模型。

在给定置信水平为 0.95，即 $\alpha = 0.05$，固定样本量 $n = 10$，观测到的成功次数 $x$ 可能为 $0,1,\cdots,10$。对于给定的 $p$，不同 $x$ 值出现的机率 $P(X = x)$ 由 $(p + q)^{10}$ 二项展开式的项给出，这里 $q = 1-p$，即：

$$
P(X = x) = \binom{n}{x}p^x(1-p)^{n-x}
$$ {#eq-binom-probability}

在给定 $p = 0.1$ 的情况下，求二项分布的上 $\alpha/2 = 0.025$ 分位点，尾项之和不应超过 $\alpha/2$，最大的 $x$ 值可有如下方程给出：

$$
\sum_{r = x}^{n}\binom{n}{x}p^x(1-p)^{n-x} = \frac{\alpha}{2}
$$ {#eq-binom-quantile}

在 R 语言中，函数 `qbinom()` 可以计算上述二项分布的上分位点 $x = 3$，即

```{r}
qbinom(0.025, size = 10, prob = 0.1, lower.tail = F)
```

反过来，若已知上分位点为 $x = 3$，则概率为

$$
P(X > 3) = \sum_{x > 3}^{10}\binom{10}{x}0.1^x*(1-0.1)^{10-x}
$$ {#eq-binom-prob}

在 R 语言中，函数 `pbinom()` 可以计算上述二项分布的上分位点对应的概率为 $0.0127952$。

```{r}
pbinom(q = 3, size = 10, prob = 0.1, lower.tail = F)
```

首先简单回顾一下置信区间，在学校和教科书里，有两种说法如下：

1.  $1-\alpha$ 的把握确定区间包含真值。
2.  区间包含真值的概率是 $1-\alpha$。

为什么要采纳第一种说法而不是第二种呢？这其实涉及到置信区间的定义问题，历史上 E. S. Pearson 和 R. A. Fisher 曾有过争论。和大多数以正态分布为例介绍参数的置信估计不同，下面以二项分布为例展开介绍。我们知道二项分布是 N 个伯努利分布的卷积，而伯努利分布又称为 0-1 分布，最形象的例子要数抛硬币了，反复投掷硬币，将正面朝上记为 1，反面朝上记为 0，记录正反面出现的次数，正面朝上的总次数又叫成功次数。

1934 年 C. J. Clopper 和 E. S. Pearson 在给定置信水平 $1- \alpha = 0.95$ 和样本量 $n = 10$ 的情况下，给出二项分布 $B(n, p)$ 参数 $p$ 的区间估计（即所谓的 Clopper-Pearson 精确区间估计）和置信带 [@Clopper1934]，如 @fig-confidence-belt 所示，横坐标为观测到的成功次数，纵坐标为参数 $p$ 的置信限。具体来说，固定样本量为 10，假定观测到的成功次数为 2，在置信水平为 0.95 的情况下，Base R 内置的二项精确检验函数 `binom.test()`，可以获得参数 $p$ 的精确区间估计为 $(p_1, p_2) = (0.025, 0.556)$，即：

```{r}
# 精确二项检验 p = 0.2
binom.test(x = 2, n = 10, p = 0.2)
```

值得注意，这个估计的区间与函数 `binom.test()` 中参数 `p` 的取值无关，也就是说，当 $p = 0.4$，区间估计结果是一样的，如下：

```{r}
# 精确二项检验 p = 0.4
binom.test(x = 2, n = 10, p = 0.4)
```

由此，也可以看出区间估计与假设检验的一些关系。

```{=html}
<!-- 
给定置信水平 $1- \alpha = 0.95$ 和样本量 $n = 10$ 的情况下，二项分布 $\mathrm{Bin}(n,p)$ 参数 $p$ 的置信带，如 @fig-confidence-belt 所示。 
-->
```
```{r}
#| label: fig-confidence-belt
#| fig-cap: "二项分布参数的置信带"
#| echo: !expr knitr::is_html_output()
#| code-fold: true
#| fig-width: 5
#| fig-height: 5
#| fig-showtext: true
#| par: true

library(rootSolve) # uniroot.all
options(digits = 4)
# r 为上分位点
p_fun <- function(p, r = 9) qbinom(0.025, size = 10, prob = p, lower.tail = F) - r # 上分位点
l_fun <- function(p, r = 9) qbinom(0.025, size = 10, prob = p, lower.tail = T) - r # 下分位点

# 计算每个分位点对应的最小的概率 p
p <- sapply(0:10, function(x) min(uniroot.all(p_fun, lower = 0, upper = 1, r = x)))

# 计算每个分位点对应的最大的概率 l
l <- sapply(0:10, function(x) max(uniroot.all(l_fun, lower = 0, upper = 1, r = x)))

plot(
  x = seq(from = 0, to = 10, length.out = 11),
  y = seq(from = 0, to = 1, length.out = 11),
  type = "n", ann = FALSE, family = "sans", panel.first = grid()
)
title(xlab = "成功次数", ylab = "置信限", family = "Noto Serif CJK SC")
lines(x = 0:10, y = p, type = "s") # 朝下的阶梯线
lines(x = 0:10, y = p, type = "l") # 折线
# points(x = 0:10, y = p, pch = 16, cex = .8) # 散点

# abline(a = 0, b = 0.1, col = "gray", lwd = 2, lty = 2) # 添加对称线
text(x = 5, y = 0.5, label = "置信带", cex = 1.5, srt = 45, family = "Noto Serif CJK SC")
# points(x = 5, y = 0.5, col = "black", pch = 16) # 中心对称点
# points(x = 5, y = 0.5, col = "black", pch = 3) # 中心对称点

lines(x = 0:10, y = l, type = "S") # 朝上的阶梯线
lines(x = 0:10, y = l, type = "l") # 折线
# points(x = 0:10, y = l, pch = 16, cex = .8) # 散点

points(x = c(2, 2), y = c(0.03, 0.55), pch = 8, col = "black")
text(x = 2, y = 0.55, labels = expression(p[2]), pos = 1)
text(x = 2, y = 0.03, labels = expression(p[1]), pos = 3)
```

## 解释置信区间及其覆盖概率 {#sec-coverage-probability}

<!-- 统计理论、方法的可视化，理论 -->

统计图形很重要的一个作用是解释统计概念，这就要求不拘泥于抽象的严格数学表达，借助数值模拟，可视化等手段帮助读者发散思维，加深理解复杂的逻辑概念，建立统计直觉，正如顾恺之所言「以形写神，形神兼备」。下面仅以二项分布为例讲讲区间估计及其覆盖概率。众所周知，在置信水平为 $1 - \alpha$ 的情况下，二项分布 $\mathrm{Bin}(n,p)$ 的参数 $p$ （也叫成功概率）的 Wald 区间估计为

$$
(\hat{p} - Z_{1-\alpha/2} \sqrt{\hat{p}*(1-\hat{p})/n}, \hat{p} + Z_{1-\alpha/2} \sqrt{\hat{p}*(1-\hat{p})/n})
$$ {#eq-wald-ci}

其中，$n$ 为样本量，$Z_{1-\alpha/2}$ 为标准正态分布 $\mathcal{N}(0,1)$ 在 $1-\alpha/2$ 处的分位点。 $\alpha$ 一般取 0.05，进而 $Z_{1-\alpha/2} \approx 1.96$。用通俗的话说，有 $1 - \alpha$ 的把握确定参数真值 $p$ 在该估计区间内。可见区间估计的意义是解决点估计可靠性问题，但是可靠性和精度往往不能兼得。统计上，通常的做法是先给定可靠性，去尽可能提升精度，即给定置信水平，使估计区间的长度尽可能短，这就涉及到区间估计的方法问题。

下面通过数值模拟的方式辅助说明 Wald 和 Agresti-Coull 两种区间估计方法，现固定样本量 $n = 10$ 或 $n = 100$，重复抽样 1000 次，将参数 $p$ 以 0.01 的间隔离散化，从 0.01 取值到 0.99。已知给定参数 $p$，每次抽样都可以得到参数 $p$ 的估计值 $\hat{p}$ 及其置信区间，1000 次的重复抽样可以计算出来 1000 个置信区间，每个区间要么覆盖真值，要么没有覆盖真值，覆盖的比例可以近似为覆盖概率。

如 @fig-coverage 所示，从上往下分别代表 Wald、 Agresti-Coull、 Wilson 和 Clopper-Pearson 区间估计，纵坐标是覆盖概率，横坐标是参数 $p$ 的真值，图中黑虚线表示置信水平 $1-\alpha=0.95$，红、蓝点线分别表示样本量 $n=10$ 和 $n=100$ 的模拟情况。不难看出，Wald 区间估计方法在小样本情况下表现很差，覆盖概率很少能达到置信水平的，而 Agresti-Coull 区间估计在 Wald 基础上添加了修正后，情况得到显著改善。更多区间估计方法的详细比较见文献 @Blyth1960;@Lawrence2001;@Geyer2005 。

```{r}
#| label: fig-coverage
#| fig-cap: "二项分布参数的几种区间估计：覆盖概率随成功概率的变化"
#| fig-width: 6
#| fig-height: 8.5
#| fig-showtext: true
#| echo: !expr knitr::is_html_output()
#| code-fold: true

# 计算覆盖概率
# Wald 覆盖
coverage_wald <- function(p = 0.1, n = 10, nsim = 1000) {
  phats <- rbinom(nsim, prob = p, size = n) / n
  ll <- phats - qnorm(1 - 0.05 / 2) * sqrt(phats * (1 - phats) / n)
  ul <- phats + qnorm(1 - 0.05 / 2) * sqrt(phats * (1 - phats) / n)
  mean(ll < p & ul > p)
}
# Agresti-Coull 覆盖
coverage_agresti <- function(p = 0.1, n = 10, nsim = 1000) {
  phats <- (rbinom(nsim, prob = p, size = n) + 2) / (n + 4)
  ll <- phats - qnorm(1 - 0.05 / 2) * sqrt(phats * (1 - phats) / n)
  ul <- phats + qnorm(1 - 0.05 / 2) * sqrt(phats * (1 - phats) / n)
  mean(ll < p & ul > p)
}
# Clopper and Pearson (1934)
# 与 binom.test() 计算结果一致
coverage_clopper <- function(p = 0.1, n = 10, nsim = 1000) {
  nd <- rbinom(nsim, prob = p, size = n)
  ll <- qbeta(p = 0.05 / 2, shape1 = nd, shape2 = n - nd + 1)
  ul <- qbeta(p = 1 - 0.05 / 2, shape1 = nd + 1, shape2 = n - nd)
  mean(ll < p & ul > p)
}
# Wilson (1927)
# 与 prop.test(correct = FALSE) 计算结果一致
coverage_wilson <- function(p = 0.1, n = 10, nsim = 1000) {
  phats <- rbinom(nsim, prob = p, size = n) / n
  lambda <- qnorm(1 - 0.05 / 2)
  ll <- phats + lambda^2 / (2 * n) - lambda * sqrt(phats * (1 - phats) / n + lambda^2 / (4 * n^2))
  ul <- phats + lambda^2 / (2 * n) + lambda * sqrt(phats * (1 - phats) / n + lambda^2 / (4 * n^2))
  mean(ll / (1 + lambda^2 / n) < p & ul / (1 + lambda^2 / n) > p)
}

sim_dat <- transform(expand.grid(
  p = seq(0.01, 0.99, by = 0.01),
  n = c(10, 100),
  nsim = 1000,
  methods = c("Wald", "Agresti-Coull", "Wilson", "Clopper-Pearson")
), prob = ifelse(methods == "Wald",
  Vectorize(coverage_wald)(p = p, n = n, nsim = nsim),
  ifelse(methods == "Agresti-Coull",
    Vectorize(coverage_agresti)(p = p, n = n, nsim = nsim),
    ifelse(methods == "Wilson",
      Vectorize(coverage_wilson)(p = p, n = n, nsim = nsim),
      Vectorize(coverage_clopper)(p = p, n = n, nsim = nsim)
    )
  )
), nsample = ifelse(n == 10, "n=10", "n=100"))

ggplot(data = sim_dat, aes(x = p, y = prob, color = nsample)) +
  geom_hline(yintercept = 0.95, linetype = 2, 
             linewidth = 1, color = "gray60") +
  geom_point() +
  geom_path() +
  # annotate(geom = "text", x = 0, y = 0.95, label = "0.950",
  #          fontface = "bold", hjust = 2, size = 3.5) +
  # scale_color_grey() +
  scale_color_brewer(palette = "Set1") +
  facet_wrap(facets = ~methods, ncol = 1, scales = "free_y") +
  labs(x = "成功概率", y = "覆盖概率", color = "样本量") +
  theme_bw(base_size = 13, base_family = "sans") +
  theme(title = element_text(family = "Noto Serif CJK SC")) + 
  coord_cartesian(clip = 'off')
```

通过 @fig-coverage 一看就明白了几种区间估计方法的优劣，以及为什么软件普遍默认采用 Wilson 估计方法？因为它又稳定又准确。 Wilson 区间估计用的更加广泛的，Base R 内置的比例检验函数 `prop.test()` 在不启用 Yates 修正时，就是用该方法获得比例 $p$ 的区间估计 [@Wilson1927]。Clopper-Pearson 区间估计特别适合小样本情形，它是精确区间估计方法，Base R 内置的二项比例检验函数 `binom.test()` 就是用该方法获得比例 $p$ 的区间估计[@Clopper1934]。

::: callout-tip
请读者再思考两个问题： @fig-coverage 为什么呈现对称的形式，泊松分布会和二项分布有类似的现象吗？如果有的话，连续分布，如正态分布和指数分布也会有吗？
:::

## 习题 {#sec-exercise-practice}

1.  1888 年，瑞士开始进入一个人口转变的阶段，从发展中国家的高出生率开始下滑。分析生育率和经济指标的关系。数据集 swiss 记录了 1888 年瑞士 47 个说法语的省份的生育率和社会经济指标数据。Fertility（生育率，采用常见的标准生育率统计口径）、Agriculture（男性从事农业生产的比例）、Examination（应征者在军队考试中获得最高等级的比例）、Education（应征者有小学以上教育水平的比例）、Catholic（信仰天主教的比例）、Infant.Mortality（婴儿死亡率，仅考虑出生一年内死亡），各个指标都统一标准化为百分比的形式。其中，Examination 和 Education 是 1887 年、1888 年和 1889 年的平均值。瑞士 182 个地区 1888 年及其它年份的数据可从[网站](https://opr.princeton.edu/archive/pefp/switz.aspx)获得。

    ```{r}
    #| eval: false
    #| echo: false

    library(GGally)
    ggpairs(swiss, columns = 1:4, aes(color = Catholic > 50),
            upper = "blank", progress = FALSE)
    ```