[R] extend summary.lm for hccm?

[Frank E Harrell Jr]

John Fox wrote:
> Dear Frank,
> 
> If I remember Freedman's recent paper correctly, he argues that sandwich
> variance estimator, though problematic in general, is not problematic in the
> case that White described -- an otherwise correctly specified linear model
> with heteroscedasticity estimated by least-squares. 
> 
> Regards,
>  John

That's right John although the precision of the variance estimators is 
still worth studying in that case.

Best regards,
Frank

> 
> --------------------------------
> 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 Frank 
>> E Harrell Jr
>> Sent: Sunday, December 24, 2006 10:29 AM
>> To: John Fox
>> Cc: r-help at stat.math.ethz.ch; 'ivo welch'; 'Dirk Eddelbuettel'
>> Subject: Re: [R] extend summary.lm for hccm?
>>
>> 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