WIP: fit_constrained #1714

Merged
merged 24 commits into from Jun 5, 2014

Projects

None yet

2 participants

@josef-pkt
Member

add constrained fitting to models

starting with generic setup and Poisson as example.
OLS was test example for transformation, but has separate RestrictedGLS implementation, PR #1673

  • needs offset for all models except linear, see #1640
    currently only count models and GLM have offset
  • add method to models
    • Poisson
    • other discrete, NegativeBinomial has offset but an extra parameter that needs to be handled
    • GLM
    • ? linear models, has other implementation, that avoids auxiliary model creation, delegate
    • RLM I'm pretty sure it works almost out of the box, but haven't tried yet,
      using correction to endog instead of offset is not implemented in the function
      (I have a separate example code for that, linear case is more related to RestrictedGLS PR)
  • does not return results instance, currently returns params and bse, calculates cov_params
    Those two need to be overwritten in the creation of the results instance.
    I would prefer generic results creation see #1615, but will start with a full test case for Poisson.
  • unit test currently only verify params and bse, tvalues, pvalues, confint, llf, and llnull (if available) are verified
  • start_params for transformed model: I got convergence failures sometimes in bfgs, sometimes in newton. I didn't get them in both at the same time for the poisson test cases. Needed to add start_params for TravisCI to converge in test cases 1a and 2a
  • still unclear (waiting for results instance to inherit all results to be verified):
    • df_resid, df_model, small sample correction, I assume those need to be adjusted unless we only use asymptotic normal version, no use_t=True
    • aic, bic, ... ? not yet in results module, works without changes but see #1733
  • problem: functions or results methods that want to access score_obs, hessian, get the wrong method, doesn't take constraints into account, ???
    maybe no problem: Likelihood function doesn't change, we just evaluate it at the constrained parameters.
    need to check what kshedden did for score test in GEE.
  • summary, other
    • summary shows converged=False from the wrong fit, need to adjust mle_fit results
    • summary doesn't mention constrained estimation

code duplication:
base.TransformRestriction has the same purpose as Kerby's version in GEE. I haven't checked if the result is identical.

  • API decisions, naming, ...
  • bonus question: Can I get the robust covariances just from the transformed model without any extra work? Need to be wired up.

Unit tests
I have added more test results from Stata that haven't been use so far, unconstraint, and constraint with vce(robust) option
I don't know what vce(robust) does in Stata, and I haven't figured out yet how to add sandwich robust covariances to GLM. Stata also allows HAC and cluster robust vce in GLM.

follwow-up

  • extension to other classes
  • clean up and find pattern and naming for Transformation classes (missing get_offset(exog) method)
  • unused transform_params_constraint that obtains restricted parameters from the unrestricted parameters (approximate in non-linear models, similar to Wald test)
  • development script examples/try_fit_constrained.py contains example with application to linear model (OLS) and example for transform_params_constraint
@coveralls

Coverage Status

Coverage increased (+0.05%) when pulling 7cff3eb on josef-pkt:fit_constrained into e4326b7 on statsmodels:master.

@coveralls

Coverage Status

Coverage increased (+0.06%) when pulling b62b36b on josef-pkt:fit_constrained into e4326b7 on statsmodels:master.

@josef-pkt
Member

test failures: python 3: print in junk function (draft for another example based on randhi)

TestPoissonConstrained 1a and 2a look like convergence failures on one python 2.7
still need start_params

@josef-pkt
Member

The last commit added start_params, but used it only in one of the two non-converging test cases.
Still need to figure out why two examples don't converge on TravisCI. Tests pass on my current version of dependencies, but optimizer was chosen so they pass for me.

@coveralls

Coverage Status

Coverage increased (+0.07%) when pulling b583394 on josef-pkt:fit_constrained into e4326b7 on statsmodels:master.

@josef-pkt
Member

no StringIO import on python 3
setting the start_params for the constant worked in the first example,
if we want to look at convergence issues later, link to #1649

@coveralls

Coverage Status

Coverage increased (+0.07%) when pulling edb986f on josef-pkt:fit_constrained into e4326b7 on statsmodels:master.

@coveralls

Coverage Status

Coverage increased (+0.09%) when pulling 7c901b9 on josef-pkt:fit_constrained into e4326b7 on statsmodels:master.

@josef-pkt
Member

TravisCI failure: I attached the auxiliary, transformed results to the regular results instance.
pickle_shrink test fails on all python.

I didn't attach it in Poisson.fit_regularized

@josef-pkt
Member

