ggplot2 中稳健的标准错误

发布于 2025-01-04 16:06:52 字数 513 浏览 1 评论 0原文

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

我可以告诉 ggplot2 使用我的 VCOV，或者，我可以以某种方式强制 Predict.lm 使用我的 VCOV 矩阵吗？一个虚拟示例：

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)

如果我能够根据增强的 VCOV 矩阵获得“正确”的预测，那么添加到 ggplot 中将相对容易。

原文

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.

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

沧笙踏歌 2025-01-11 16:06:52

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

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

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

回复收藏 0 原文

~没有更多了~