[R] Logistic regression and robust standard errors

[Achim Zeileis]

On Fri, 1 Jul 2016, Faradj Koliev wrote:

> Dear Achim Zeileis, 
> Many thanks for your quick and informative answer. 
> 
> I?m sure that the vcovCL should work, however, I experience some problems. 
> 
> 
> > coeftest(model, vcov=vcovCL(model, cluster=mydata$ID))
> 
> First I got this error: 
> 
> Error in vcovCL(model, cluster = mydata$ID) : 
>   length of 'cluster' does not match number of observations

This probably means that not all rows of "mydata" have been used for the 
estimation of the model, e.g., due to missing values or something like 
that. Hence, the mismatch seems to occur.

> After checking the observations I got this error: 
> 
> Error in vcovCL(model, cluster = mydata$ID) : object 'tets' not found
> Called from: vcovCL(model, cluster = mydata$ID)
> Browse[1]> 

The variable "tets" is presumably one of the regressors in your "model" 
and apparently it cannot be found when extracting the model scores. Maybe 
you haven't stored the model frame and deleted the data.

Hard to say without a simple reproducible example.
Z

> What can I do to fix it? What am I doing wrong now? 
> 
> 
> 
> 
>
>       1 jul 2016 kl. 14:57 skrev Achim Zeileis
>       <Achim.Zeileis at uibk.ac.at>:
> 
> 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.
> 
> 
> 
>