Allow for custon `vcov` matrices

[pdeffebach]

trafficstars

Currently, RegressionTables.jl doesn't support custom covariance matrices. This is a pain, because GLM doesn't allow the user to cluster errors or your heteroskedastic errors in the glm command itself. People have to use CovarianceMatrices.jl to do this

It would be nice if the regtable command supported adding vcov as a keyword argument to make this work.

[nilshg]

This discussion was had a while ago on Discourse, including an explanation of how to do this here: https://discourse.julialang.org/t/ann-regressiontables-jl-produces-publication-quality-regression-tables/7516/33

The problem with the approach outlined there I think is that one can't just build up a linear model in GLM like that, as vcov is defined as

vcov(x::LinPredModel) = rmul!(invchol(x.pp), dispersion(x, true))

So afaict one would have to adapt the pp field of the LinPredModel in a way to get the desired standard errors (and I believe changing this would also update p-values and confidence intervals, although I haven't checked this in detail!). Doing this manually is probably asking a bit too much of users, but I'm also not sure a printing/formatting package like RegressionTables is the right place for this - maybe CovarianceMatrices should provide convenience functions to take model objects and return them with an updated vcov matrix?

[pdeffebach]

I would just imagine a keyword argument that accepts a matrix. Wouldn't that solve the problem?

I think that, to really solve this issue, GLM should allow for rubust SEs and clustering at the @formula level. But people are resistant to this because there are so many ways to go about calculating standard errors I think, and because of the notion that the betas are what really defines a model and not the SEs.

[nilshg]

But wouldn't you then have to implement p-values, standard errors, and confidence intervals as well, as you can't query them from the underlying model but need to recompute them with the given covariance matrix? Seems a bit heavy for a printing package?

[jmboehm]

I'm inclined to agree with @nilshg. I see the package as printing statistics from estimated models that you pass as arguments. I think passing extra arguments with information that the model type does not contain is not ideal (ideally, the model should contain it), but a compromise that I can live with for a few things (labels etc). But things that are really at the core of the model statistics, like covariance matrices and related statistics, should be part of the model and not some extra kwargs.

Of course, I agree that the issue needs to be addressed; I just feel that patching things up here by adding kwargs is not the "right" solution. But I'm open to suggestions.

I should say that I've been meaning to expand the RegressionTables interface to take arguments that are <: RegressionModel (but haven't gotten around to this due to lack of time). That would allow CovarianceMatrices.jl to define its own type and extend the methods from StatsBase.

[pdeffebach]

I hadn't thought about p-values and t-stats. You are right, that would not be an easy solution.

Maybe the best solution is a really good, well tested, clean PR to lm in GLM to get the ball rolling.

[IljaK91]

@nilshg , it is not obvious to me at the moment, is it necessary to adjust the pp object, or is it necessary to adjust the vcov function? It seems to me that the former requires a different approach than (i.e. adjusting the residuals) than computing the vcov matrix through CovarianceMatrices.jl and then substituting the vcov matrix in the regression object.

[jariji]

There is now also Vcov.jl which provides simple, robust, and cluster errors for some models.