在predict.lm（）中使用聚類協方差矩陣

我正在分析一個數據集，其中數據聚集在幾個組中（城鎮中的區域）。該數據集的樣子：在predict.lm（）中使用聚類協方差矩陣

R> df <- data.frame(x = rnorm(10), 
        y = 3*rnorm(x), 
        groups = factor(sample(c('0','1'), 10, TRUE))) 
R> head(df) 
     x  y groups 
1 -0.8959 1.54  1 
2 -0.1008 -2.73  1 
3 0.4406 0.44  0 
4 0.0683 1.62  1 
5 -0.0037 -0.20  1 
6 -0.8966 -2.34  0

我希望我的LM（）估計佔同類相關的羣體，併爲此我使用的是功能cl()接受一個lm()並返回可靠的集羣協方差矩陣（原here ）：

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

現在，

output <- lm(y ~ x, data = df) 
clcov <- cl(output, df$groups) 
coeftest(output, clcov, nrow(df) - 1)

給我我需要的估計。現在的問題是我想使用該模型進行預測，並且需要用新的協方差矩陣clcov來計算預測的標準誤差。也就是說，我需要

predict(output, se.fit = TRUE)

但使用clcov，而不是vcov(output)。像vcov() <-就像是完美的。

當然，我可以編寫自己的函數來做預測，但我只是想知道是否有一個更實用的方法，允許我使用方法簽名lm（如arm :: sim）。

來源

2010-09-24 griverorz

您需要指定更多一點。那個集羣函數是以什麼開始的？爲什麼從lm（）中產生的標準錯誤無效？我無法真正跟隨你想要做的事情。很可能你需要一個更廣義的模型，例如glm，glmm或gam/gamm。除非在完全不同的環境中使用它們，否則對於簡單的lm函數的標準錯誤幾乎沒有做任何事情。但後來我們需要上下文... – 2010-09-24 23:16:29

@Joris我編輯了這個問題。希望現在更清楚。請注意，我明確避免使用「glmm」模型。 – griverorz 2010-09-25 00:12:47

預測中的se.fit不是使用vcov矩陣計算的，而是使用qr分解和殘差方差。這同樣適用於vcov（）函數：它將summary.lm（）中的非標度cov矩陣與剩餘方差一起使用，並使用這些矩陣。並且未分級的cov矩陣 - 再次從QR分解中計算。

所以恐怕答案是「不，除了編寫自己的功能沒有別的選擇」。您無法真正設置vcov矩陣，因爲它在需要時會被重新計算。然而，編寫自己的函數是相當微不足道的。

predict.rob <- function(x,clcov,newdata){ 
    if(missing(newdata)){ newdata <- x$model } 
    m.mat <- model.matrix(x$terms,data=newdata) 
    m.coef <- x$coef 
    fit <- as.vector(m.mat %*% x$coef) 
    se.fit <- sqrt(diag(m.mat%*%clcov%*%t(m.mat))) 
    return(list(fit=fit,se.fit=se.fit)) 
}

我沒有使用predict（）函數來避免不必要的計算。無論如何，它不會縮短代碼。

在一個側面說明，類似這樣的問題，更好地問上stats.stackexchange.com

來源

2010-09-25 13:15:46

我修改了上面的代碼稍微要與預測功能更加一致的 - 這樣你是不是有望進入值新數據數據幀中的結果

predict.rob <- function(x,clcov,newdata){ 
if(missing(newdata)){ newdata <- x$model } 
tt <- terms(x) 
Terms <- delete.response(tt) 
m.mat <- model.matrix(Terms,data=newdata) 
m.coef <- x$coef 
fit <- as.vector(m.mat %*% x$coef) 
se.fit <- sqrt(diag(m.mat%*%clcov%*%t(m.mat))) 
return(list(fit=fit,se.fit=se.fit))}

來源

2012-01-27 21:18:32 Michael

在predict.lm（）中使用聚類協方差矩陣

回答

相關問題