Statsmodels: statistical modeling and econometrics in Python

@josef-pkt I am starting a PR . At the moment there is the gam file that contains the gam penalty class. Smooth_basis contains some functions to get bsplines and polynomial basis. This file will be removed when we will be able to use directly patsy.

There are also 2 files that contains examples or small scripts. We will remove them later.

Let me know what do you think about that.

Todo

• [ ] predict errors (stateful transform, patsy ?), Note fittedvalues are available only a problem in pirls example, predict after fit works, requires spline basis values
• [ ] get_prediction also errors, maybe consequence of predict error currently errors because weights is None
• [ ] check test coverage for offset and exposure

Interface

• [ ] partial_values plot_partial has inconvenient arguments (smoother and mask) instead of column index or term name or similar
• [ ] formula-like interface for predict (create spline basis values internally)
• [ ] adjust inherited methods like plot_partial_residuals (after GSOC?)
• [ ] default param names for splines (should be more informative than "xi" (i is range(k)), but less (?) verbose than patsy's), need them for test_terms

Here's a simple branch with the code added into statsmodels.tsa.statespace. A couple of thoughts

• At least in the dev process, I thought it might be nicer to keep it in its own module, rather than putting it with the kalmanf. I don't know what makes the most sense in the long run.
• I have unit tests that rely on the statespace model, but I'm rewriting them so the KF pull request can be done on its own without other dependencies, especially since the statespace model is likely to change.

Edit: Original line comments.

L72

Question: do we want to keep the single-precision version (and the complex single precision, below)? I don't really see a use case, and it appears from preliminary tests that the results tend to overflow. Maybe I'll post a unit test to demonstrate and we can go from there.

L332

Question: there are a bunch of ways to initialize the KF, depending on the theory underlying the data. This one is only valid for stationary processes. Probably best to move out to the Python wrapper?

L444 Question: I think we'll want to add the ability to specify whether or not to check for convergence and alter the tolerance.

L414

This inversion is using an LU decomposition, but I think in this case I can rely on f to be positive definite since it's the covariance matrix of the forecast error, in which case I could use the much faster Cholesky decomposition approach. This is something I'm looking into, but if you happen to know one way or the other, that would be great too.

Related to this is the idea that you should "never take inverses", and I guess I need to look into replacing this with a linear solver routine if possible, in the updating stage below.

This is a new version of the kernel density estimation. The purpose is to provide an implementation that is faster in the case of grid evaluation, and also works on bounded domains.

There is still some work to do, in particular, I need to add tests for the multi-dimensional and discrete densities.

This model include following models:

1. Generic Zero Inflated model
2. Zero-Inflated Poisson model
3. Zero-Inflated Generalized Poisson model (ZIGP-P)
4. Zero-Inflated Generalized Negative Binomial model (ZINB-P)

Each model include this parts:

• [x] LLF
• [x] Score
• [x] Hessian
• [x] Predict
• [x] Fit
• [x] Docs
• [x] Tests

Status: - reviewing, need to implement better way to generate start_params Last commit for end of GSoC17: Changed way to find start params

This PR include implementation of Generalized Poisson model This model include this parts:

• [x] Log-likelihood function
• [x] Score
• [x] Hessian
• [x] Fit
• [x] Result
• [x] Tests
• [x] Docs

Status - merged #3795

Hi, I have started implementing Kim Filter, outlined a basic functionality, as described in Kim-Nelson book (see diagram on p. 105). I didn't even run the code yet to check for errors. Coding style and class interface bother me more for the moment, as well as the possible ways to test it without implementing models.

xref #3773

I apologize for this taking longer than had hoped, but my busy period was exacerbated by some unforeseen events. But I'm somewhat back now.

This is the first step to solving the above issue.

Still a lot is missing... I've only really gotten loglike to work for all the families except Gamma, Binomial, Gaussian, and Tweedie (which has never had loglike). I also need to re-work the docstrings. I'm thinking its more logical to have a loglike_obs function called in each family and then have a loglike function in the Family class so that loglike will simply be inherited. I think you the doctrings could be elegantly written to handle this too.

I'm thinking I will work on deviance next and then ciricle back to loglike... I feel like R takes some computational shortcuts...

rebased version of #4575 which was rebased version of #2744 The original GSOC PR with most of the development discussion is #2435

rebase conflict in compat.python: This has now unneeded itertools.combinations import but I dropped the py 2.6 compat code.

Hello,

xref #3371

This (should) get var_weights to work for GLM. I added a Tweedie usecase against R.

I can't seem to get the Poisson with var_weights to HC0 test that was added (but disabled because of the lack of functionality) to work. The test is here:

https://github.com/thequackdaddy/statsmodels/blob/var_weight/statsmodels/genmod/tests/test_glm_weights.py#L142

It fails on the following assert

assert_allclose(res1.bse, corr_fact * res2.bse, atol= 1e-6, rtol=2e-6)

Its relatively close... consistently off by a factor of 0.98574... corr_fact brings it much closer... I'm wondering if another adjustment is necessary?

Honestly, I don't really understand much about sandwiches (except PBJ, meatball, and bánh mì).

Rest of the test seem to work pretty well.

I'd be happy to have this merged relatively soon, so thanks for the feedback and review!

PR as discussed in #3274. Tested with a SAS example and two R examples, produced the same results.

To-do:

• [X] Core stats computation
• [x] api
• [x] automatic create dummy variable and hypothesis testing for categorical type independent variables.
• [x] Input validation
• [x] More examples to be tested

references: [1] https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_introreg_sect012.htm) [2] GLM Algorithms ftp://public.dhe.ibm.com/software/analytics/spss/documentation/statistics/20.0/en/client/Manuals/IBM_SPSS_Statistics_Algorithms.pdf [3] http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.278.6976&rep=rep1&type=pdf

