# Pushing Julia/statistics development

**URL:** <https://discourse.julialang.org/t/pushing-julia-statistics-development/80111>\
**Category:** Statistics\
**Created:** [April 27, 2022, 8:23am UTC](https://discourse.julialang.org/t/pushing-julia-statistics-development/80111 "2022-04-27T08:23:22Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ajay\_Shah](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ajay_shah/32/10520_2.png) [@Ajay\_Shah](https://discourse.julialang.org/u/Ajay_Shah)\
**Post date:** [April 27, 2022, 8:23am UTC](https://discourse.julialang.org/t/pushing-julia-statistics-development/80111/1 "2022-04-27T08:23:22Z")

</div>

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](https://github.com/JuliaStats/GLM.jl/pull/480)) and one new link function ([Power](https://github.com/JuliaStats/GLM.jl/pull/474)) and some small improvements ([PR](https://github.com/JuliaStats/GLM.jl/pull/481), [PR](https://github.com/JuliaStats/GLM.jl/pull/482)). **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](https://github.com/xKDR/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](https://github.com/xKDR/TSx.jl)”): 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](https://github.com/xKDR/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](https://github.com/xKDR/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](https://github.com/xKDR/DtD.jl) to measure the credit risk to a firm.

7. Small improvements: We have a good [Lowess.jl](https://github.com/xKDR/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](https://www.itl.nist.gov/div898/strd/general/bkground.html), 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](https://github.com/xKDR) . The persons in this project are:

From [XKDR Forum](https://xkdr.org/): [Susan Thomas](https://scholar.google.com/citations?user=0QRD90YAAAAJ&hl=en), [Ayush Patnaik](https://github.com/ayushpatnaikgit), [Ajay Shah](https://www.mayin.org/ajayshah/).

Independent researchers: [Mousum Dutta](https://www.linkedin.com/in/mousumdutta/), [Chirag Anand](https://github.com/chiraganand).

From [Chennai Mathematical Institute](https://www.cmi.ac.in/): [Sourish Das](https://www.cmi.ac.in/~sourish/), undergraduate students ([Siddhant Chaudhary](https://github.com/codetalker7), [Harsh Arora](https://github.com/harsharora21), [Naman Kumar](https://github.com/matcauthon49)), masters students ([Arnab Sen](https://github.com/ArnabbLank), [Anisha Saha](https://github.com/anisha0325), [Tanuj Sur](https://github.com/surtantheta), [Sumeet Suley](https://github.com/sumeetsuley), [Gudeet Siyan](https://github.com/16bitmood)).

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](https://julialang.org/jsoc/gsoc/juliastats/).

---

<div class="post-metadata">

**Author:** ![nalimilan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nalimilan/32/147_2.png) [@nalimilan](https://discourse.julialang.org/u/nalimilan)\
**Post date:** [April 27, 2022, 10:15am UTC](https://discourse.julialang.org/t/pushing-julia-statistics-development/80111/2 "2022-04-27T10:15:04Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![palday](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/palday/32/12640_2.png) [@palday](https://discourse.julialang.org/u/palday)\
**Post date:** [April 27, 2022, 4:39pm UTC](https://discourse.julialang.org/t/pushing-julia-statistics-development/80111/3 "2022-04-27T16:39:41Z")

</div>

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](https://beacon.bio). In addition to our contributions to the JuliaStats organization, we’ve also open-sourced several statistical packages developed in our professional roles:

- [Effects.jl](https://github.com/beacon-biosignals/Effects.jl) Similar in spirit to R’s `effects` package, but quite different (more Julian) in interface
- [StandardizedPredictors.jl](https://github.com/beacon-biosignals/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. 🙂

---

<div class="post-metadata">

**Author:** ![nalimilan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nalimilan/32/147_2.png) [@nalimilan](https://discourse.julialang.org/u/nalimilan)\
**Post date:** [April 27, 2022, 6:26pm UTC](https://discourse.julialang.org/t/pushing-julia-statistics-development/80111/4 "2022-04-27T18:26:10Z")

</div>

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

> [@Ajay\_Shah](#):
>
> Improvements to GLM: We have done one new distribution ([Geometric](https://github.com/JuliaStats/GLM.jl/pull/480)) and one new link function ([Power](https://github.com/JuliaStats/GLM.jl/pull/474)) and some small improvements ([PR](https://github.com/JuliaStats/GLM.jl/pull/481), [PR](https://github.com/JuliaStats/GLM.jl/pull/482)). **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](https://github.com/JuliaStats/GLM.jl/issues/426) 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 [https://github.com/JuliaStats/GLM.jl/pull/340](https://github.com/JuliaStats/GLM.jl/pull/340)? It’s almost ready and could easily be finished and merged.

> [@Ajay\_Shah](#):
>
> A framework for statistical models ([CRRao.jl](https://github.com/xKDR/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](https://github.com/xKDR/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](https://github.com/xKDR/TSx.jl)”): 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.

---

<div class="post-metadata">

**Author:** ![chiraganand](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chiraganand/32/32787_2.png) [@chiraganand](https://discourse.julialang.org/u/chiraganand)\
**Post date:** [April 29, 2022, 6:59am UTC](https://discourse.julialang.org/t/pushing-julia-statistics-development/80111/5 "2022-04-29T06:59:42Z")

</div>

> 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.

---

<div class="post-metadata">

**Author:** ![mousum](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mousum/32/35843_2.png) [@mousum](https://discourse.julialang.org/u/mousum)\
**Post date:** [April 30, 2022, 4:37am UTC](https://discourse.julialang.org/t/pushing-julia-statistics-development/80111/6 "2022-04-30T04:37:40Z")

</div>

> [@nalimilan](#):
>
> As @palday noted GLM already supports the QR decomposition but it’s not the default. We discussed this [in the past](https://github.com/JuliaStats/GLM.jl/issues/426) 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](https://github.com/JuliaStats/GLM.jl/pull/340) to add the multicollinearity feature in the existing `GLM`
2. A fullfledged QR decomposition in `GLM`

---

<div class="post-metadata">

**Author:** ![alonsoC1s](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alonsoc1s/32/25819_2.png) [@alonsoC1s](https://discourse.julialang.org/u/alonsoC1s)\
**Post date:** [May 5, 2022, 12:51pm UTC](https://discourse.julialang.org/t/pushing-julia-statistics-development/80111/7 "2022-05-05T12:51:11Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![ayushpatnaikgit](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ayushpatnaikgit/32/25229_2.png) [@ayushpatnaikgit](https://discourse.julialang.org/u/ayushpatnaikgit)\
**Post date:** [May 5, 2022, 6:34pm UTC](https://discourse.julialang.org/t/pushing-julia-statistics-development/80111/8 "2022-05-05T18:34:41Z")

</div>

Hi @alonsoC1s, email me at [ayushpatnaik@gmail.com](mailto:ayushpatnaik@gmail.com)  
Then we’ll see 🙂

---

<div class="post-metadata">

**Author:** ![ParadaCarleton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paradacarleton/32/20005_2.png) [@ParadaCarleton](https://discourse.julialang.org/u/ParadaCarleton)\
**Post date:** [May 9, 2022, 1:47am UTC](https://discourse.julialang.org/t/pushing-julia-statistics-development/80111/9 "2022-05-09T01:47:43Z")

</div>

Can I ask – is there a reason why [HypothesisTests.jl’s Anderson-Darling](https://juliastats.org/HypothesisTests.jl/stable/nonparametric/#Anderson-Darling-test) 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.

---

<div class="post-metadata">

**Author:** ![BioTurboNick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bioturbonick/32/6380_2.png) [@BioTurboNick](https://discourse.julialang.org/u/BioTurboNick)\
**Post date:** [May 25, 2022, 5:31pm UTC](https://discourse.julialang.org/t/pushing-julia-statistics-development/80111/10 "2022-05-25T17:31:14Z")

</div>

FWIW, I think AD has to have custom probability distributions computed for each non-Normal distribution you want to test. So it can in theory test many distribtuions, but only a subset are supported.

---

<div class="post-metadata">

**Author:** ![ParadaCarleton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paradacarleton/32/20005_2.png) [@ParadaCarleton](https://discourse.julialang.org/u/ParadaCarleton)\
**Post date:** [May 28, 2022, 8:52pm UTC](https://discourse.julialang.org/t/pushing-julia-statistics-development/80111/11 "2022-05-28T20:52:31Z")

</div>

No need for that, the Anderson-Darling statistic is independent of the distribution you’re testing.

---

<div class="post-metadata">

**Author:** ![BioTurboNick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bioturbonick/32/6380_2.png) [@BioTurboNick](https://discourse.julialang.org/u/BioTurboNick)\
**Post date:** [May 28, 2022, 10:17pm UTC](https://discourse.julialang.org/t/pushing-julia-statistics-development/80111/12 "2022-05-28T22:17:27Z")

</div>

Do you have a source for that?

[https://www.itl.nist.gov/div898/handbook/eda/section3/eda35e.htm](https://www.itl.nist.gov/div898/handbook/eda/section3/eda35e.htm)

> The Anderson-Darling test makes use of the specific distribution in calculating critical values. This has the advantage of allowing a more sensitive test and the disadvantage that critical values must be calculated for each distribution.

---

<div class="post-metadata">

**Author:** ![ParadaCarleton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paradacarleton/32/20005_2.png) [@ParadaCarleton](https://discourse.julialang.org/u/ParadaCarleton)\
**Post date:** [May 28, 2022, 10:30pm UTC](https://discourse.julialang.org/t/pushing-julia-statistics-development/80111/13 "2022-05-28T22:30:32Z")

</div>

Technically it’s only asymptotically that the dependency vanishes:  
[https://www.jstor.org/stable/pdf/2237910.pdf](https://www.jstor.org/stable/pdf/2237910.pdf)

---

<div class="post-metadata">

**Author:** ![BioTurboNick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bioturbonick/32/6380_2.png) [@BioTurboNick](https://discourse.julialang.org/u/BioTurboNick)\
**Post date:** [May 29, 2022, 3:36am UTC](https://discourse.julialang.org/t/pushing-julia-statistics-development/80111/14 "2022-05-29T03:36:44Z")

</div>

After some poking around, I think it’s just discrete univariate distributions that really fail.

```julia
julia> OneSampleADTest.([rand(Poisson(2), 100) for i ∈ 1:10], Poisson(2)) .|> pvalue
10-element Vector{Float64}:
 0.001966853989250339
 6.3489859959409145e-6
 0.0005871355756167373
 3.645708576505147e-5
 6.000000011940898e-6
 6.000241096759673e-6
 0.00012823416159524204
 0.0005448765568399905
 1.0625078049386616e-5
 4.366874769401452e-5

```

---

<div class="post-metadata">

**Author:** ![Sourish](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sourish/32/34156_2.png) [@Sourish](https://discourse.julialang.org/u/Sourish)\
**Post date:** [August 8, 2022, 6:52pm UTC](https://discourse.julialang.org/t/pushing-julia-statistics-development/80111/15 "2022-08-08T18:52:09Z")

</div>

Shapiro-Wilk test was developed to test particularly for testing Normality. If you have a bunch of samples, say X1, X2,…,Xn ~ f(x), and you want to test H0: f follows Normal vs Ha: f does not follow Normal.

On the other hand, the purpose of the Anderson-Darling test is completely different. The purpose of the AD test is to check if you have two different samples. If both the samples are generated from the same distributions or not. That is if you have one sample say {X1,X2,…Xn} ~ f(x) and you have another sample {Y1,Y2,…,Ym} ~ g(y) you would like to test:

H0: f(x) = g(y) vs Ha: f(x) =!= g(y)

It does not require any distributional assumption. It is a non-parametric test. Now what you can do is something like this. Suppose you have only {X1, X2,…Xn} ~ f(x), and you want to check if  
H0: f(x) = N(0,1) vs Ha: f(x) =!= N(0,1)

You can simulate {y1,y2,…,ym} from N(0,1) and plug it in AD test. This is the practice that is typically followed. I believe the statistical power varies significantly depending on the distributions.