Pushing Julia/statistics development

[Ajay_Shah]

Julia is great, but there are many limitations in putting it to work for doing statistics. With the help of Viral, Alan, and @bkamins, a bunch of us got going on a systematic attack on this problem. We have gone after this top-down and bottom-up:

Top down: We form opinions about important gaps in the Julia package ecosystem and set about filling these.

Bottom up: We have begun on full implementations of applied statistics research papers using Julia, thus discovering limitations along the way. We file bug reports, tests and feature requests, and also try to build the requisite code either as PRs or as new packages.

We have done:

1. Improvements to GLM: We have done one new distribution (Geometric) and one new link function (Power) and some small improvements (PR, PR). Coming up: (1) A choice between the (existing) fast+imprecise Cholesky decomposition vs. (a new) slower+precise QR decomposition in the iterative least squares, (2) Improved handling of collinear data, and (3) A paper with comparisons between R/Julia/SAS GLM implementations on features, performance and correctness.

2. A framework for statistical models (CRRao.jl): Applied statisticians value a consistent API for a wide variety of statistical models. The package embeds a consistent API, and a group of models that are ready to use. Coming up: We will write more models, we hope others will also build new models in this framework. We will build out the API to support a sensible workflow for statistical modelling.

3. A time series class (“TSx”): For working with time series data, a set of metaphors and operations are required. Our TSx package is syntactic sugar on top of the powerful capabilities of DataFrames.jl and is easily maintainable. In previous years, we have been intensive users of zoo and xts in R, for working with financial and macro data, and we have brought these experiences to bear on the design of TSx while keeping the design flexible enough to incorporate use cases from other fields.

4. Working with survey data (Survey.jl): We have a small set of much-used functions for working with survey data.

5. Working with the VIIRS night lights data (NighttimeLights.jl): Satellite imagery of nighttime lights is a valuable path to observing economic prosperity at high-frequency and high resolution. This package is a complete set of steps for cleaning and bias correcting the raw data that’s released by NASA/NOAA.

6. Distance-to-Default: We have implemented the Merton Model in DtD.jl to measure the credit risk to a firm.

7. Small improvements: We have a good Lowess.jl.

8. Coming up: Measuring the precision of simple statistical calculations (NISTTests.jl): We are packaging a group of test cases from the US NIST, as functions that measure the precision of a supplied Julia function.

We are keen to make these good; please do criticise our work and help us make it better.

Our main page is https://github.com/xKDR . The persons in this project are:

From XKDR Forum: Susan Thomas, Ayush Patnaik, Ajay Shah.

Independent researchers: Mousum Dutta, Chirag Anand.

From Chennai Mathematical Institute: Sourish Das, undergraduate students (Siddhant Chaudhary, Harsh Arora, Naman Kumar), masters students (Arnab Sen, Anisha Saha, Tanuj Sur, Sumeet Suley, Gudeet Siyan).

We welcome your interest and involvement in carrying this work forward. Folks in India: the team is located in Bombay, Delhi, Madras and Pune: we can readily meet up in any of these places. We are part of the 2022 Google Summer of Code.

[nalimilan]

Thanks for all of this! I was starting to think a larger project was underlying all these cool PRs.

I’m happy to discuss any lacks you may identify in the existing API in StatsAPI, DataAPI or other packages so that we can develop a consistent interface across the ecosystem.

[palday]

This is great to hear!

Re GLM.jl: The QR decomposition is already implemented, it’s just not exposed in a convenient user-facing way (e.g. via a keyword argument to glm()/lm()).

Several of the JuliaStats contributors (myself, @ararslan, @dave.f.kleinschmidt) have day jobs at Beacon Biosignals. In addition to our contributions to the JuliaStats organization, we’ve also open-sourced several statistical packages developed in our professional roles:

• Effects.jl Similar in spirit to R’s effects package, but quite different (more Julian) in interface
• StandardizedPredictors.jl Do scaling, centering, other standardization as part of your model formula. Provides a nice display of relevant bits of info in your coef table (e.g. center point) and plays nicely with things like predict.

We also have a few more contributions planned in the foreseeable future.

[nalimilan]

A few more specific comments on some of the topics you raised:

@Ajay_Shah