Unit test in test_MultivariateOLS.py

compare_r_output_dogs_data()

It reproduces results from the following R code:

library(car)
Drug = c('Morphine', 'Morphine', 'Morphine', 'Morphine', 'Morphine', 'placebo', 'placebo', 'placebo', 'placebo', 'placebo', 'Trimethaphan', 'Trimethaphan', 'Trimethaphan', 'Trimethaphan', 'Trimethaphan')
Depleted = c('N', 'N', 'N', 'N', 'Y', 'Y', 'Y', 'N', 'N', 'N', 'N', 'Y', 'Y', 'Y', 'Y')
subject  = c(1, 2, 3, 4, 5, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16)
Histamine0 = c(-3.218876, -3.912023, -2.65926, -1.771957, -2.302585, -2.65926, -2.995732, -3.506558, -3.506558, -2.65926, -2.407946, -2.302585, -2.525729, -2.040221, -2.813411)
Histamine1 = c(-1.609438, -2.813411, 0.336472, -0.562119, -2.407946, -2.65926, -2.65926, -0.478036, 0.04879, -0.18633, 1.141033, -2.407946, -2.407946, -2.302585, -2.995732)
Histamine3 = c(-2.302585, -3.912023, -0.733969, -1.049822, -2.040221, -2.813411, -2.813411, -1.171183, -0.314711, 0.067659, 0.722706, -2.407946, -2.407946, -2.120264, -2.995732)
Histamine5 = c(-2.525729, -3.912023, -1.427116, -1.427116, -1.966113, -2.65926, -2.65926, -1.514128, -0.510826, -0.223144, 0.207014, -2.525729, -2.302585, -2.120264, -2.995732)
data = data.frame(Histamine0, Histamine1, Histamine3, Histamine5)
hismat = as.matrix(data[,1:4])
result = lm(hismat ~ Drug * Depleted)
linearHypothesis(result, c(1, 0, 0, 0, 0, 0)) 
linearHypothesis(result, t(cbind(c(0, 1, 0, 0, 0, 0), c(0, 0, 1, 0, 0, 0))))
linearHypothesis(result, c(0, 0, 0, 1, 0, 0)) 
linearHypothesis(result, t(cbind(c(0, 0, 0, 0, 1, 0), c(0, 0, 0, 0, 0, 1))))
# Or ManRes <- Manova(result, type="III")

test_affine_hypothesis()

It reproduces results from the following R code:

result = lm(hismat ~ Drug*Depleted)
fml = t(cbind(c(0, 1.1, 1.2, 1.3, 1.4, 1.5), c(0, 2.1, 3.2, 3.3, 4.4, 5.5)))
linearHypothesis(result, fml, t(cbind(c(1, 2, 3, 4), c(5, 6, 7, 8))), verbose=TRUE)

Ok, I wanted to make a PR for this, there are still a couple of things that I need to fix but things are pretty close to done and I'm happy with the current approach.

The key part of the current approach is a function distributed_estimation. This function works by taking a generator for endog and exog, endog_generator and exog_generator, as well as a series of functions and key word arguments to be run on each machine and then used to recombine the results. The generator approach allows for a variety of use cases and can handle a lot of the ideas discussed in the initial proposal. For each data set yielded by the generators a model is initialized using model_class and init_kwds and then the function estimation_method is applied to the model along with the key words fit_kwds and estimation_kwds. Finally, the results are recombined from each data set using join_method.

