Stata+R：你了解什么是抽样吗？Sampling

👇 连享会 · 推文导航 | www.lianxh.cn

• 🍎 Stata： Stata基础 | Stata绘图 | Stata程序 | Stata新命令
• 📘 论文： 数据处理 | 结果输出 | 论文写作 | 数据分享
• 💹 计量： 回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归
• ⛳ 专题： SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析
• 🔃 因果： DID | RDD | 因果推断 | 合成控制法 | PSM-Matching
• 🔨 工具： 工具软件 | Markdown | Python-R-Stata
• 🎧 课程： 最新专题 | 计量专题 | 关于连享会

作者： 杨凡佳 (中央财经大学)
邮箱： [email protected]

Source： Chester Ismay and Albert Y. Kim, 2022. Statistical Inference via Data Science: A ModernDive into R and the Tidyverse. -Link- . Chapter 7 Sampling

温馨提示： 文中链接在微信中无法生效。请点击底部 「阅读原文」 。或直接长按/扫描如下二维码，直达原文：

目录

• 1. 背景介绍

• 2. 手动抽样

• 2.1 单次抽样

• 2.3 重复抽样

• 3. 计算机抽样

• 3.1 单次抽样

• 3.2 重复抽样

• 4. 抽样的理论基础

• 4.1 抽样的基础概念

• 4.2 什么是好的估计值

• 4.3 中心极限定理

• 5. 用 Stata 模拟抽样

• 5.1 生成总体

• 5.2 单次抽样

• 5.3 重复抽样

• 6. 小结

• 7. 相关推文

1. 背景介绍

抽样 (sampling) 是统计学中最为基础的概念之一，在难以确知总体分布的情况下，我们都需要通过抽样对样本数据进行研究。但是，你真的了解什么是抽样吗？下面，我们借鉴一本英文教材 Statistical Inference via Data Science 上的讲解和案例，通过一些最基本的问题，重新建构对抽样的认知，从底层原理认识统计学中利用样本推断总体的方法论。教材中附有 R 语言的代码和操作，我们用 Stata 代码进行了同样过程的复现，这也将有助于读者熟悉和使用 R 和 Stata 对统计抽样的模拟与编程。

我们需要先下载并调用如下 R 程辑包 (packages)：

library (tidyverse) 
library (moderndive) 
library (ggplot2) 
library (dplyr) 
library (tidyr) 
library (readr) 
library (purrr) 
library (tibble) 
library (stringr) 
library (forcats) 

2. 手动抽样

现在我们从最简单的统计学问题出发：如下图所示，有一只 碗 (bowl) ，里面装有若干个红球和若干个白球，红球占全部球的 比例 (proportion) 是多少？

2.1 单次抽样

一个最直接的方法是进行挨个的计数：取出所有的球，数清红球的数量和白球的数量，然后用红球的数量除以球的总数。然而，这将是一个漫长而乏味的过程。下面，我们将尝试一种更简单的办法。 “人类与动物的最大区别就是制造工具和使用工具”。使用下图的 铲子 (shovel) ，我们每次可以取出 50 个红球和白球。第一次的结果中，我们取出了 17 个红球，计算得到红球的比例是 17÷50 = 34%。 现在，我们可以把 34% 看作是对整个碗中红球比例的一个猜测。虽然这不像数清碗里所有球那么精确，但我们得到 34% 这一猜测值花费的时间和精力要少得多。

2.3 重复抽样

然而，上述单次抽样的方法可靠吗？显然，如果我们每一次用铲子手动抽取 50 个球，这一比例不可能永远都是 34%。因此，我们 (在教材中，是 33 组同学) 重复抽取 33 次。每完成一次抽取，将球放回碗里并混合碗里的红球和白球，从而避免上一次的抽样结果干扰下一次的抽样结果。 如果我们用直方图记录这 33 次手动抽样中红球的比例，可以将这一变量的数值分布可视化。如下图所示，我们可以发现这种分布类似于钟形。

使用 R 语言 可以对这一分布得到更为美观的可视化结果。

tactile_prop_red

