Utilisation de la matrice de covariance en cluster dans predict.lm()

Question

J'analyse un ensemble de données dans lequel les données sont regroupées en plusieurs groupes (villes dans les régions). L'ensemble de données ressemble à:

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 

Je veux que mon estime lm() pour tenir compte de la corrélation intraclasse en groupes et à cette fin, je me sers d'une fonction cl() qui prend un lm() et retourne la matrice de covariance cluster robuste (original 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) 
} 

maintenant,

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

me donne les estimations dont j'ai besoin. Le problème est maintenant que je veux utiliser le modèle pour la prédiction, et j'ai besoin que l'erreur standard de la prédiction soit calculée avec la nouvelle matrice de covariance clcov. Qui est, je dois

predict(output, se.fit = TRUE) 

mais en utilisant clcov au lieu de vcov(output). Quelque chose comme un vcov() <- serait parfait.

Bien sûr, je pourrais écrire ma propre fonction pour faire des prédictions, mais je me demande s'il existe une méthode plus pratique qui me permet d'utiliser des méthodes pour la signature lm (comme arm :: sim).

Answer 1

Le paramètre se.fit de prédiction n'est pas calculé à l'aide de la matrice vcov, mais en utilisant la décomposition de qr et la variance résiduelle. Cela vaut aussi pour la fonction vcov(): elle prend la matrice cov non-scalée du summary.lm() avec la variance résiduelle, et utilise celles-ci. Et la matrice cov non mise à l'échelle est - encore - calculée à partir de la décomposition QR.

J'ai donc peur que la réponse soit "non, il n'y a pas d'autre choix que d'écrire sa propre fonction". Vous ne pouvez pas vraiment définir la matrice vcov, car elle est recalculée en cas de besoin. Pourtant, écrire votre propre fonction est plutôt trivial.

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

Je n'ai pas utilisé la fonction prédire() pour éviter les calculs inutiles. Cela ne raccourcirait pas le code de toute façon.

---

Sur une note de côté, des questions comme celle-ci sont mieux posées sur stats.stackexchange.com

Answer 2

J'ai modifié le code ci-dessus légèrement pour être plus conforme à la prédire la fonction - de cette façon, vous n'êtes pas censé entrer des valeurs pour le résultat dans la nouvelle dataframe

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