[R] Robust SE for lrm object

Achim Zeileis Achim.Zeileis at uibk.ac.at
Sat Mar 6 07:22:08 CET 2010


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


> -- 
>
> David Winsemius, MD
> West Hartford, CT
>
> ______________________________________________
> R-help at r-project.org 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