如何为分组稳健回归（lm_robust）添加分组观测数？

先解决tidy函数报错问题 #

你遇到的tidy函数找不到的错误，是因为lm_robust（来自estimatr包）的tidy格式化方法由broom包提供，需要先加载broom才能正常使用。

添加子组观测数的两种可行方法 #

方法一：分组回归时直接从模型提取观测数

在do块里先拟合模型，同时提取模型的观测数并合并到结果中：

library(tidyverse)
library(estimatr)
library(broom)

data(iris)

example <- iris %>% 
  group_by(Species) %>% 
  do({
    # 先拟合当前组的稳健回归模型
    model <- lm_robust(Petal.Length ~ Petal.Width + Sepal.Length + Sepal.Width, data = .)
    # 将tidy结果与模型观测数合并
    tidy(model) %>% mutate(nobs = nobs(model))
  }) %>% 
  filter(term == "(Intercept)") %>% 
  select(Species, term, estimate, conf.low, conf.high, p.value, nobs)

# 查看结果
example

方法二：先计算组内观测数，再与回归结果合并

先单独统计每个子组的观测数，再通过关联字段合并到回归结果中：

library(tidyverse)
library(estimatr)
library(broom)

data(iris)

# 第一步：统计每个物种的观测数
group_nobs <- iris %>% 
  group_by(Species) %>% 
  summarise(nobs = n())

# 第二步：执行分组回归并提取系数结果
reg_results <- iris %>% 
  group_by(Species) %>% 
  do(tidy(lm_robust(Petal.Length ~ Petal.Width + Sepal.Length + Sepal.Width, data = .))) %>% 
  filter(term == "(Intercept)") %>% 
  select(Species, term, estimate, conf.low, conf.high, p.value)

# 第三步：合并观测数到回归结果
example <- reg_results %>% left_join(group_nobs, by = "Species")

# 查看结果
example

关于nobs(example)报错的说明 #

nobs()函数是用来提取回归模型对象的观测数，而你之前的example是整理后的tibble数据框，并非模型对象，所以会提示无可用方法。必须从每个子组的回归模型中提取观测数，或者单独统计组内样本量再合并。

内容的提问来源于stack exchange，提问作者skylight