第 21 章 经济统计综合案例
前面章节分别介绍了数值计算、优化、随机模拟、重采样、模型计算和贝叶斯计算。本章用一个完整的模拟经济统计案例把这些内容串起来。案例主题是家庭收入与消费支出分析。我们关心三个问题：
- 家庭收入和就业状态如何影响消费支出？
- 消费收入比这一统计指标的不确定性有多大？
- 不同模型在预测消费支出时表现如何？
本章所有数据均为模拟数据，不对应任何真实调查来源。这样做的好处是数据生成过程透明，便于教学中检查模型和算法是否按预期工作。
21.1 数据生成与变量说明
设家庭所在地区分为东部、中部和西部。收入服从右偏分布，就业状态与收入相关，消费支出由收入、就业和地区共同决定，并加入随机扰动。
set.seed(2026)
n <- 800
region <- sample(c("东部", "中部", "西部"), n, replace = TRUE,
prob = c(0.42, 0.32, 0.26))
region_effect <- ifelse(region == "东部", 0.20,
ifelse(region == "中部", 0, -0.15))
income <- rlnorm(n, meanlog = log(8) + region_effect, sdlog = 0.55)
employed <- rbinom(n, 1, prob = plogis(-0.5 + 0.18 * income))
consumption <- 1.4 + 0.58 * income + 0.75 * employed +
ifelse(region == "东部", 0.5, ifelse(region == "西部", -0.3, 0)) +
rnorm(n, sd = 1.6)
consumption <- pmax(consumption, 0.2)
survey <- data.frame(region, income, employed, consumption)
head(survey)
#> region income employed consumption
#> 1 <U+4E2D><U+90E8> 16.689 1 10.552
#> 2 <U+4E2D><U+90E8> 10.842 1 10.244
#> 3 <U+4E1C><U+90E8> 5.465 0 4.156
#> 4 <U+4E1C><U+90E8> 8.801 1 9.603
#> 5 <U+4E2D><U+90E8> 4.852 1 3.904
#> 6 <U+4E1C><U+90E8> 3.251 0 5.749变量含义如下：
region：家庭所在地区；income：家庭收入，单位可以理解为万元；employed：家庭主要劳动力是否就业，1 表示就业；consumption：家庭消费支出，单位可以理解为万元。
21.2 数据检查与预处理
经济统计数据常见问题包括缺失值、极端值和变量尺度差异。为了贴近实际工作流，我们人为制造少量缺失，并进行简单处理。
set.seed(1)
miss_id <- sample(seq_len(n), size = round(0.06 * n))
survey$income_obs <- survey$income
survey$income_obs[miss_id] <- NA
c(missing_rate = mean(is.na(survey$income_obs)),
observed_mean = mean(survey$income_obs, na.rm = TRUE),
complete_mean = mean(survey$income))
#> missing_rate observed_mean complete_mean
#> 0.06 10.14 10.06这里缺失是随机制造的，因此完整样本均值和观测样本均值差异不大。真实调查中，如果高收入家庭更倾向于不报告收入，简单删除缺失样本会产生偏差。此处为了聚焦计算流程，使用完整案例分析：
dat <- survey[complete.cases(survey[, c("income_obs", "consumption", "employed", "region")]), ]
dat$log_income <- log(dat$income_obs)
summary(dat[, c("income_obs", "log_income", "consumption")])
#> income_obs log_income consumption
#> Min. : 1.29 Min. :0.253 Min. : 0.20
#> 1st Qu.: 5.83 1st Qu.:1.764 1st Qu.: 5.11
#> Median : 8.70 Median :2.164 Median : 7.17
#> Mean :10.14 Mean :2.152 Mean : 7.95
#> 3rd Qu.:12.29 3rd Qu.:2.509 3rd Qu.: 9.80
#> Max. :63.98 Max. :4.159 Max. :39.33收入右偏，因此建立模型时同时考虑收入水平和对数收入两种设定。
21.3 描述性统计指标
先计算平均收入、平均消费和消费收入比：
\[ R=\frac{\bar C}{\bar Y}, \]
其中 \(\bar C\) 为平均消费，\(\bar Y\) 为平均收入。
ratio_stat <- function(d) mean(d$consumption) / mean(d$income_obs)
c(mean_income = mean(dat$income_obs),
mean_consumption = mean(dat$consumption),
consumption_income_ratio = ratio_stat(dat))
#> mean_income mean_consumption consumption_income_ratio
#> 10.1367 7.9506 0.7843
aggregate(cbind(income_obs, consumption) ~ region, data = dat, mean)
#> region income_obs consumption
#> 1 <U+4E1C><U+90E8> 11.496 9.191
#> 2 <U+4E2D><U+90E8> 9.737 7.493
#> 3 <U+897F><U+90E8> 8.152 6.275描述性统计不是建模的替代品，但它能帮助我们理解变量尺度、组间差异和后续模型是否合理。
21.4 Bootstrap 评估指标不确定性
消费收入比是一个比率统计量，其标准误不如样本均值那样直观。我们使用非参数 Bootstrap 评估其不确定性。
set.seed(2)
B <- 1000
boot_ratio <- numeric(B)
n_dat <- nrow(dat)
for (b in seq_len(B)) {
id <- sample(seq_len(n_dat), size = n_dat, replace = TRUE)
boot_ratio[b] <- ratio_stat(dat[id, ])
}
c(estimate = ratio_stat(dat),
boot_se = sd(boot_ratio),
q025 = quantile(boot_ratio, 0.025),
q975 = quantile(boot_ratio, 0.975))
#> estimate boot_se q025.2.5% q975.97.5%
#> 0.784334 0.007256 0.770653 0.799043这个结果回答的是：如果当前样本可视为总体的代表性样本，那么消费收入比的抽样不确定性大约有多大。
21.5 回归模型设定与计算
考虑两个候选模型：
\[ M_1:\quad C_i=\beta_0+\beta_1Y_i+\beta_2E_i+\gamma_{r(i)}+\varepsilon_i, \]
\[ M_2:\quad C_i=\beta_0+\beta_1\log Y_i+\beta_2E_i+\gamma_{r(i)}+\varepsilon_i. \]
其中 \(C_i\) 是消费支出，\(Y_i\) 是收入，\(E_i\) 是就业状态，\(\gamma_{r(i)}\) 是地区效应。
fit_level <- lm(consumption ~ income_obs + employed + region, data = dat)
fit_log <- lm(consumption ~ log_income + employed + region, data = dat)
summary(fit_level)$coefficients
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 2.0021 0.147377 13.585 1.007e-37
#> income_obs 0.5819 0.009107 63.896 5.550e-305
#> employed 0.6601 0.130109 5.073 4.934e-07
#> region<U+4E2D><U+90E8> -0.6369 0.127463 -4.997 7.254e-07
#> region<U+897F><U+90E8> -0.9503 0.145482 -6.532 1.200e-10
summary(fit_log)$coefficients
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -4.9176 0.3610 -13.621 6.731e-38
#> log_income 5.9391 0.1557 38.156 1.354e-177
#> employed 0.6495 0.1950 3.331 9.077e-04
#> region<U+4E2D><U+90E8> -0.6107 0.1891 -3.229 1.296e-03
#> region<U+897F><U+90E8> -0.7811 0.2177 -3.588 3.545e-04
c(AIC_level = AIC(fit_level),
AIC_log = AIC(fit_log))
#> AIC_level AIC_log
#> 2762 3353如果使用收入水平，收入系数可以解释为收入增加 1 万元时消费支出的平均变化；如果使用对数收入，系数表示收入按比例变化时消费支出的变化。模型选择不仅取决于 AIC，也取决于解释目标。
21.6 交叉验证预测评估
下面使用 5 折交叉验证比较两个模型的预测均方误差。
cv_lm_case <- function(formula, data, K = 5) {
n <- nrow(data)
fold <- sample(rep(seq_len(K), length.out = n))
mse <- numeric(K)
for (k in seq_len(K)) {
train <- data[fold != k, ]
valid <- data[fold == k, ]
fit <- lm(formula, data = train)
pred <- predict(fit, newdata = valid)
mse[k] <- mean((valid$consumption - pred)^2)
}
mean(mse)
}
set.seed(3)
c(level_model = cv_lm_case(consumption ~ income_obs + employed + region, dat),
log_model = cv_lm_case(consumption ~ log_income + employed + region, dat))
#> level_model log_model
#> 2.304 5.108如果预测是主要目标，应重视交叉验证结果；如果政策解释是主要目标，应同时考虑变量形式是否符合经济含义。
21.7 结果解释与报告
一份简洁的结果报告可以按以下顺序组织：
- 数据说明：本例为模拟家庭调查数据，包含收入、就业、地区和消费支出。
- 描述性发现：报告平均收入、平均消费和地区差异。
- 指标不确定性：用 Bootstrap 给出消费收入比的标准误和区间。
- 模型设定：说明为什么比较收入水平和对数收入模型。
- 计算诊断：报告样本量、缺失处理、模型比较指标和交叉验证误差。
- 经济解释：把系数和预测结果翻译成收入、就业与消费之间的关系。
- 局限性：说明数据为模拟数据，真实研究还需要调查设计、权重、缺失机制和因果识别。