Add Pearsons correlation to _column_associations

[reshamas]

Towards #1176

Add Pearsons correlation to column associations.

[reshamas]

This is the output, which I excluded for now, due to errors:

       left_column_name  left_column_idx right_column_name  right_column_idx  cramer_v  pearson
    0               c_3                3             c_str                 5    0.8215      NaN
    1               c_1                1               c_4                 4    0.8215   0.1597
    2               c_0                0               c_1                 1    0.8215   0.1123
    3               c_2                2             c_str                 5    0.7551      NaN
    4               c_0                0             c_str                 5    0.7551      NaN
    5               c_0                0               c_3                 3    0.7551   0.3212
    6               c_1                1               c_3                 3    0.6837  -0.1887
    7               c_0                0               c_4                 4    0.6837  -0.3202
    8               c_4                4             c_str                 5    0.6837      NaN
    9               c_3                3               c_4                 4    0.6053  -0.0150
    10              c_2                2               c_3                 3    0.6053   0.1757
    11              c_1                1             c_str                 5    0.6053      NaN
    12              c_0                0               c_2                 2    0.6053  -0.0578
    13              c_2                2               c_4                 4    0.5169  -0.2885
    14              c_1                1               c_2                 2    0.4122  -0.4986

Not sure why Pearson is showing up as NaN for the diagonal correlations, they should be 1.000

       left right  pearson
    0   c_0   c_0   1.0000
    1   c_0   c_1   0.1123
    2   c_0   c_2  -0.0578
    3   c_0   c_3   0.3212
    4   c_0   c_4  -0.3202
    5   c_1   c_0   0.1123
    6   c_1   c_1   1.0000

[Vincent-Maladiere]

Hey @reshamas, thank you for starting this!

In skrub, we use methods defined in skrub._dataframe._common.py to dispatch operations on polars or pandas dataframes. We do this to avoid copying dataframes from one backend to the other, and later to enable features like lazy evaluation. Therefore, the input of the function skrub._column_associations.column_associations() could be either a pandas or a polars DataFrame.

So, here, we would like to use something like sbd.corr() ("sbd" stands for "skrub dataframe"), instead of pandas corr. Since this method doesn't exist yet, you'll have to add it to skrub._dataframe._common.py 😁

Alternatively, you could turn your dataframe into a numpy array with skb.to_numpy() and then use np.corrcoef. I slightly prefer the former suggestion because it adds corr() in skrub, which I think would be also useful later.

If you go with the first option, I suggest converting the output to a numpy 2d array, so that you can directly use the _stack_symmetric_associations() function to melt the symmetric matrix into a long format. You need to select only numeric columns when passing your dataframe to _stack_symmetric_associations().

Finally, you can apply the merge operation with skrub._join_utils.left_join.

We know that developing skrub is fairly involved at the moment, and we plan on making a simple documentation for all our internals. In the meantime, I'm here to help you to move this PR forward!

[Vincent-Maladiere]

Hey @reshamas, do you plan to continue working on this?

[reshamas]

@Vincent-Maladiere No, I won't be continuing work on this.

[Vincent-Maladiere]

@jeromedockes, I have an error on the TruncatedSVD for Windows; have you experienced this recently?

[jeromedockes]

I don't think so 🤔

[Vincent-Maladiere]

Ok, it's gone!

[jeromedockes]

works great 🎉

[jeromedockes]

Suggested change

	<td> {{ association["pearson_corr"] | format_number }} </td>
	<td>
	{%- if not (association['pearson_corr']) | is_null -%}
	{{ association["pearson_corr"] | format_number }}
	{%- endif -%}
	</td>

maybe we could do this, so that we don't print the "nan"s

[Vincent-Maladiere]

