标准误与聚类稳健标准误的理解

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()



1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136

stata中

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

【传送门】