[R] extend summary.lm for hccm?

John Fox jfox at mcmaster.ca
Sun Dec 24 17:30:57 CET 2006


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

--------------------------------
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
> >>
> >> --
> 
> 
> 
> -- 
> Frank E Harrell Jr   Professor and Chair           School of Medicine
>                       Department of Biostatistics   
> Vanderbilt University
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.



More information about the R-help mailing list