ggplot2：如何为geom_smooth中的预测获得稳健的置信区间？

2023年2月16日 493次阅读

考虑这个简单的例子

dataframe <- data_frame(x = c(1,2,3,4,5,6),
                        y = c(12,24,24,34,12,15))
> dataframe
# A tibble: 6 x 2
      x     y
  <dbl> <dbl>
1     1    12
2     2    24
3     3    24
4     4    34
5     5    12
6     6    15    

dataframe %>% ggplot(., aes(x = x, y = y)) + 
geom_point() + 
geom_smooth(method = 'lm', formula = y~x)

这里使用默认选项计算标准错误.但是,我想使用包三明治和lmtest中提供的鲁棒方差 – 协方差矩阵

也就是说,使用vcovHC(mymodel,“HC3”)

有没有办法使用geom_smooth()函数以简单的方式获得它？

《ggplot2：如何为geom_smooth中的预测获得稳健的置信区间？》

最佳答案 HC强大的SE(简单)

由于estimatr package及其lm_robust函数系列,现在可以轻松完成.例如.

library(tidyverse)
library(estimatr)

dataframe %>% 
  ggplot(aes(x = x, y = y)) + 
  geom_point() + 
  geom_smooth(method = 'lm_robust', formula = y~x, fill="#E41A1C") + ## Robust (HC) SEs
  geom_smooth(method = 'lm', formula = y~x) + ## Just for comparison
  theme_minimal()

《ggplot2：如何为geom_smooth中的预测获得稳健的置信区间？》

HAC强大的SES(更多的腿部工作)

一个警告是,估计does not仍然支持HAC(即异方差性和自相关一致)SEs a lawewey-West.但是,可以使用三明治包手动获得这些(这是原始问题无论如何都要问),然后使用geom_ribbon()进行绘图.

我会说HAC SEs对于这个特定的数据集没有多大意义,但是这里有一个如何做到这一点的例子,在相关主题上回答this excellent SO回答.

reg1 <- lm(y~x, data = dataframe)

## Generate a prediction DF
pred_df <-
  data.frame(predict(reg1, se.fit = T, interval="confidence")) %>% 
  as_tibble() 
## Clean up a little bit (optional)
colnames(pred_df) <- gsub("fit.", "", colnames(pred_df))

## Get the design matrix
X_mat <- model.matrix(reg1)

## Get HAC VCOV matrix and calculate SEs
library(sandwich)
v_hac <- NeweyWest(reg1, prewhite = F, adjust = T) ## HAC VCOV (adjusted for small data sample)
var_fit_hac <- rowSums((X_mat %*% v_hac) * X_mat)  ## Point-wise variance for predicted mean

## Add these to pred_df
pred_df <-
  pred_df %>%
  mutate(se_fit_hac = sqrt(var_fit_hac)) %>%
  mutate(
    lwr_hac = fit - qt(0.975, df=df)*se_fit_hac,
    upr_hac = fit + qt(0.975, df=df)*se_fit_hac
    )

bind_cols(
  dataframe,
  pred_df
  ) %>%
  ggplot(aes(x = x, y = y, ymin=lwr_hac, ymax=upr_hac)) + 
  geom_point() + 
  geom_ribbon(fill="#E41A1C", alpha=0.3, col=NA) + ## Robust (HAC) SEs
  geom_smooth(method = 'lm', formula = y~x) + ## Just for comparison
  theme_minimal()

《ggplot2：如何为geom_smooth中的预测获得稳健的置信区间？》

请注意,如果您愿意,也可以使用此方法手动计算并绘制其他强大的SE预测(例如HC1,HC2等).您需要做的就是使用相关的三明治估算器.例如,使用vcovHC(reg1,type =“HC2”)而不是NeweyWest(reg1,prewhite = F,adjust = T)将为使用estimatetr包的第一个示例提供相同的HC-robust CI.

点赞