Currently, this defaults to the distributed regularized approach discussed here:

http://arxiv.org/pdf/1503.04337v3.pdf

but the way I've set things up means that the user should be able to apply any number of procedures here.

The current todo list:

• [x] Fix hess_obs
• [x] Fix joblib fit
• [x] Add dask fit
• [x] Add WLS/GLS for debiased regularized fit
• [x] Add likelihood result
• [x] Add data tests

Let me know if there are any comments/questions/criticisms, as I've said, it certainly isn't complete but I wanted to get this out there so I can start integrating any changes as I finish it up.

GEE does not have use_t option, use_t = False is hardcoded in summary, no fit option

Monte Carlo comparison. they use both t and z/normal based tests and pvalues t-test can be conservative with some of the bias corrections and small number of clusters

Also we might add more options for small sample adjusted, bias reduced cov_types related to variants of cluster correlation robust sandwiches, CR3, ...

Li, Peng, and David T. Redden. “Small Sample Performance of Bias-Corrected Sandwich Estimators for Cluster-Randomized Trials with Binary Outcomes.” Statistics in Medicine 34, no. 2 (2015): 281–96. https://doi.org/10.1002/sim.6344

reference found in Huang, Shuang, Mallorie H Fiero, and Melanie L Bell. “Generalized Estimating Equations in Cluster Randomized Trials with a Small Number of Clusters: Review of Practice and Simulation Study.” Clinical Trials 13, no. 4 (August 1, 2016): 445–49. https://doi.org/10.1177/1740774516643498.

also compares several small sample corrections

this article has pseudo-hatmatrix for small sample corrections for HC2, HC3, ...

moment condition z (y - g(x, beta)) i.e. nonlinear IV, but not for generic nonlinear GMM

Lin, Eric S., and Ta-Sheng Chou. “Finite-Sample Refinement of GMM Approach to Nonlinear Models under Heteroskedasticity of Unknown Form.” Econometric Reviews 37, no. 1 (January 2, 2018): 1–28. https://doi.org/10.1080/07474938.2014.999499.

I did not find any HC3 equivalent for generic GMM, and nothing for more general non-gaussian MLE. (Poisson, GLM use the IRLS, WLS hatmatrix)

parking a reference for a useful approach to cluster bootstrap, "resampling" clusters

Instead of actually resampling clusters which would result in different nobs in replication samples, we can use weights by cluster.

Without computational shortcut, we can just use freq_weights in GLM or cluster weights in GEE and compute full model for each replication sample. I don't know how much we gain in computational efficiency if we just adjust weights in a single model instance. We might only safe on model __init__ setup cost.

Cheng, Guang, Zhuqing Yu, and Jianhua Z. Huang. “The Cluster Bootstrap Consistency in Generalized Estimating Equations.” Journal of Multivariate Analysis 115 (March 1, 2013): 33–47. https://doi.org/10.1016/j.jmva.2012.09.003.

parking references

we already have several issues for variation on GEE/GLM including variance modeling #2898, #5674

Here the extension is to use "standard" GEE but with underlying score functions that are not coming from a LEF family, e.g. zero-inflated count models.

How do we reuse the pattern for GEE estimation with working correlation for score functions or moment conditions that do not come from GLM (mean, variance functions)? I don't know what the theoretical properties are. Non-LEF, non-GLM do not have the consistency properties robust to misspecification.

Hall, Daniel B, and Zhengang Zhang. “Marginal Models for Zero Inflated Clustered Data.” Statistical Modelling 4, no. 3 (October 1, 2004): 161–80. https://doi.org/10.1191/1471082X04st076oa.

Kong, Maiying, Sheng Xu, Steven M. Levy, and Somnath Datta. “GEE Type Inference for Clustered Zero-Inflated Negative Binomial Regression with Application to Dental Caries.” Computational Statistics & Data Analysis 85 (May 1, 2015): 54–66. https://doi.org/10.1016/j.csda.2014.11.014.

Is your feature request related to a problem? Please describe

The statsmodels.tsa.seasonal.seasonal_decompose method throws errors for not specifying a frequency and having NaN values, which can be conveniently addressed within the method by inferring frequencies and dropna() respectively.

Describe the solution you'd like

Quality of life improvements: include the relevant pandas transformations within the seasonal_decompose method. Include detailed error messages when DateTimeIndex is broken, etc.

Describe alternatives you have considered

NIL

Additional context

I can work on this, but I would like to understand reasons behind the current implementation first.