Надежные стандартные ошибки в ggplot2
Я хотел бы построить модель с ggplot2. Я оценил надежную матрицу дисперсии-ковариации, которую я хотел бы использовать при оценке доверительного интервала.
Могу ли я сказать ggplot2 использовать мой VCOV или, альтернативно, могу ли я каким-то образом принудительно заставить предикат.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)
Было бы относительно легко добавить к ggplot, если бы я мог получить "правильные" прогнозы, учитывая мою расширенную VCOV-матрицу.
1 ответ
Решение
Должны измениться только стандартные ошибки, а не прогнозы - верно?
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)))