ggplot (tactile_prop_red， aes (x = prop_red！ +
  geom_histogram (binwidth = 0.05, boundary = 0.4, color = "white") +
  labs (x = "Proportion of 50 balls that were red", 
       title = "Distribution of 33 proportions red")

输出图如上，可以看出直方图具有一定的钟形特征。

3. 计算机抽样

3.1 单次抽样

很明显，手动抽样受限于人的时间和精力，用计算机进行虚拟 (virtual) 抽样则不存在这一问题。我们使用 R 语言对上述手动抽样过程可以进行模拟 首先，我们定义一个 “碗”，格式是一个 数据框 (dataframe)。

bowl

# A tibble: 2,400 × 2
   ball_ID color
      
 1       1 white
 2       2 white
 3       3 white
 4       4 red  
 5       5 white
 6       6 white
 7       7 red  
 8       8 white
 9       9 red  
10      10 white
# ℹ 2,390 more rows

能看到这个 “碗” 里包括了 2400 个球，其中包括红球和白球，每一个球有一个 序号 (ball_ID)。 接下来，在这个碗里，我们使用一个 “铲子” 进行抽样，抽出 50 个球。

virtual_shovel % 
  rep_sample_n(size = 50)
virtual_shovel

# A tibble: 50 × 3
# Groups:   replicate [1]
   replicate ball_ID color
           
 1         1     177 red  
 2         1    1340 white
 3         1    2220 red  
 4         1    2245 red  
 5         1    1999 white
 6         1     732 white
 7         1    1650 red  
 8         1    1938 white
 9         1    1400 white
10         1       8 white
# ℹ 40 more rows

可以看出，这 50 个球已经被抽出，分别是编号为 177 的红球，编号为 1340 的白球…… 最后，我们对抽样结果进行统计，从而得到单次抽样中的红球比例。

virtual_shovel %>% 
  summarize(num_red = sum(color == "red")) %>% 
  mutate(prop_red = num_red / 50)

# A tibble: 1 × 3
  replicate num_red prop_red
             
1         1      20      0.4

通过简单的计算，在 50 个球中红球比例是 40%。这是单次抽样得到的统计结果。

3.2 重复抽样

3.2.1 reps = 33

利用 R 语言中的循环函数 (在 rep_sample_n ! 中已经内置)，我们将上述单次抽样的过程重复 33 次，得到新的直方图，这一直方图同样具有钟形特征。

virtual_samples % 
  rep_sample_n(size = 50, reps = 33)

输出直方图如下：

3.2.2 reps = 1000

如果我们将上述单次抽样的过程重复 1000 次呢，只需要将 reps = 33 修改成 reps = 1000 。

virtual_samples % 
  rep_sample_n(size = 50, reps = 1000)

输出直方图如下

此时，直方图的钟形特征更为明显，分布形态的 “钟形曲线” 也更为光滑。

3.2.3 size = 25，size = 50，size = 100

如果我们用不同大小的铲子，也会有不同的抽样结果。例如，我们使用 size = 25，size = 50，size = 100 分别进行 1000 次抽样，得到三个不同的分布。

# Segment 1: sample size = 25 ------------------------------
# 1.a) Virtually use shovel 1000 times
virtual_samples_25 % 
  rep_sample_n(size = 25, reps = 1000)

# 1.b) Compute resulting 1000 replicates of proportion red
virtual_prop_red_25 % 
  group_by(replicate) %>% 
  summarize(red = sum(color == "red")) %>% 
  mutate(prop_red = red / 25)

# 1.c) Plot distribution via a histogram
ggplot(virtual_prop_red_25, aes(x = prop_red)) +
  geom_histogram(binwidth = 0.05, boundary = 0.4, color = "white") +
  labs(x = "Proportion of 25 balls that were red", title = "25") 


# Segment 2: sample size = 50 ------------------------------
# 2.a) Virtually use shovel 1000 times
virtual_samples_50 % 
  rep_sample_n(size = 50, reps = 1000)

# 2.b) Compute resulting 1000 replicates of proportion red
virtual_prop_red_50 % 
  group_by(replicate) %>% 
  summarize(red = sum(color == "red")) %>% 
  mutate(prop_red = red / 50)

# 2.c) Plot distribution via a histogram
ggplot(virtual_prop_red_50, aes(x = prop_red)) +
  geom_histogram(binwidth = 0.05, boundary = 0.4, color = "white") +
  labs(x = "Proportion of 50 balls that were red", title = "50")  


# Segment 3: sample size = 100 ------------------------------
# 3.a) Virtually using shovel with 100 slots 1000 times
virtual_samples_100 % 
  rep_sample_n(size = 100, reps = 1000)

# 3.b) Compute resulting 1000 replicates of proportion red
virtual_prop_red_100 % 
  group_by(replicate) %>% 
  summarize(red = sum(color == "red")) %>% 
  mutate(prop_red = red / 100)

# 3.c) Plot distribution via a histogram
ggplot(virtual_prop_red_100, aes(x = prop_red)) +
  geom_histogram(binwidth = 0.05, boundary = 0.4, color = "white") +
  labs(x = "Proportion of 100 balls that were red", title = "100") 

输出的三个分布的直方图如下：

基于三个分布，我们还可以计算出分布的标准差。

virtual_prop_red_25 %>% 
  summarize(sd = sd(prop_red))

# A tibble: 1 × 1
      sd
   
1 0.0985

# n = 50
virtual_prop_red_50 %>% 
  summarize(sd = sd(prop_red))

# A tibble: 1 × 1
      sd
   
1 0.0669

# n = 100
virtual_prop_red_100 %>% 
  summarize(sd = sd(prop_red))

# A tibble: 1 × 1
      sd
   
1 0.0487

能看出，随着单次抽样个体数量 (size) 的增加，标准差逐渐缩小。

4. 抽样的理论基础

回顾上述的手动抽样和计算机抽样过程，我们提出以下关于抽样的知识体系。

4.1 抽样的基础概念

第一组概念关于总体：总体 (population)、总体参数 (population Parameters) 和普查 (census)。

• 总体 是指研究对象的全部个体或观测值的集合，是我们想要研究和了解的完整群体。
• 总体参数 是用来描述总体特征的数值度量。常见的总体参数包括总体均值 (平均值)、总体方差、总体标准差等。这些参数是固定的，但通常未知，需要通过样本统计量来估计。
• 普查 是指对整个总体进行的调查或观测。在普查中，我们收集并分析总体中每一个个体的数据。由于普查可以提供最全面的信息，因此它的结果通常被认为是非常准确的，但成本和时间消耗也很高。

在我们之前的手动抽样过程中，总体是碗中的所有球，总体参数是我们关心的红球比例，普查是穷尽计数碗中的所有球并获取红球比例。

在我们之前的计算机抽样过程中，总体是数据框 bowl 中带有编号 (ball_ID) 和颜色 (color) 的 2400 个数据，总体参数是我们每次抽样计算出的红球比例，普查是对数据框 bowl 中 2400 个数值直接计算红球比例。

第二组概念关于样本：样本 (sample) 和样本统计量 (sample Statistics) / 点估计量 (point estimate)。

• 样本 是从总体中选取的一部分个体或观测值。样本的目的是代表总体，通过对样本的研究来推断总体的特征。
• 样本统计量 是根据样本数据计算得到的数值度量，用来估计总体参数。例如，样本均值、样本方差和样本标准差等。样本统计量是随机变量，因为它们是基于随机抽取的样本计算的。

在我们之前的手动抽样过程中，样本是每次利用铲子取出的 50 个球，样本统计量是 50 个球中的红球比例。

在我们之前的计算机抽样过程中，样本是每次使用函数抽出的数据 (25 个、50 个或 100 个)，样本统计量是基于抽出的样本计算出的红球比例。

第三组概念关于抽样的方法论：代表性 (representative)、推广性 (generalizable) 和无偏性 (unbiased)。

• 代表性 是指样本能够准确反映总体特征的程度。一个具有代表性的样本应该在结构和分布上与总体相似，这样从样本得到的结论才能有效推广到总体。
• 推广性 是指从样本中得到的结论可以被推广到整个总体或其他相似总体的程度。在统计学中，一个具有良好推广性的抽样设计意味着样本能够代表总体，从而使得基于样本数据的统计推断在总体上也是有效的。
• 无偏性 是指样本统计量的平均值等于总体参数的真实值。一个无偏的估计量可以确保随着样本容量的增加，样本统计量会趋近于总体参数。

在我们之前的手动抽样过程中，代表性是用铲子取出的球能够 “大致类似于” 碗里的其他球，推广性是用铲子取出的红球比例可以推广到碗里所有球的红球比例。如果能够保证上述过程，那么抽样是无偏的。然而，如果我们在抽样过程中人为 “作弊” (cheat)，例如在使用铲子取球时 (希望取出更高比例的红球) 故意将一些白球替换成红球，那么这种抽样结果不具有代表性和推广性，也是有偏的。我们抽样得到的估计值也不是一个好的猜测值 (good guess)。

在我们之前的计算机抽样过程中，只要我们的抽样代码是正确无误的，这种代表性、推广性和无偏性天然地得以满足。

第四组概念关于抽样的目标：随机抽样 (Random Sampling) 和统计推断 (Statistical Inference)。

• 随机抽样 是一种抽样方法，其中每个样本成员被选中的概率是已知的，通常是相等的。这种方法可以确保样本的代表性，是进行有效统计推断的基础。
• 统计推断 是指使用样本数据来估计总体参数或测试关于总体参数的假设的过程，包括参数估计 (如点估计和区间估计) 和假设检验。

在之前的手动抽样过程中，如果抽样过程中没有作弊，并且不断地将碗中的球搅拌均匀，那么就是随机抽样。基于抽样结果估计红球比例的过程就是统计推断。

假如机器并不能作弊，也不存在系统性偏差，那它的抽样过程可以视为随机抽样。基于每一次的抽样结果计算红球比例，就是一个统计推断 (点估计) 的过程。

最后，还要用到两个概念， 抽样分布 和 标准误 。

• 点估计的 抽样分布 是指从总体中抽取样本，并根据样本数据计算出的统计量 (如样本均值、样本方差等) 的分布。这些统计量是基于样本数据的，因此它们本身也具有随机性，即每次抽取的样本可能不同，导致计算出的统计量也会有所不同。点估计的抽样分布描述了这些统计量在多次抽样中的变化情况。
• 点估计的 标准误 (Standard Error，简称 SE) 是衡量点估计量精确度的一个指标。它表示在多次抽样中，点估计量与其真实总体参数 (如总体均值) 之间的平均偏差的平方根。标准误可以用来评估样本统计量的代表性，即它越小，表示样本统计量越接近总体的真实参数。

4.2 什么是好的估计值

我们每次的抽样不同，抽样统计量的计算结果也就不同。什么是一个好的估计值？让我们思考两个问题。

第一个问题：无偏性 (准确性，accuracy)。尽管每一次抽样的计算结果可能高于真实值，也可能低于真实值。但平均而言，估计值以真实值为中心变化，是 “平均” 正确的。事实上，这就是我们上面所介绍的无偏性。

在我们之前的抽样过程中，总体 2400 个球中，有 900 个球是红球，真实的红球比例是 37.5%。虽然每一次估计的结果并不完全等于 37.5%，但是始终围绕 37.5% 上下波动。

第二个问题：精确性 (precision)。对比三个分布可以看出，随着样本量 n 的增加，红球比例点估计值越来越集中在 37.5% 周围。这种集中的趋势可以用样本的标准误减小来衡量，这反映了抽样估计的结果越来越精确。

用一个四象限图来理解抽样的无偏性和精确性。严格的随机抽样能够保证点估计量具有无偏性，大样本量能够使点估计量更具有精确性。就像是向靶上射击，保证每次射击不能偏离靶中心是无偏性。而保证射击的结果尽可能每次集中是精确性。我们希望的结果是抽样既具有无偏性，又具有精确性。

4.3 中心极限定理

在实际生活中，我们很难对总体进行普查，只能通过抽样计算样本统计量。抽样统计的方式是否是可靠的呢？这决定了我们当前统计学的底层方法论是否可靠。对此解答的理论就是著名的中心极限定理，中心极限定理的表达包括以下两条。

(1) 随着样本量的增加，点估计的抽样分布越来越遵循均值为总体真实值正态分布。 (2) 随着样本量的增加，点估计抽样分布的标准误变得越来越小。

中心极限定理成为我们利用样本研究总体的桥梁。因为我们每一次抽样得到的样本统计量都是正态分布中的一个结果，而这个正态分布始终围绕总体的真实值，并且随着样本量的扩大越来越集中于总体的真实值。正因如此，我们的抽样才是有意义的！

5. 用 Stata 模拟抽样

本节用 Stata 来简要演示常用抽样过程。对此部分内容感兴趣的读者可以预先阅读如下推文：

• 陈勇吏, 2020, Stata程序：Monte-Carlo-模拟之产生符合特定分布的随机数 , 连享会 No.343.
• 陈勇吏, 2020, Stata：蒙特卡洛模拟分析 (Monte Carlo Simulation) , 连享会 No.47.
• 滕泽鹏, 2023, Stata：随机抽样命令介绍-gsample , 连享会 No.1140.
• 侯新烁, 2020, Stata：蒙特卡洛模拟A-(Monte-Carlo-Simulation)没那么神秘 , 连享会 No.405.
• 李坤, 2022, Stata：蒙特卡洛模拟新命令-simulate2 , 连享会 No.931.

5.1 生成总体

clear all
set obs 2400
gen ball_ID = _n
gen color = ""

forvalues i = 1/2400 {
    local r = runiform (0， 1) 
    if `r' < 0.5 {
        replace color = "Red" in `i'
    }
    else {
        replace color = "White" in `i'
    }
	}
	
tab color
save "bowl.dta", replace

输出结果如下：