在predict.lm()中使用聚类协方差矩阵

Question

我正在分析一个数据集，其中的数据被分成了几个组（地区内的城镇）。这个数据集看起来像：

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()并返回强健的聚类协方差矩阵（原始文献在这里）：

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）。

- griverorz

---

你需要更具体地说明一下。首先，要使用哪个聚类函数？为什么从lm()中输出的标准误差无效？我真的不太明白你试图做什么。很可能你需要一个更广义的模型，例如glm、glmm或gam/gamm。对于简单的lm函数的标准误差，除非你在完全不同的上下文中使用它们，否则几乎没有什么可做的了。但是这时我们需要上下文... - Joris Meys

@Joris 我编辑了问题。希望现在更清楚了。请注意，我明确避免使用 glmm 模型。 - griverorz

Answer 1

我稍微修改了上面的代码，使之与预测函数更加一致——这样你就不需要在新数据框中输入结果值。

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

- Michael

Answer 2

在predict中，se.fit并不是使用vcov矩阵计算的，而是使用QR分解和残差方差计算的。vcov()函数也是如此：它将summary.lm()中的未缩放协方差矩阵与残差方差一起使用，并使用这些值。未缩放的协方差矩阵再次是从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上询问。

- Joris Meys