very interesting to see in your iris dataset example that these correlations are not concordant (third line: 0.395 vs 0.818, compared to the two lines above where the Cramer's is much higher)

[jeromedockes]

Suggested change

	<th scope="col">Pearson's Correlation</a></th>
	<th scope="col"><a href="https://en.wikipedia.org/wiki/Pearson_correlation_coefficient">Pearson's Correlation</a></th>

to be consistent with the other column we can add the link even though most people know what pearson correlation is. if you prefer not to add the link you need to remove the closing </a> tag, as there is no corresponding opening <a>

[Vincent-Maladiere]

woops, html is forgiving

[jeromedockes]

Suggested change

	pandas.
	pandas or polars, depending on the input dataframe.

[jeromedockes]

if we have the link to pandas doc would it make sense to ahve the polars one?

[Vincent-Maladiere]

Only if we use pl.DataFrame.corr() I'd say, otherwise we should mention what we use for polars inputs

[jeromedockes]

very interesting to see in your iris dataset example that these correlations are not concordant (third line: 0.395 vs 0.818, compared to the two lines above where the Cramer's is much higher)

yeah that's true... the cramer thing is quite rough and does aggressive binning and subsampling... it will be interesting to investigate a bit more

[jeromedockes]

very interesting to see in your iris dataset example that these correlations are not concordant (third line: 0.395 vs 0.818, compared to the two lines above where the Cramer's is much higher)

yeah that's true... the cramer thing is quite rough and does aggressive binning and subsampling... it will be interesting to investigate a bit more

it should be a bit more similar to spearman's I guess

[jeromedockes]

I guess polars and pandas have a different way of handling nans, pandas drops them whereas in polars they result in NaN correlations. WDYT would be the right behavior here? I think the pandas one because we already have other ways of indicating when there are nans so even if it should be interpreted cautiously the non-nan correlation is more informative. it will require a bit more work so could also be done in a separate pr

[jeromedockes]

and finally, have you checked the computation time for largish dataframes, does it change a lot the total time of the report?

[jeromedockes]

I think we could simplify this a bit

• select numeric columns only
• use sbd.pearson_corr on those
• use pandas stack().reset_index() to convert to tall format
• or in polars case insert a column like with_column(left_col=<the column names>) then use df.unpivot(index=left_col) (only we would use a random name instead of left_col to have no name clash)
• then do the left join as you are doing in the caller

we would no longer use stack_symmetric_associations for this one I guess, but we would no longer need to do the filling with np.nan or the loops over columns. WDYT?

[Vincent-Maladiere]

Yes, I hesitated with this approach, but I saw avoiding manually dispatching between polars and pandas as an advantage, and reusing stack_symmetric_associations also makes things easier to understand IMO. Note that you wouldn't need to select numeric because sbd.pearson_corr does it for you. I agree that adding NaNs columns is a bit convoluted, though. WDYT?

[jeromedockes]

ok I see what you mean. I still think a nice advantage of using a left join as you did is that we don't care if we have different columns in the 2 pairwise matrices, and it's not that bad to do the dispatch with the decorator we will just have 2 helpers like

@_stack.specialize('pandas')
def _stack_pandas(df, ...):
    # one-liner

(private to this module and with only the parameters we need here, not a general one in skrub._dataframe), and we don't need to insert the numpy nans and then convert all the columns to numpy and do the hstack

but both options make sense

[Vincent-Maladiere]

Ok, will do this option

[Vincent-Maladiere]

I guess polars and pandas have a different way of handling nans, pandas drops them whereas in polars they result in NaN correlations. WDYT would be the right behavior here? I think the pandas one because we already have other ways of indicating when there are nans so even if it should be interpreted cautiously the non-nan correlation is more informative. it will require a bit more work so could also be done in a separate pr

Right, I guess polars .drop_nulls().corr() do the trick.

[jeromedockes]

Right, I guess polars .drop_nulls().corr() do the trick.

I think it's a bit more complicated than that because you want to drop nulls for each pair of columns. when computing the correlation between columns A and B I don't want to drop rows that have a null in column C but have values for A and B:

>>> df = pd.DataFrame({'a': [None, 1, 2, 3], 'b': [0, None, 2, 3], 'c': [0, 1, None, 3]})
>>> df.corr()
     a    b    c
a  1.0  1.0  1.0
b  1.0  1.0  1.0
c  1.0  1.0  1.0
>>> df.dropna().corr()
    a   b   c
a NaN NaN NaN
b NaN NaN NaN
c NaN NaN NaN

so it might be easier to leave that for a different PR

[Vincent-Maladiere]

Yes right, I realized that just now. Not a huge deal to compute corr pairwise, though. Worse case scenario we can convert everything to pandas and forget about it

[Vincent-Maladiere]

Ok, I did some benchmarking with both solutions (pairwise or converting to pandas). On the figures below, "module" indicates the input dataframe module, i.e. the figures display the pandas corr() solution against one the two solutions for polars input (not polars corr()).
Some conclusions:

• The pairwise solution doesn't scale well with the number of columns.
• The "converting to pandas" solution follows closely the performances of pandas corr, and doesn't lead to an bottleneck, computing time-wise.
• Pandas corr scales well with the number of rows. At 1e6 rows and 100 columns, it takes ~4s so we might want to subsample at some point
• However, Pandas corr gets very slow when the number of columns > 1000, and takes up to 120s for 10k columns. When the number of columns is this high, we might want to avoid computing Pearson's correlation on it.
  WDYT?

Pairwise

@pearson_corr.specialize("polars", argument_type="DataFrame")
def _pearson_corr_polars(df):
    col_names = df.columns
    n_cols = len(col_names)

    out = np.eye(n_cols)

    # Compute pairwise correlations for upper triangle
    for i, col1 in enumerate(col_names):
        for j in range(i + 1, n_cols):
            col2 = col_names[j]
            out[i, j] = out[j, i] = df.select([col1, col2]).drop_nulls().corr()[0, 1]

    return pl.DataFrame(
        out, schema={col: pl.Float64 for col in col_names}, orient="col"
    )

Converting to pandas

@pearson_corr.specialize("polars", argument_type="DataFrame")
def _pearson_corr_polars(df):
    return pl.from_pandas(_pearson_corr_pandas(df.to_pandas()))

Row-wise performance:

[Vincent-Maladiere]

For completeness, here is the same benchmark for Cramer's V. For 1000 rows and 10k columns, it crashed my machine 😬

Reproducer code:

import numpy as np
import pandas as pd
import polars as pl
from skrub._column_associations import _cramer_v_matrix
from tqdm import tqdm
from time import time
import seaborn as sns

results = []
n_cols_span = 10 ** np.arange(1, 4)
for px in pd, pl:
    for n_cols in tqdm(n_cols_span):
        X = px.DataFrame(np.random.randn(1000, n_cols))
        tic = time()
        _cramer_v_matrix(X)
        duration = time() - tic
        results.append(
            dict(
                duration=duration,
                n_cols=n_cols,
                module=px.__name__,
            )
        )

ax = sns.barplot(
    pd.DataFrame(results),
    x="n_cols",
    y="duration",
    hue="module",
)
ax.set_ylabel("duration (s)")
ax.set_xlabel("n cols (when n rows = 1000)")

[jeromedockes]

thanks a lot for doing this investigation @Vincent-Maladiere . I guess convert to pandas is ok for now, though I would like to limit the number of things we do that require pandas to work in case one day we want to make it optional. another option would be to keep using polars.corr() and if it ends up full of nan that's ok, and maybe one day polars will add an option to handle nans differently. as you prefer

I think that as for plotting, we should skip computing associations when there are many columns by default, but that's for another PR. as you show anyway the cramer statistic is much more expensive to compute than correlations, which is expected, so this PR doesn't really weigh on that

[Vincent-Maladiere]

As I feared, dispatching was not straightforward, and I ran into tons of edge cases

[jeromedockes]

As I feared, dispatching was not straightforward, and I ran into tons of edge cases

ouch sorry about that @Vincent-Maladiere I did not have this foresight and was more optimistic :/ . thanks so much for pushing through, though!! I'll do a final pass tomorrow but this looks like it is ready 🚀

[Vincent-Maladiere]

No worries, I learned a lot while doing it

[jeromedockes]

great, thank you very much @reshamas and @Vincent-Maladiere ! seeing correlations is one of the first things we want so this will be appreciated