PAI Calculation Issue

[norbertgerena]

Good day. Thank you so much for this library.

I just want to clarify the PAI calculation used here. In the original paper, they were predicting X using Mxy, however, in your code, you are predicting Y using Mxy.

Original code inside pai function:

Mₓy = crossmapembed(x, y, d, τ, PAIEmbedding()) 

# Compute `ỹ | Mₓy` and return `correspondence_measure(ỹ | Mₓy, y).
crossmap_basic(Mₓy, y, d, τ; correspondence_measure = correspondence_measure, r = r)

It should be:

Mₓy = crossmapembed(x, y, d, τ, PAIEmbedding()) 

# Compute `x | Mₓy` and return `correspondence_measure(x | Mₓy, x).
crossmap_basic(Mₓy, x, d, τ; correspondence_measure = correspondence_measure, r = r)

Also, If you look closely at crossmap_basic under utils, the corresponding y[t] used is off by (d-1)tau. It should be ` ỹ[i] = sum(w[k]y[j+(d-1)τ] for (k, j) in enumerate(J))` since the indices in J will start from 1 even if it is really referring to (d-1)tau to L timesteps.

[kahaaga]

Hi there, @norbertgerena! Thanks for letting us know, and for the detailed feedback. I will have a look at this immediately.

[norbertgerena]

Thank you @kahaaga. I was trying to replicate the result of the original PAI paper, specifically figure 6b and 6d, but I could not replicate it, even after using the codes of the authors.

Hopefully, they respond to my pull request.

[kahaaga]

As part of addressing this issue, I'll try to reproduce the same figure myself. That would be nice to have in the documentation anyways.

[kahaaga]

@norbertgerena It's been some time since I looked at the source code for these methods, so I just implemented both the cross mapping (CCM) and pairwise asymmetric inference (PAI) methods from scratch again, instead of trying to debug.

• For PAI, I use the embedding given at the bottom of the page under section IV in McCracken & Weigel (2014), i.e. the extended mixed embedding is the regular source embedding, but with one extra variable added: a non-lagged version of the target variable.
• Neither Sugihara et al nor McCracken & Weigel outline how they deal with short time series segments, nor what they consider an "acceptable" lower-limit time series length. I don't do any extra considerations in this edge case, so you'll see some erratic behaviour in the plot below for the shortest library sizes.

Reproducing original figures

In the attached screenshot, I use the new implementation to reproduce figure 3 in Sugihara et al (2012) (first row), and McCracken & Weigel (2014) (second row). I use the same data for both methods. Here, I use embedding dimension d = 2 for CCM, and embedding dimension d = 2 + 1 for PAI (+1 due to the extra variable included in the embedding).

Note:

• Sugihara et al. (2012) give parameters used for their figure 3 in the supplement, but erroneously refer to it as Figure 2. This was a bit confusing.
• Sugihara et al. (2012) state that Figfs 3C-3D are for an asymmetrically coupled system, but the system they present is strictly unidirectional (βxy = 0.0). I guess you can argue that unidirectional is also asymmetrical, but with one-way coupling, I'm not able to reproduce their figures. In what you see below, I use βxy = 0.03 and βyx = 0.32.
• McCracken and Weigel also state in their paper that they use βxy = 0.0 when they reproduce Fig 3 from Sugihara et al. I don't know how to reconcile this with the fact that I can't reproduce Fig 3 myself with those parameters. Anyways, here I just adjusted the system to have weak bidirectional coupling.

Note: in the upper row, I used source-only embeddings. In the lower row, I use mixed embeddings. I just forgot to change the notation

As expected, for PAI the correlation is much stronger when the target variable is also included. Not do the predictions match the observations much more closely (the last two columns of the figure), but convergence is also much faster (first column of the figure).

A note on CCM/PAI

It's been quite a few years since I did a deep literature dive into the CCM & friends world, so I don't know what current state of the art is. However, I'd like to note a few informal things regarding the relationship between PAI and CCM.

Let's reiterate what PAI does. Denote ŷ | X̂ the predictions we make by finding neighbors in the embedding X, and ŷ | X̂Y the predictions we make by finding neighbors in the mixed embedding XY. As I see it, the take-away from the McCracken & Weigel paper is that

