test_glm_weights commented-out/unused

[jbrockmendel]

https://github.com/statsmodels/statsmodels/blob/master/statsmodels/genmod/tests/test_glm_weights.py#L213

    def test_basic(self):
        super(TestGlmPoissonPwNr, self).test_basic()

    @pytest.mark.xfail(reason='Known to fail')
    def test_compare_optimizers(self):
        super(TestGlmPoissonPwNr, self).test_compare_optimizers()

"Known to fail" is not a helpful error message.

https://github.com/statsmodels/statsmodels/blob/master/statsmodels/genmod/tests/test_glm_weights.py#L231

        cls.res1 = GLM(cpunish_data.endog, cpunish_data.exog,
                        family=sm.families.Poisson(), freq_weights=fweights
                        ).fit(cov_type='HC0') #, cov_kwds={'use_correction':False})

Is the cov_kwds={'use_correction': False} needed? As usual: should be deleted or explained. This shows up in a couple of places in this file.

https://github.com/statsmodels/statsmodels/blob/master/statsmodels/genmod/tests/test_glm_weights.py#L275

        # no wnobs yet in sandwich covariance calcualtion
        cls.corr_fact = 1 / np.sqrt(n_groups / (n_groups - 1))   #np.sqrt((wsum - 1.) / wsum)

Does the "no wnobs..." comment relate to the commented-out sqrt?

https://github.com/statsmodels/statsmodels/blob/master/statsmodels/genmod/tests/test_glm_weights.py#L903

    y = y_true + 2 * np.random.randn(nobs)

y defined but unused. What is it for?

https://github.com/statsmodels/statsmodels/blob/master/statsmodels/genmod/tests/test_glm_weights.py#L201
aweights is defined but unused. Should it be? This shows up a couple of times in this file.

[thequackdaddy]

In addition, I probably should have better documented this... Essentially, the comparisons came from R, and R took shortcuts for a couple families, so matching isn't possible/sensible.

        if isinstance(self, TestGlmGaussianWLS):
            # This won't work right now either
            return None
        if not isinstance(self, (TestGlmGaussianAwNr, TestGlmGammaAwNr)):
            # Matching R is hard
            assert_allclose(res1.llf, res2.ll, atol=1e-6, rtol=1e-7)

To answer your other questions, "Known to fail" might be replaced with "prob_weights not yet implemented". The cov_kwds={'use_correction': False} is only meaningful with cluster covariance. There was a lot of copypasta making these cases, so IMHO there are extraneous things that could be removed harmlessly.

[jbrockmendel]

@thequackdaddy thanks for helping clear this up. Double-checking to make sure all the bases are covered:

• "known to fail" I'll replace.

• commented out occurrences of #,cov_kwds={'use_correction': False} in TestGlmPoissonFwHC and TestGlmPoissonAwHC I'll delete.

• delete unused aweights in TestGlmPoissonFwClu, TestGlmPoissonFwHC, TestGlmPoissonPwNr. (pls confirm this one)

• what about the #np.sqrt((wsum - 1.) / wsum) in TestGlmPoissonFwClu?

• test_poisson_residuals has a y that is defined but never used. Any thoughts there?

Are you also a good person to ask about test_glm or other test files (results_glm would be a big win)?

[thequackdaddy]

• aweights - I think its safe to delete... but Josef is the ultimate arbitrator.

• I think the #np.sqrt((wsum - 1.) / wsum) was something we needed in some families/cov_params to match stata. I presume it could be removed.

• test_poisson_resids and y ... the way that looks, my guess is it was related to the binomial family... not sure if it its necessary.

• test_glm - I've edited a few parts of this when we were first adding weights and a few other items. I might know something...

• results_glm that's older than Methuselah... Probably legitimately before linting was as well-defined as it is now. I'm gonna walk away from that one...

[jbrockmendel]

I'm gonna walk away from that one...

Totally fair, you've done your good deed for the day.

@josef-pkt the last unresolved commented-out stuff in results_glm is within reach:
https://github.com/statsmodels/statsmodels/blob/master/statsmodels/genmod/tests/results/results_glm.py#L1151

class InvGauss(object):
    '''
    Usef
    Data was generated by Hardin and Hilbe using Stata.
    Note only the first 5000 observations are used because
    the models code currently uses np.eye.
    '''
#        np.random.seed(54321)
#        x1 = np.abs(stats.norm.ppf((np.random.random(5000))))
#        x2 = np.abs(stats.norm.ppf((np.random.random(5000))))
#        X = np.column_stack((x1,x2))
#        X = add_constant(X)
#        params = np.array([.5, -.25, 1])
#        eta = np.dot(X, params)
#        mu = 1/np.sqrt(eta)
#        sigma = .5
#       This isn't correct.  Errors need to be normally distributed
#       But Y needs to be Inverse Gaussian, so we could build it up
#       by throwing out data?
#       Refs: Lai (2009) Generating inverse Gaussian random variates by
#        approximation
# Atkinson (1982) The simulation of generalized inverse gaussian and
#        hyperbolic random variables seems to be the canonical ref
#        Y = np.dot(X,params) + np.random.wald(mu, sigma, 1000)
#        model = GLM(Y, X, family=models.family.InverseGaussian(link=\
#            models.family.links.identity))

    def __init__(self):

