在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（）中使用聚类协方差矩阵

回答

相关问题