I'm trying to get clustered standard errors on an ordered logit regression using the POLR function from MASS. The regression has about 2,900 fixed effects, (U.S. counties and years.)
When I run the regression without the fixed effects, I'm able to use the vcovCL function from sandwich with no problems, (using the procedure outlined in this answer.) However, when I include the fixed effects, I get this error:
Error in xmat %*% x$coefficients : non-conformable arguments
I would post the data, but the data set is quite large and the regression takes about 2 days to run with all the fixed effects.
Has anyone come across a similar problem and have you come up with a solution?
EDIT: Following what @ben-bolker commented below, I tried the regression with a random sample of 1,000 lines. It appears that some of the fixed effects produce a square matrix throwing the following errors:
1: glm.fit: fitted probabilities numerically 0 or 1 occurred
2: design appears to be rank-deficient, so dropping some coefs
When I try using the vcovCL function to cluster the standard errors, I get the same "non-conformable arguments" error. So I suspect that since some coefficients are NA's, the matrix vcovCL is a different size from the POLR object and that's what's causing the error.
Is there a way to excise these missing fixed-effects? I won't be displaying them in the regression output anyhow.