标准误与聚类稳健标准误的理解，基于R实现和Stata实现

1 标准误

1.1 定义

标准误（Standard Error）是用来衡量统计样本估计量（如均值、回归系数等）与总体参数之间的差异的一种统计量。标准误衡量了样本估计量的变异程度，提供了对总体参数的估计的不确定性的度量。标准误越小，表示样本估计量与总体参数的估计越接近，估计越稳定。

1.2 计算公式

$$SE=\frac{\sigma }{\sqrt{n}}$$

2 聚类稳健标准误

聚类稳健标准误的计算方法通常涉及到对观察数据进行分组，然后在每个组内计算残差平方，并最终将这些残差平方加总起来。具体步骤如下：

• 将数据分组： 将观察数据按照聚类结构分成不同的组。

• 计算每个组内的残差平方和： 在每个组内进行回归分析，得到每个观察值的残差（观察值与回归线的差异），然后将这些残差平方加总得到每个组的残差平方和。

• 计算聚类稳健标准误： 将每个组内的残差平方和相加，然后除以总观察数减去组数得到均值，最后取平方根即得到聚类稳健标准误。

2.1 为何聚类之后能降低估计误差？

使用聚类稳健标准误可以降低估计误差，主要是因为它纠正了数据的聚类结构可能导致的异方差性（heteroscedasticity）问题。异方差性是指误差项的方差不是恒定的，而是随着自变量的变化而变化。在具有聚类结构的数据中，观察值往往在同一个聚类内更加相似，这可能导致同一聚类内的观察值之间的误差方差较小，而不同聚类之间的误差方差较大。

在传统的普通最小二乘（OLS）回归中，如果忽略了这种异方差性，估计的标准误可能会被低估。也就是说，估计结果看起来比实际更加精确，而这种低估会使得统计检验的结果产生误导，导致错误的显著性结论。聚类稳健标准误通过将数据分成聚类组并纠正组内相关性，更准确地估计了总体误差的方差，从而避免了异方差性引起的估计误差。

代码实现

# 安装并加载必要的包
# install.packages("boot")
library(boot)

# 模拟数据
set.seed(123)  # 设置随机种子以确保结果的可重复性
n <- 100  # 样本数量
x <- rnorm(n)  # 自变量
y <- 2 * x + rnorm(n)  # 因变量（带有误差项）

# 最小二乘回归
lm_model <- lm(y ~ x)
print('----------------------不加聚类稳健标准误-----------------------------------------')
summary(lm_model)
print('----------------------不加聚类稳健标准误(见上)-----------------------------------------')
# 自助法计算回归系数的标准误
# 自助法计算回归系数的P值、系数、标准误和统计量
boot_results <- boot(data = data.frame(x = x, y = y), statistic = function(data, indices) {
  sampled_data <- data[indices, ]
  lm_result <- lm(y ~ x, data = sampled_data)
  
  # 提取回归系数
  coefficients <- coef(lm_result)
  
  # 计算标准误
  se <- summary(lm_result)$coefficients[, "Std. Error"]
  
  # 计算t统计量
  t_stat <- coefficients / se
  
  # 计算P值
  p_values <- 2 * (1 - pt(abs(t_stat), df = nrow(sampled_data) - length(coefficients)))
  
  # 返回回归系数、P值、标准误和t统计量
  result <- cbind(coefficients, p_values, se, t_stat)
  return(result)
}, R = 1000)  # 进行1000次自助法抽样

print('----------------------加聚类稳健标准误-----------------------------------------')
# 查看回归系数的P值、系数、标准误和t统计量估计结果
head(boot_results$t)
print('----------------------加聚类稳健标准误（见上）-----------------------------------------')




# 加载所需的库
library(ggplot2)
library(dplyr)
library(cluster)

# 生成模拟数据
set.seed(123)
data <- mtcars %>% sample_n(50, replace = TRUE)
if (nrow(data) > 32) {
  warning("Sample size is larger than the number of rows in the data frame.")
} else {
  print(data)
}

# 使用OLS进行回归分析
model <- lm(mpg ~ wt + disp, data = data)
summary(model)

# 计算聚类标准误
kmeans_result <- kmeans(data[, c("wt", "disp")], centers = 3)
data$cluster <- as.factor(kmeans_result$cluster)
ols_by_cluster <- lapply(unique(data$cluster), function(x) {
  cluster_data <- data[data$cluster == x, ]
  model <- lm(mpg ~ wt + disp, data = cluster_data)
  summary(model)$coefficients["(Intercept)"]
})

# 将聚类标准误转换为数据框
ols_by_cluster <- do.call(rbind, ols_by_cluster)
names(ols_by_cluster) <- paste0("Cluster", unique(data$cluster))
ols_by_cluster <- as.data.frame(ols_by_cluster)

# 可视化结果并进行对比
ggplot(data, aes(x = wt, y = mpg)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  labs(title = "OLS回归与聚类回归对比", x = "wt", y = "mpg") +
  theme_minimal() +
  geom_line(aes(y = ols_by_cluster$Cluster1), group = 1, color = "blue", linetype = "dashed") +
  geom_line(aes(y = ols_by_cluster$Cluster2), group = 1, color = "green", linetype = "dashed") +
  geom_line(aes(y = ols_by_cluster$Cluster3), group = 1, color = "orange", linetype = "dashed")




# 加载所需的库
library(ggplot2)
library(dplyr)
library(cluster)

# 生成随机数据
set.seed(123)
data <- mtcars %>% sample_n(100, replace = TRUE)

# 使用OLS进行回归分析
model <- lm(mpg ~ wt + qsec, data = data)
summary(model)
ols_results <- summary(model)$coefficients

# 计算聚类标准误
kmeans_result <- kmeans(data[, c("wt", "qsec")], centers = 3)
cluster_centers <- kmeans_result$centers
cluster_labels <- kmeans_result$cluster

regression_by_cluster <- lapply(1:3, function(cluster_label) {
  cluster_data <- data[cluster_labels == cluster_label, ]
  model <- lm(mpg ~ wt + qsec, data = cluster_data)
  return(summary(model)$coefficients)
})

cluster_se <- lapply(1:3, function(cluster_label) {
  regression_result <- regression_by_cluster[[cluster_label]]
  residuals <- resid(model)
  cluster_residuals <- residuals[cluster_labels == cluster_label]
  return(sqrt(sum((cluster_residuals - mean(cluster_residuals))^2) / (length(cluster_residuals) - 3)))
})

# 可视化结果并进行对比
data_for_plot <- data.frame(wt = data$wt, qsec = data$qsec, mpg = data$mpg, cluster = cluster_labels)
ggplot(data_for_plot, aes(x = wt, y = qsec, color = cluster)) +
  geom_point() +
  geom_line(aes(y = ols_results[1], group = 1), color = "red") +
  geom_line(aes(y = ols_results[2], group = 1), color = "blue") +
  geom_line(aes(y = regression_by_cluster[[1]][1], group = 1), color = "green", linetype = "dashed") +
  geom_line(aes(y = regression_by_cluster[[1]][2], group = 1), color = "purple", linetype = "dashed") +
  labs(title = "OLS回归与聚类之后回归的对比", x = "wt", y = "qsec") +
  theme_minimal()

stata中

*-双向固定效应模型
xtset id year
xtreg y x i.year, fe robust
xtreg y x i.year, fe vce(cluster id)  //与上一条命令等价

【传送门】