Używanie klastrowej macierzy kowariancji w predict.lm()

Question

Analizuję zbiór danych, w którym dane są zgrupowane w kilku grupach (miejscowości w regionach). Zbiór danych wygląda następująco:

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 

chcę mój lm() szacuje, w celu uwzględnienia intraclass korelacji w grupach i do tego celu używam funkcji cl() że bierze lm() i zwraca solidne klastrowego macierz kowariancji (oryginalna 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) 
} 

teraz

output <- lm(y ~ x, data = df) 
clcov <- cl(output, df$groups) 
coeftest(output, clcov, nrow(df) - 1) 

daje mi szacunki mi potrzeba. Problem polega na tym, że chcę użyć modelu do prognozowania i potrzebuję, aby standardowy błąd prognozy został obliczony przy użyciu nowej macierzy kowariancji clcov. Oznacza to, że muszę

predict(output, se.fit = TRUE) 

ale stosując clcov zamiast vcov(output). Coś takiego jak vcov() <- byłoby idealne.

Oczywiście, mógłbym napisać własną funkcję do przewidywania, ale zastanawiam się, czy istnieje bardziej praktyczna metoda, która pozwala mi używać metod do podpisu lm (jak ramię :: sim).

Answer 1

Sformułowanie se.fit nie jest obliczane przy użyciu macierzy vcov, ale przy użyciu rozkładu qr i wariancji resztkowej. Dotyczy to również funkcji vcov(): pobiera nieskalowaną macierz cęgową z summary.lm() wraz z wariancją resztkową i używa tych. I nieskalowana macierz cov jest - znowu - obliczana z dekompozycji QR.

Obawiam się, że odpowiedź brzmi "nie, nie ma innej możliwości niż napisanie własnej funkcji". Naprawdę nie można ustawić macierzy vcov, ponieważ jest ona ponownie obliczana w razie potrzeby. Jednak pisanie własnej funkcji jest dość trywialne.

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

Nie użyłem funkcji predict(), aby uniknąć niepotrzebnych obliczeń. W każdym razie nie skróciłoby to zbytnio kodu.

---

Na marginesie, takie pytania są zadawane na lepsze stats.stackexchange.com

Answer 2

I zmodyfikowany powyższy kod nieznacznie się być bardziej zgodne z funkcją przewidywania - w ten sposób, że nie oczekuje się, aby wprowadzić wartości dla wynik w nowej ramce danych danych

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