[R] extend summary.lm for hccm?

[Frank E Harrell Jr]

John Fox wrote:
> Dear Dirk and Ivo,
> 
> It's true that the sandwich package has more extensive facilities for
> sandwich variance estimation than the hccm function in the car package, but
> I think that the thrust of Ivo's question is to get a model summary that
> automatically uses the adjusted standard errors rather than having to
> compute them and use them "manually."  Here's one approach, which could be
> modified slightly if Ivo wants "hc0" as the default; it could also be
> adapted to use the sandwich package.
> 
> I hope this helps,
>  John

Another approach:

library(Design)  # also requires library(Hmisc)
f <- ols(y ~ x1 + x2, x=TRUE, y=TRUE)
f <- robcov(f)   # sandwich; also allows clustering.  Also see bootcov
anova(f)         # all later steps use sandwich variance matrix
summmary(f)
contrast(f, list(x1=.5), list(x1=.2))

BUT note that sandwich covariance matrix estimators can have poor mean 
squared error (a paper by Bill Gould in Stata Technical Bulletin years 
ago related to logistic regression showed an example with a 100-fold 
increase in the variance of a variance estimate) and can give you the 
right estimate of the wrong quantity (reference below).

Frank Harrell

@Article{free06so,
   author =               {Freedman, David A.},
   title =                {On the so-called ``{Huber} sandwich 
estimator'' and
``robust standard errors''},
   journal =      The American Statistician,
   year =                 2006,
   volume =               60,
   pages =                {299-302},
   annote =               {nice summary of derivation of sandwich
estimators;questions why we should be interested in getting the right
variance of the wrong parameters when the model doesn't fit}
}


> 
> ----------- snip --------------
> 
> summaryHCCM <- function(model, ...) UseMethod("summaryHCCM")
> 
> summaryHCCM.lm <- function(model, type=c("hc3", "hc0", "hc1", "hc2", "hc4"),
> 
>     ...){
>     if (!require(car)) stop("Required car package is missing.")
>     type <- match.arg(type)
>     V <- hccm(model, type=type)
>     sumry <- summary(model)
>     table <- coef(sumry)
>     table[,2] <- sqrt(diag(V))
>     table[,3] <- table[,1]/table[,2]
>     table[,4] <- 2*pt(abs(table[,3]), df.residual(model), lower.tail=FALSE)
>     sumry$coefficients <- table
>     print(sumry)
>     cat("Note: Heteroscedasticity-consistant standard errors using
> adjustment",
>         type, "\n")
>     }    
> 
> --------------------------------
> John Fox
> Department of Sociology
> McMaster University
> Hamilton, Ontario
> Canada L8S 4M4
> 905-525-9140x23604
> http://socserv.mcmaster.ca/jfox 
> -------------------------------- 
> 
>> -----Original Message-----
>> From: r-help-bounces at stat.math.ethz.ch 
>> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Dirk 
>> Eddelbuettel
>> Sent: Saturday, December 23, 2006 11:16 PM
>> To: ivo welch
>> Cc: r-help at stat.math.ethz.ch
>> Subject: Re: [R] extend summary.lm for hccm?
>>
>>
>> On 23 December 2006 at 20:46, ivo welch wrote:
>> | I wonder whether it is possible to extend the summary 
>> method for the 
>> | lm function, so that it uses an option "hccm" (well, model 
>> "hc0").  In 
>> | my line of work, it is pretty much required in reporting of 
>> almost all 
>> | linear regressions these days, which means that it would be 
>> very nice 
>> | not to have to manually library car, then sqrt the diagonal, and 
>> | recompute T-stats; instead, I would love to get everything 
>> in the same 
>> | format as the current output---except errors heteroskedasticity 
>> | adjusted.
>> | 
>> | easy or hard?
>>
>> Did you consider the 'sandwich' package? A simple
>>
>> 	> install.packages("sandwich")
>> 	> library(sandwich)
>> 	> ?vcovHC
>> 	> ?vcovHAC
>> 	
>> should get you there.
>>
>> Dirk
>>
>> --



-- 
Frank E Harrell Jr   Professor and Chair           School of Medicine
                      Department of Biostatistics   Vanderbilt University