[R] Robust SE for lrm object

[Frank E Harrell Jr]

Achim Zeileis wrote:
> On Sat, 6 Mar 2010, David Winsemius wrote:
> 
>> On Mar 5, 2010, at 11:54 PM, Patrick Shea wrote:
>>
>>>
>>> I'm trying to obtain the robust standard errors for a multinomial 
>>> ordered logit model:
>>>
>>> mod6 <- lrm(wdlshea ~   initdesch + concap + capasst + qualrat + 
>>> terrain,data=full2)
>>>
>>> The model is fine but when I try to get the RSE I get an error.
>>>
>>> coeftest(mod6, vcov = vcovHAC(mod6))
>>>
>>> Error in match.arg(type) :
>>> 'arg' should be one of ?ordinary?, ?score?, ?score.binary?,  
>>> ?pearson?, ?deviance?, ?pseudo.dep?, ?partial?, ........etc.
>>>
>>> I'm a novice R user and am not sure how to address this problem. I 
>>> have also tried to use alternatives   (zelig, polr) but have had no 
>>> luck. Any assistance on generating RSE for a multinomial order logit 
>>> model would be appreciated
>>
>> Have you loaded the library that contains the vcovHAC function?
> 
> That is in the "sandwich" package. However, I doubt that it makes sense 
> in this context. Using HAC covariances would imply that you have time 
> series data and want to correct for heteroskedasticity and 
> autocorrelation. I'm not even sure whether sandwich standard errors 
> would be terribly useful. Both would require that you correctly 
> specified the estimating functions of your proportional odds logistic 
> regression but misspecified a few other aspects of the remaining 
> likelihood. Not sure whether that can be obtained for an ordinal 
> multinomial response.
> 
>> (And do you know whether coeftest works with Design/rms objects?)
> 
> It does (unlike its own summary function in some situations):
> 
> library("rms")
> library("lmtest")
> data("BankWages", package = "AER")
> fm <- lrm(job ~ ., data = BankWages)
> summary(fm)
> coeftest(fm)
> 
> The reason why vcovHAC() or sandwich() do not work is that bread() and 
> estfun() methods would need to be available for "lrm" objects which is 
> currently not the case (dito for "polr" objects). In principle they 
> could be written, see
>   vignette("sandwich-OOP", package = "sandwich")
> but as I said above I'm not sure whether it would be very useful.
> Z

Excellent posts.  I'll just add that the Huber-White sandwhich estimator 
is obtain in the rms and Design packages by using the robcov function on 
a fit object.  You can also use the bootstrap with bootcov.

Frank

> 
> 
>> -- 
>>
>> David Winsemius, MD
>> West Hartford, CT
>>
>> ______________________________________________


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