Improvements to GLM: We have done one new distribution (Geometric ) and one new link function (Power ) and some small improvements (PR , PR ). Coming up : (1) A choice between the (existing) fast+imprecise Cholesky decomposition vs. (a new) slower+precise QR decomposition in the iterative least squares, (2) Improved handling of collinear data, and (3) A paper with comparisons between R/Julia/SAS GLM implementations on features, performance and correctness.

As @palday noted GLM already supports the QR decomposition but it’s not the default. We discussed this in the past and I agree it would make sense to use it by default as Cholesky seems to be less robust in ill-conditioned cases (including nearly-collinear variables) where R succeeds.

Regarding collinear data, are you aware of Add support for rank deficient `GeneralizedLinearModel` by jason-xuan · Pull Request #340 · JuliaStats/GLM.jl · GitHub? It’s almost ready and could easily be finished and merged.

@Ajay_Shah

A framework for statistical models (CRRao.jl ): Applied statisticians value a consistent API for a wide variety of statistical models. The package embeds a consistent API, and a group of models that are ready to use. Coming up : We will write more models, we hope others will also build new models in this framework. We will build out the API to support a sensible workflow for statistical modelling.

I don’t completely understand the goal of this package. Thanks to StatsModels (and StatsAPI) modeling packages are already supposed to use an unified API. This can certainly be improved, but I would think this requires adding new elements to the common API (keyword argument names…) and ensuring that each package implements it, rather than having One Package to Rule Them All. Is there something I’m missing? (I must say I don’t really get the principle of Zelig either. )

@Ajay_Shah

Working with survey data (Survey.jl ): We have a small set of much-used functions for working with survey data.

This is an area that is particularly useful for a social scientist like me, so a big +1!

@Ajay_Shah

A time series class (“TSx ”): For working with time series data, a set of metaphors and operations are required. Our TSx package is syntactic sugar on top of the powerful capabilities of DataFrames.jl and is easily maintainable. In previous years, we have been intensive users of zoo and xts in R, for working with financial and macro data, and we have brought these experiences to bear on the design of TSx while keeping the design flexible enough to incorporate use cases from other fields.

Nice. Wrapping a DataFrame is indeed an interesting alternative to TimeArray when you need to store heterogeneous columns. It would be good to try using consistent interfaces across TSx, TimeSeries and DataFrames where it makes sense.

[chiraganand]

Nice. Wrapping a DataFrame is indeed an interesting alternative to TimeArray when you need to store heterogeneous columns. It would be good to try using consistent interfaces across TSx, TimeSeries and DataFrames where it makes sense.

Yes, it would be good to have consistent interfaces across TSx and TimeSeries, though, we are not considering DataFrames because, IMHO, DataFrames syntax isn’t the best for timeseries-like data manipulation. But, if users do want to use DataFrames syntax then they can always use TS.coredata property to directly manipulate the underlying DataFrame albeit with certain precautions.

As of now, TSx interfaces somewhat resemble the zoo and xts syntax in R but I think we can create method aliases to provide users with an option of using TimeSeries interfaces.

[mousum]

@nalimilan

As @palday noted GLM already supports the QR decomposition but it’s not the default. We discussed this in the past and I agree it would make sense to use it by default as Cholesky seems to be less robust in ill-conditioned cases (including nearly-collinear variables) where R succeeds.

Thanks a lot, @palday and @nalimilan.
Yes, we already have QR decomposition in place, we need to add a few more delbeta! functions to accommodate multicollinearity/rank deficiency in lm and glm models.
So our plan is something like the following:

1. Move forward with PR #340 to add the multicollinearity feature in the existing GLM
2. A fullfledged QR decomposition in GLM

[alonsoC1s]

I have personaly noticed a lot of features you would expect to use on a multivariate statistics course are missing. For instance a variety of normality tests and others that you can readily find in R.
I would love to be part of this effort to revamp Julia’s stats capabilities. Is there a way I or anyone else could join this concerted effort?

[ayushpatnaikgit]

Hi @alonsoC1s, email me at ayushpatnaik@gmail.com
Then we’ll see

[Lime]

Can I ask – is there a reason why HypothesisTests.jl’s Anderson-Darling doesn’t work for you? It’s probably the most common test for normal distributions, alongside Shapiro-Wilk and Kolmogorov-Smirnov. I believe the AD test was chosen over Shapiro-Wilk because SW does poorly against thick tails and datasets with many ties, and is also nongeneralizable – it only works for testing the normal distribution, unlike Anderson-Darling which can be used to test any model. KS has very bad power in general.