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

Answer 1

预测中的拟合不是使用vcov矩阵计算的，而是使用qr分解和残差的。这也适用于vcov（）函数：它从summary.lm（）中获取未缩放的cov矩阵以及残差，然后使用那些残差。再次，根据QR分解计算未缩放的cov矩阵。

因此，恐怕答案是“不，除了编写自己的函数，别无选择”。您无法真正设置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上问类似这样的问题

Answer 2

我对上面的代码进行了一些修改，使其与预测函数更加一致-这样一来，您就不必在ne​​wdata数据框中输入结果的值

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