2012-02-13 15 views
3

ggplot2でモデルをプロットしたいと思います。私は、信頼区間を推定する際に使用したいと考える堅牢な分散 - 共分散行列を推定しました。ggplot2の堅牢な標準エラー

私のVCOVを使用するようにggplot2に指示することもできますし、代わりに何らかの理由で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に追加することは比較的容易であろう。

答えて

4

予測ではなく標準エラーのみが変更される必要があります。

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

これは素晴らしいです。ベンに感謝します。コアRの機能はすごく強力ですが、私はしばしば完全な可能性を利用するのに十分なほどよく分かっていません。 – Rasmus

関連する問題