我正在分析一个数据集,其中数据聚集在多个组(区域中的城镇)中。数据集如下所示:
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 http://people.su.se/~ma/clustering.pdf):
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)
}
Now,
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)。