Clustered Standard Errors for MLE Coefficients

[Damien Moore]

Hi,

I had a need to generate clustered standard errors for a paper that I have been working on and saw the note here: https://github.com/statsmodels/statsmodels/issues/863

It didn't look like anyone had done this work so I took a crack at it for the MLE case following the reference implementation in the CGM paper.  The implementing code is in the file clustered_se.py here: https://github.com/spillz/sci-comp/tree/master/statsmodels-fun

The key functions are:

clustered_se_from_model 

multiway_clustered_se_from_model

The functions take as arguments a maximum likelihood-based model, the estimated parameters and the clustering variable(s). They should work if the model has at least defined loglike and loglikeobs (information and scoreobs would also be nice!!) -- unfortunately only some of them do.

test_clusered_se.py has examples that could be turned into tests (one of the tests replicates estimates computed by Mitch Petersen that have been tested on R, Stata and Sas).

Cheers,

dm

[josef...@gmail.com]

Good you are working on and interested in this.

I need to read your code soon in detail, just a quick look so far.

I don't know yet how to get robust (patterned) covariances for MLE, so this is very helpful.

Overall, I think it's better to reuse generic code

https://github.com/statsmodels/statsmodels/blob/master/statsmodels/stats/sandwich_covariance.py

and make it more user friendly.

For example, I just used those function for HAC Newey West for GMM moment conditions (it uses only the inner S, the sandwich itself is in GMM.)

I haven't looked at incorporating cluster robust standard errors, but it's next on my GMM todo list (plus comparison with GEE todo)

import statsmodels.stats.sandwich_covariance as smcov

       elif method == 'hac':

            maxlag = wargs['maxlag']

            if 'kernel' in wargs:

                weights_func = wargs['kernel']

            else:

                weights_func = smcov.weights_bartlett

                wargs['kernel'] = weights_func

            

            moms_centered = moms - moms.mean()

            w = smcov.S_hac_simple(moms_centered, nlags=maxlag, 

                                   weights_func=weights_func)

I also recently opened https://github.com/statsmodels/statsmodels/issues/1099

for other robust cov_params

for example cluster HAC for large T (which is already in robust_covariance), 

Thank you,

Josef

 

Cheers,

dm

[Damien Moore]

Some comments inserted below...

On Tuesday, October 8, 2013 4:05:44 PM UTC-4, josefpktd wrote:

Overall, I think it's better to reuse generic code

https://github.com/statsmodels/statsmodels/blob/master/statsmodels/stats/sandwich_covariance.py

and make it more user friendly.

For example, I just used those function for HAC Newey West for GMM moment conditions (it uses only the inner S, the sandwich itself is in GMM.)

I haven't looked at incorporating cluster robust standard errors, but it's next on my GMM todo list (plus comparison with GEE todo)

I'm not sure what you have in mind here. Feel free to use or adapt what I have there as you see fit.

In a sense the equations in the code I provided are very general because they support any model that can be expressed as a set of MLE equations with additively separable log likelihood terms for each observation. (CGM provide an even more general approach for m-estimators, encompassing both non-linear and linear. but additively separable, GMM and MLE models.) 

There are probably a variety of other diagnostics that apply to all MLE estimators (or all m-estimators) that could be applied across a wide variety of models if the models provide gradient(by obs)+hessian or score(by obs)+information matrix. For example, if the OLS model, which does derive from Likelihood model, implemented methods for the score (by obs) and information, then the clustered standard errors could be computed very easily and efficiently from the likelihood (or m)-based standard error formulas. 

I also recently opened https://github.com/statsmodels/statsmodels/issues/1099

for other robust cov_params

for example cluster HAC for large T (which is already in robust_covariance), 

 

That looks pretty interesting (I'm not familiar with the paper). There's an interesting design problem here in that you have these more general classes of models that nest a bunch of interesting special cases (e.g. m-class estimators > MLE estimators > Generalized Linear models > OLS) and in theory you would probably want to be able to supply diagnostics that apply to the widest class of models. So, for example, any model that can be expressed as MLE, you should be able to get the standard set of model tests (Wald, LR, LM, AIC? ...). For more specialized models there will be more efficient ways to calculate some of the needed quantities (e.g. standard closed forms for the score and info matrices for OLS) or more model-specific tests. 

Anyway, I guess it would be good to make sure that all models correctly derive from the more general model classes and supply as many of the base class methods as possible. That's not currently the case (I got many NotImplementedErrors in trying the get my code working)

Thanks,

dm

[josef...@gmail.com]

Just replying to this part for now.

What you are describing with the design problem is exactly what we are trying to get at.

However, we have two main limitations to get to this:

One is to figure out to which class or model a statistical result, like hypothesis tests or robust standard error, applies.

For example, I only recently found some articles that argued that we can use robust sandwich covariance matrices after WLS or GLS. 

I'm currently looking into what diagnostic tests can be used after instrumental variables estimation (IV2SLS). We have some tickets for a generic hat matrix, so we could generalize influence and outlier diagnostics to other models than OLS.

Without a good reference or an expert in the area it is difficult to figure this out, and then we are often conservative and don't apply it to general model class. The additional part is that we should write unit test that verify that those extensions or generalizations apply and are correct. (*)

The second reason is that we **make it up** as we add more models. The top level classes currently need a restructuring.

We don't have yet a generic M-Estimator class. In my current rewrite, IV2SLS subclasses LikelihoodModel and GMM results subclasses LikelihoodModelResults even though they don't have a log-likelihood, but that's where the goodies are. QuantileRegression is misclassified as RegressionModel even though it's just an M-estimator.

It would be great if we could find a classification tree for all possible models that show up in statistics and econometrics.

Overall our inheritance works great. When we subclass the "right" model, we suddenly get a lot of additional functionality without having to code anything additional. And we will get more of the "right" model classes.

(*) example: case that I didn't realize before writing unit tests for IV2SLS: 

https://github.com/statsmodels/statsmodels/issues/1104 

calculating fvalue in regressionmodel is not "generic" and uses a shortcut formula based on sum of squares instead of proper Wald test.

Josef

 

Thanks,

dm