Использование кластеризованной ковариационной матрицы в предикате.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() и возвращает устойчивую кластеризованную ковариационную матрицу (оригинал здесь * 1007) *):

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

Я немного изменил приведенный выше код, чтобы он больше соответствовал функции прогнозирования - таким образом, вы не должны вводить значения для результата в фрейме данных newdata

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

Answer 2

Значение se.fit в прогнозировании рассчитывается не с использованием матрицы vcov, а с использованием разложения qr и остаточной дисперсии. Это относится и к функции vcov (): она берет немасштабированную матрицу cov из 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))
}

Я не использовал функцию предсказания (), чтобы избежать ненужных вычислений. В любом случае, это не слишком укоротит код.

---

Кстати, такие вопросы лучше задать на stats.stackexchange.com