[R] Logistic regression and robust standard errors

[Achim Zeileis]

On Fri, 1 Jul 2016, Faradj Koliev wrote:

> Dear all, 
>
>
> I use ?polr? command (library: MASS) to estimate an ordered logistic regression.
>
> My model:   summary( model<- polr(y ~ x1+x2+x3+x4+x1*x2 ,data=mydata, Hess = TRUE))
>
> But how do I get robust clustered standard errors? 
>
> I??ve tried coeftest(resA, vcov=vcovHC(resA, cluster=lipton$ID))

The vcovHC() function currently does not (yet) have a "cluster" argument. 
We are working on it but it's not finished yet.

As an alternative I include the vcovCL() function below that computes the 
usual simple clustered sandwich estimator. This can be applied to "polr" 
objects and plugged into coeftest(). So

coeftest(resA, vcov=vcovCL(resA, cluster=lipton$ID))

should work.

> and summary(a <- robcov(model,mydata$ID)).

The robcov() function does in principle what you want by I'm not sure 
whether it works with polr(). But for sure it works with lrm() from the 
"rms" package.

Hope that helps,
Z

vcovCL <- function(object, cluster = NULL, adjust = NULL)
{
   stopifnot(require("sandwich"))

   ## cluster specification
   if(is.null(cluster)) cluster <- attr(object, "cluster")
   if(is.null(cluster)) stop("no 'cluster' specification found")
   cluster <- factor(cluster)

   ## estimating functions and dimensions
   ef <- estfun(object)
   n <- NROW(ef)
   k <- NCOL(ef)
   if(n != length(cluster))
     stop("length of 'cluster' does not match number of observations")
   m <- length(levels(cluster))

   ## aggregate estimating functions by cluster and compute meat
   ef <- sapply(levels(cluster), function(i) colSums(ef[cluster == i, ,
     drop = FALSE]))
   ef <- if(NCOL(ef) > 1L) t(ef) else matrix(ef, ncol = 1L)
   mt <- crossprod(ef)/n

   ## bread
   br <- try(bread(object), silent = TRUE)
   if(inherits(br, "try-error")) br <- vcov(object) * n

   ## put together sandwich
   vc <- 1/n * (br %*% mt %*% br)

   ## adjustment
   if(is.null(adjust)) adjust <- class(object)[1L] == "lm"
   adj <- if(adjust) m/(m - 1L) * (n - 1L)/(n - k) else m/(m - 1L)

   ## return
   return(adj * vc)
}


> Neither works for me. So I wonder what am I doing wrong here?
>
>
> All suggestions are welcome ? thank you!
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.