Verwenden der Cluster-Kovarianzmatrix in predict.lm()

Question

Ich analysiere ein Dataset, in dem Daten in mehreren Gruppen (Städte in Regionen) gruppiert sind. Der Datensatz wie folgt aussieht:

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 

ich meine lm() will für Intra Korrelation in Gruppen zu berücksichtigen schätzt, und zu diesem Zweck verwende ich eine Funktion cl(), die ein lm() und gibt die robuste gruppierten Kovarianzmatrix (original here nimmt):

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

Jetzt

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

gibt mir die Schätzungen ich brauche. Das Problem ist nun, dass ich das Modell für die Vorhersage verwenden möchte und dass der Standardfehler der Vorhersage mit der neuen Kovarianzmatrix clcov berechnet werden muss. Das heißt, ich brauche

predict(output, se.fit = TRUE) 

aber mit clcov statt vcov(output). Etwas wie ein vcov() <- wäre perfekt.

Natürlich könnte ich meine eigene Funktion schreiben, um Vorhersagen zu machen, aber ich frage mich nur, ob es eine praktischere Methode gibt, die mir Methoden zur Signatur lm (wie arm :: sim) erlaubt.

Answer 1

Die se.fit in Vorhersage wird nicht mit der vcov-Matrix berechnet, sondern mit der qr-Zerlegung und der Restvarianz. Dies gilt auch für die Funktion vcov(): Sie nimmt die unskalierte Cov-Matrix aus summary.lm() zusammen mit der Restvarianz und verwendet diese. Und die unskalierte cov-Matrix wird - wiederum - aus der QR-Zerlegung berechnet.

Also ich fürchte die Antwort ist "Nein, es gibt keine andere Möglichkeit, als Ihre eigene Funktion zu schreiben". Sie können die VCOV-Matrix nicht wirklich festlegen, da sie bei Bedarf neu berechnet wird. Ihre eigene Funktion zu schreiben ist jedoch eher trivial.

Ich habe die Funktion predict() nicht verwendet, um unnötige Berechnungen zu vermeiden. Es würde den Code sowieso nicht zu sehr verkürzen.

---

Auf einer seitlichen Anmerkung, sind wie diese Fragen besser auf stats.stackexchange.com fragte

Answer 2

ich den obigen Code geändert etwas mehr im Einklang mit der vorherzusagen Funktion zu sein - Werte auf diese Weise Sie sind nicht zu erwarten eingeben, für das Ergebnis im newdata Datenrahmen

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