Uso de la matriz de covarianzas agrupadas en predict.lm()

Question

Estoy analizando un conjunto de datos en el que los datos se agrupan en varios grupos (ciudades en regiones). El conjunto de datos se parece a:

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 

Quiero que mi lm() calcula a cuenta de la correlación intraclase en grupos y para ello estoy usando una función cl() que toma un lm() y devuelve la matriz de covarianza agrupada robusta (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) 
} 

Ahora,

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

me da las estimaciones que necesito. El problema ahora es que quiero usar el modelo para la predicción, y necesito que el error estándar de la predicción se calcule con la nueva matriz de covarianza clcov. Es decir, necesito

predict(output, se.fit = TRUE) 

pero utilizando en lugar de clcovvcov(output). Algo así como un vcov() <- sería perfecto.

Por supuesto, podría escribir mi propia función para hacer predicciones, pero me pregunto si hay un método más práctico que me permita usar métodos para la firma lm (como arm :: sim).

Answer 1

El se.fit in predict no se calcula utilizando la matriz vcov, pero utilizando la descomposición qr y la varianza residual. Esto también se aplica a la función vcov(): toma la matriz de cov sin escalar del summary.lm() junto con la varianza residual, y usa esos. Y la matriz de CV no escalada se calcula, nuevamente, a partir de la descomposición QR.

Me temo que la respuesta es "no, no hay otra opción que escribir su propia función". Realmente no puede establecer la matriz vcov, ya que se vuelve a calcular cuando sea necesario. Sin embargo, escribir tu propia función es bastante 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)) 
} 

no hice uso de la función de predecir() para evitar cálculos innecesarios. No acortaría el código demasiado de todos modos.

---

En una nota lateral, este tipo de preguntas se les pide mejor en stats.stackexchange.com

Answer 2

he modificado el código de arriba ligeramente para ser más coherente con la función de predecir - de esta forma no se espera que introduzca valores para el resultado en el dataframe 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))}