It looks like there might be a couple of usable references in there, but most of this is adding more confusion than its resolving. Thoughts?

[josef-pkt]

This are comments from Skipper's initial implementation.
I never even looked at or ran an example of invgauss.

This is part of the maintainer and developer documentation until inverse gaussian gets a proper review. Until then, I wouldn't have any idea about which information is useful and which is not.

Similarly, gamma has numerical or convergence problems with the canonical link, but I never looked closely enough to learn anything about it. (gamma works fine with exp link which does not run into corner cases at zero, so figuring out the canonical inverse link is indefinitely postponed.)

[jbrockmendel]

@jseabold Can you help clear up this last-unsolved-mystery in glm_results?

@josef-pkt It's totally reasonable to put off dealing with this until the right person shows up, but this is distinctly not documentation.

[josef-pkt]

@jbrockmendel It's maintainer and developer documentation! It's distinctly not user documentation and it's not part of it.
I never looked at how invgauss random variables are generated, but it provides two references and a buggy (?) code that just needs to be fixed or completed.
Or check scipy.stats.distributions which just uses self._random_state.wald(mu, 1.0, size=self._size)
Wikipedia mentions the alias "inverse Gaussian distribution (also known as the Wald distribution)"

So it just takes a bit of effort to figure this out, if somebody takes the time.

[jbrockmendel]

The word you're looking for is "technical debt". If it a) explains/documents nothing and b) needs to be explained/documented, then it is not documentation of any variety.

I acknowledged the paper references as being usable. Its the commented-out code and the "This isn't correct" that are unhelpful.

[josef-pkt]

so I guess the corrected code should be something like
Y = np.random.wald(np.dot(X,params), sigma, 1000)
but I don't know what sigma is supposed to be,
var(y) = mu**3 for invgauss
and scale would be a dispersion parameter a, like var(y) = a * mu**3
and scale for scipy.stats.distribution is sqrt of our scale, i.e. the parameterization for get_distribution might correspond to Y = np.random.wald(np.dot(X, res.params), np.sqrt(res.scale), 1000)
(needs to be confirmed)

[josef-pkt]

"technical debt" is ok, it's another term for "agile software development". :)
https://en.wikipedia.org/wiki/Agile_software_development

(
with the caveat that the delay between iterations can be many years, given the size of our developer and reviewer team
"Most agile development methods break product development work into small increments that minimize the amount of up-front planning and design. Iterations, or sprints, are short time frames (timeboxes) that typically last from one to four weeks."
substitute: from one to four weeks -> from weeks to many years
)

[jbrockmendel]

so I guess the corrected code should be something like

I'm not convinced "the corrected code" needs to exist. I think its more likely this was an older attempt at producing results for glm_results that is superceded by a newer version, should be deleted entirely. I'll wait for skipper to chime in.

[jbrockmendel]

"technical debt" is ok, it's another term for "agile software development". :)

What do you think that word means? I don't think it means what you think it means.
From the Wikipedia link:

Working software is delivered frequently (weeks rather than months)
[...]
Continuous attention to technical excellence and good design
[...]
<b>Allowing technical debt to build up</b>

Focusing on delivering new functionality may result in increased technical debt. The team must allow themselves time for defect remediation and refactoring. Technical debt hinders planning abilities by increasing the amount of unscheduled work as production defects distract the team from further progress.[97]

As the system evolves it is important to refactor as entropy of the system naturally increases.[98] Over time the lack of constant maintenance causes increasing defects and development costs.[97]

And from https://www.atlassian.com/agile/software-development/technical-debt

If a feature isn't up to snuff, it doesn't ship.

Ring any bells?

[josef-pkt]

I think both wikipedia and the technical-debt link describe pretty well what we are doing, except for timing and formal organization which we lack because we are not commercial, i.e. we have volunteer work, lack of maintainer wo/man power and turnover in "team membership".

"Focusing on delivering new functionality may result in increased technical debt. The team must allow themselves time for defect remediation and refactoring."
means you allow technical debt to get into the code base but you also allocate time for reducing technical debt and for refactoring in each iteration.

The atlassian link is good: each () feature gets tested immediately and we refine and enhance in iterations.
() theoretically, there are always some gaps

"Technical debt is the difference between what was promised and what was actually delivered."

I think we are pretty good in not overselling things. Both freq_weights and var_weights have very good test coverage at least for all the basic results, But some things are not yet adjusted to handle those weights, and we say so. Prob_weights are still WIP, and it's not clear yet how various results are supposed to be defined in that case (i.e. I don't like how Stata defines things for pweights).
Our number of serious bugs is pretty low, but not zero.
Most "bugs" are more like feature requests to improve the API or handle additional cases or handle special cases. (It's a feature as long as we say so; after we change it, it would be a bug.)

I don't like rigorous adherence to any philosophy, whether software development or statistics or economics, but I think our software development is pretty test driven and agile.

To the invgauss case:
Next iteration: implement get_distribution for it, then the "junk",, I meant hints for future developers, can go into the git history.