No, unrelated to fit_regularized (unit tests don't check that case)
adding llnull triggers the calculation of null which is cached.
However, why wasn't null already called before adding llnull?

@coveralls

Coverage Status

Coverage increased (+0.11%) when pulling d332246 on josef-pkt:fit_constrained into e4326b7 on statsmodels:master.

@josef-pkt josef-pkt commented on an outdated diff Jun 3, 2014
statsmodels/discrete/tests/test_constrained.py
+
+ @classmethod
+ def setup_class(cls):
+ from statsmodels.genmod.generalized_linear_model import GLM
+ from statsmodels.genmod import families
+ from statsmodels.base._constraints import fit_constrained
+
+ cls.res2 = results.results_exposure_constraint
+ cls.idx = [6, 2, 3, 4, 5, 0] # 2 is dropped baseline for categorical
+
+ # example with offset
+ formula = 'deaths ~ smokes + C(agecat)'
+ mod = GLM.from_formula(formula, data=data,
+ family=families.Poisson(),
+ offset=np.log(data['pyears'].values))
+ mod._init_keys.append('family')
@josef-pkt
josef-pkt Jun 3, 2014 Member

don't keep monkeys in the unit tests :)

@josef-pkt
Member

It took me a while to find out why my changes don't work anymore for the GLM-logit test case.

@coveralls

Coverage Status

Coverage increased (+0.11%) when pulling 026f60d on josef-pkt:fit_constrained into e4326b7 on statsmodels:master.

@josef-pkt
Member

The last commit adds test for GLM binary (logit) , main problem for tests is that Stata uses nans for fixed parameter bse and following, we use bse=0 and following

I added a mask to ignore those differences

                 Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:                      y   No. Observations:                   32
Model:                            GLM   Df Residuals:                       29
Model Family:                Binomial   Df Model:                            2
Link Function:                  logit   Scale:                             1.0
Method:                          IRLS   Log-Likelihood:                -12.890
Date:                Tue, 03 Jun 2014   Deviance:                       25.780
Time:                        23:11:10   Pearson chi2:                     27.1
No. Iterations:                     1                                         
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
x1             2.8000          0        inf      0.000         2.800     2.800
x2             0.0958      0.138      0.693      0.488        -0.175     0.367
x3             2.3717      1.007      2.355      0.019         0.398     4.346
const        -12.9465      3.340     -3.876      0.000       -19.494    -6.399
==============================================================================
. glm grade  gpa tuce psi, constraints(1)

Iteration 0:   log likelihood = -47.088878  

Generalized linear models                          No. of obs      =        32
Optimization     : ML                              Residual df     =        29
                                                   Scale parameter =   1.22583
Deviance         =  35.54906047                    (1/df) Deviance =   1.22583
Pearson          =  35.54906047                    (1/df) Pearson  =   1.22583

Variance function: V(u) = 1                        [Gaussian]
Link function    : g(u) = u                        [Identity]

                                                   AIC             =  3.130555
Log likelihood   = -47.08887826                    BIC             = -64.95728

 ( 1)  [grade]gpa = 2.8
------------------------------------------------------------------------------
             |                 OIM
       grade |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         gpa |        2.8  (constrained)
        tuce |   -.097778   .0512957    -1.91   0.057    -.1983158    .0027598
         psi |   .3872348   .3970722     0.98   0.329    -.3910125    1.165482
       _cons |  -6.408786   1.136086    -5.64   0.000    -8.635474   -4.182097
------------------------------------------------------------------------------
@josef-pkt
Member

There are more test results from Stata that haven't been use so far, unconstraint and constraint with vce(robust) option

I don't know what vce(robust) does in Stata, and I haven't figured out yet how to add sandwich robust covariances to GLM. Stata also allows HAC and cluster robust vce in GLM.

@josef-pkt josef-pkt closed this Jun 4, 2014
@josef-pkt josef-pkt reopened this Jun 4, 2014
@josef-pkt
Member

How can they put a close button where there is supposed to be a cancel comment button.

@coveralls

Coverage Status

Coverage increased (+0.16%) when pulling 502e959 on josef-pkt:fit_constrained into 8f7d124 on statsmodels:master.

@coveralls

Coverage Status

Coverage increased (+0.16%) when pulling 7db959e on josef-pkt:fit_constrained into 8f7d124 on statsmodels:master.

@josef-pkt
Member

@jseabold @kshedden I'm about done here, or as far as I want to go in this PR

only Poisson and GLM for now
included prototype for "pre_fix" function base._constraints.fit_constrained_wrap
still some open issues, but the two new methods look stable enough, I think

If the _wrap functions that modifies an existing model work out, then I think we can add more of those that without needing to incorporate them in every model. Pinv/PCR for non-linear models (Poisson and GLM again) would be my next transform candidate.

@coveralls

Coverage Status

Coverage increased (+0.17%) when pulling 7db959e on josef-pkt:fit_constrained into 8f7d124 on statsmodels:master.

@coveralls

Coverage Status

Coverage increased (+0.18%) when pulling dd66f64 on josef-pkt:fit_constrained into 8f7d124 on statsmodels:master.

@josef-pkt
Member

Ready for the green button, when TravisCI comes back green
(only one merge behind master)

Follow-up will be extension to other models, and getting a pattern and naming for Transformation classes and better support for "prefix" functions.

@josef-pkt josef-pkt added this to the 0.6 milestone Jun 5, 2014
@josef-pkt
Member

missing: adding to release notes

@josef-pkt josef-pkt merged commit 2c3eb6d into statsmodels:master Jun 5, 2014

1 check passed

continuous-integration/travis-ci The Travis CI build passed
Details
@josef-pkt josef-pkt deleted the josef-pkt:fit_constrained branch Jun 5, 2014
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment