[R] Logistic regression and robust standard errors

[Faradj Koliev]

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

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]> 

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.


	[[alternative HTML version deleted]]