ggplot2中强大的标准错误 [英] robust standard errors in ggplot2

问题描述

我想用ggplot2绘制一个模型。我估计了一个强大的方差 - 协方差矩阵，我想在估计置信区间时使用它。

我可以告诉ggplot2使用我的VCOV，或者，也可以我以某种方式强迫predict.lm使用我的VCOV矩阵？一个虚拟的例子：

  source（http://people.su.se/~ma/clmclx.R）$ （样本（1:10,100，取代= T））=（数据帧）（x 1 = r norm（100），x 2 = r norm（100），y = ））
 lm1 < -  lm（y〜x1 + x2，data = df）
 coeftest（lm1）
 ##输出coef.test，但可以修改为输出VCOV 
 clx（lm1，1，df $ group）
  

添加到ggplot，如果我可以通过增强VCOV矩阵得到'正确'的预测结果。

只有标准错误，而不是预测，应该改变 - 对吗？

  getvcov<  -  function（fm，dfcw，cluster）{
 library三明治）;库（lmtest）
 M < - 长度（唯一（群集））
 N < - 长度（簇）
 K < -  fm $等级
 dfc （f），2，函数（x）tapply（x，cluster）（1）应用（estfun（fm） ，sum））; 
 dfc * sandwich（fm，meat = crossprod（uj）/ N）* dfcw 
} 
 
 V < -  getvcov（lm1,1，df $ group）
 X < -  as.matrix（model.frame（lm1））
 se < -  predict（lm1，se = TRUE）$ se.fit 
 se_robust<  -  sqrt（diag（X ％*％V％*％t（X）））
  

I would like to plot a model with ggplot2. I have estimated a robust variance-covariance matrix which I would like to use when estimating the confidence interval.

Can I tell ggplot2 to use my VCOV, or, alternatively, can I somehow force predict.lm to use my VCOV matrix? A dummy example:

source("http://people.su.se/~ma/clmclx.R")
df <- data.frame(x1 = rnorm(100), x2 = rnorm(100), y = rnorm(100), group = as.factor(sample(1:10, 100, replace=T))) 
lm1 <- lm(y ~ x1 + x2, data = df)
coeftest(lm1)
## outputs coef.test, but can be modified to output VCOV
clx(lm1, 1, df$group)

It would be relatively easy to add to a ggplot, if I could get 'correct' predictions given my augmented VCOV-matrix.

解决方案

Only the standard errors, not the predictions, should change -- right?

getvcov <- function(fm,dfcw,cluster) {
  library(sandwich);library(lmtest)
  M <- length(unique(cluster))   
  N <- length(cluster)           
  K <- fm$rank                        
  dfc <- (M/(M-1))*((N-1)/(N-K))  
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  dfc*sandwich(fm, meat=crossprod(uj)/N)*dfcw
}

V <- getvcov(lm1,1,df$group)
X <- as.matrix(model.frame(lm1))
se <- predict(lm1,se=TRUE)$se.fit
se_robust <- sqrt(diag(X %*% V %*% t(X)))

这篇关于ggplot2中强大的标准错误的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！