ŷ | X̂Y is more strongly correlated to the observations y than ŷ | X

But of course predictions of Y get better when I include Y in the set of variables I use to make predictions. That's like saying "does my diet contain more carrots when I explicitly bring a carrot to every meal, compared to when I just eat whatever is served" (poor analogue, sorry :D). The answer is always yes.

I think a useful way to think about the PAI statistic is in terms of degradation of prediction performance. If I used ŷ | Y, my predictions would probably be near perfect every time (controlled by number of neighbors, noise, etc, of course). If I used ŷ | X (as in CCM) the quality of predictions would, among many factors, be controlled by how strongly X and Y and coupled. Clearly, predictions are worse for ŷ | X than for ŷ | Y, unless X and Y are completely synchronized.

The moment I start including additional variables from X by doing ŷ | YX, performance starts to degrade.

• If X and Y are not coupled, then information from X will degrade predictions of Y. In the worst possible case, ŷ | YX ≈ ŷ | noise, and in any case, ŷ | YX will yield much worse predictions than ŷ | Y.
• When X and Y are coupled, including information about X will contain a mixture between relevant and irrelevant information. Therefore predictions ŷ | YX are less precise than ŷ | Y, but we can still obtain reasonably good predictions (quality(ŷ | YX) > quality(ŷ | noise).
• As X and Y become increasingly coupled, more and more relevant information information is included, so predictions degrade ŷ | YX degrade less and less. In the extreme case, when coupling gets strong enough (complete synchronization), then ŷ | YX = ŷ | Y = x̂ | XY = x̂ | X
• If coupling is stronger in one direction than the other, then this will be reflected in the directional statistic Δ_pai = ŷ | X̂Y - x̂ | YX.

For some reason, the directional statistic Δ_pai = ŷ | X̂Y - x̂ | YX performs better than Δ_ccm = ŷ | X̂ - x̂ | Y as an indicator of directional influence, for the particular systems investigated in McCracken and Weigel (2014).

This is of course not rigorous in any way, but it helped me think about the algorithm when re-implementing it now.

Summary

If you're happy with how the reproduced figures look, @norbertgerena, I'll prepare a PR with the updated code and some examples for the documentation.

[norbertgerena]

Hi @kahaaga. Thank you. I'm happy with this.

I was able to reproduce fig 3 but I'm using Python (ccm_causlity package). To be transparent, I am trying to create a simple PAI library as well (not as elegant and comprehensive as yours :D)

Due to the inconsistencies in the generation of embeddings (when looking at different libraries, including the one from the PAI authors, https://github.com/jmmccracken/PAI/pull/2#issue-1467717305), I was planning to make the embedding more general. I have not finalized it yet but here's a snapshot:

[Datseris]

why not join efforts and develop the best library of all here ;)

[norbertgerena]

Hi @Datseris. I would love to! But first, I need to study Julia :D

[Datseris]

here is a very good tutorial to get started! I am totally unbiased in my opinion of it!

https://www.youtube.com/watch?v=Fi7Pf2NveH0

[kahaaga]

@norbertgerena Given that you've been able to debug the Julia code here, and that you're already programming in Python, I don't think the transition would be difficult at all!

How about this:

When I've got my new implementation of the CCM/PAI code above ready, you can have some time to inspect/understand it. Then, we do the following:

• We do Zoom call and discuss how we can introduce completely general embeddings into the existing code
• Based on our discussion, you work on and submit a proposal for the general embedding PAI in a PR to this repo.
• Myself, and perhaps @Datseris, if he has time, review the code.
• Fix any design/code issues based on feedback
• Merge the new code main branch and release a new version of the package.

[norbertgerena]

@Datseris, subscribed! 😄

@kahaaga, I would love to contribute. Sounds like a plan.

[norbertgerena]

@kahaaga

Your library has been mentioned by the PAI authors themselves :)

[kahaaga]

Your library has been mentioned by the PAI authors themselves :)

Ah, nice! Looks like I need to get the new implementation out, so we can start fine-tuning it with the generic embeddings :D I'll try to have a basic, corrected version ready by tomorrow evening (CET time), as a last effort before Christmas holidays. If not done by then, it'll have to wait until mid next week.

Then I'll just make a new minor release which fixes the bug you found @norbertgerena , and we go from there with the generalized version, which I imagine we release as part of CausalityTools v2.

[kahaaga]

Maybe this should be in a separate issue, but just as a self-reminder for when implementing the generic embeddings:

PAI, by inclusion of the target variable in the source embedding, creates much stronger correlations (artifically, one may argue). I haven't tested it in detail, but perhaps there is a potentially new method here:

Weighting the different variables in the embedding according to some user-specified weights. That is: instead of letting the weights for the linear combinations strictly be a distance of the distance to nearest neighbor, we can make a double-weighting that also considers the "relative importance" of each variable, as specified by the user. This is a straight-forward extension of CCM/PAI that I think may be useful in practice for users that have a-priori information on the trustworthyness of each measured variable.

[kahaaga]

Hey, @norbertgerena! I think I'm just about done with the re-write of the cross mapping API (see #179).

In that process, I re-did the existing examples and added some more examples from Sugihara et al. and McCracken & Weigel. The updated docs are here.

As you pointed out, the embedding I used in the previous version was not the same as they used in the paper. I have now corrected this, so that the method does the same as in the paper (although it is in principle completely arbitrary which mixed embedding to use).

For the CCM stuff, I'm able to reproduce the original papers fairly well. But for the PAI stuff, I can't reproduce their figure 8 (upper right panel). As I understand it, they define Δ' = ρ(ŷ | yx) - ρ(x̂ | xy), where in

• ŷ | yx, ŷ are the predictions of y constructed from the embedding points (y(i), y(i-1), y(i-2), ..., x(i), and
• x̂ | xy x̂ are the predictions of x constructed from the embedding points (x(i), x(i-1), x(i-2), ..., y(i), and

Or, equivalently stated:

• ŷ | yx measures how well an embedding of y predicts y itself, if a single non-lagged component of x is added
• x̂ | xy measures how well an embedding of x predicts x itself, if a single non-lagged component of y is added

McCracken & Weigel's paper is a bit vague about this, but in the preprint version of their paper they state

"If Δ' < 0, then single time step of Y added to the delay vectors constructed from X create stronger estimators of X than the single time step of X added to the delay vectors constructed from Y do for Y"

I'm getting similar results, but with the opposite sign. That is, reproducing their Fig 8, I'm consistently getting Δ' > 0, which would imply that Y drives X, not the other way around.

The fact that the sign is consistent tells me that either 1) I've made a mistake or misunderstood their explanations, or 2) they've flipped the sign in the paper.

Have you tried to reproduce that figure using your Python implementation? It would be nice to see if we're getting the same results before I publish the new version.

[kahaaga]

For the CCM stuff, I'm able to reproduce the original papers fairly well. But for the PAI stuff, I can't reproduce their figure 8 (upper right panel).

Nevermind, I had just made a mistake in the example calculations. I'll merge the bug-fix as soon as CI tests+docs pass.

[norbertgerena]

Happy new year! Sorry I’m still on vacation. Will look into this ASAP.Regards,NorbertOn Dec 27, 2022, at 11:19 PM, Kristian Agasøster Haaga @.***> wrote: Closed #177 as completed via #179.

—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: @.***>

[kahaaga]

Hey @norbertgerena! No worries. Just leave a comment on #184 when you're back when vacation and are available. Happy new year!

[norbertgerena]

Hi @kahaaga. Here are my results when reproducing figures 3 and 8 using our Python package.

Here's for CCM (fig 3):

And here's for PAI (fig 8, using $Y_t $| ( $Yt$, $Y{t-1}$, $X_t$)):

J. McCracken updated his PAI library to address the issue I raised. I'll check it out first then work on #184.

[kahaaga]

@norbertgerena Excellent! Thank you so much for making these reproductions.

I will try to do reproduce all the same plots. Here's my reproduction of the upper-right panel of Figure 8 in McCracken & Weigel (you can find this example in the CausalityTools docs)

Before I can work on the remaining examples, I have I some urgent work that needs to be done first relating to the release of v2 of this package, so it will probably not happen until the weekend or next week. In the meantime, feel free think on how to elegantly deal with #184, and if you have some knowledge on what the new update in the original PAI repo entails, that would be excellent.