It looks like you're new here. If you want to get involved, click one of these buttons!

- All Categories 2.4K
- Chat 502
- Study Groups 21
- Petri Nets 9
- Epidemiology 4
- Leaf Modeling 2
- Review Sections 9
- MIT 2020: Programming with Categories 51
- MIT 2020: Lectures 20
- MIT 2020: Exercises 25
- Baez ACT 2019: Online Course 339
- Baez ACT 2019: Lectures 79
- Baez ACT 2019: Exercises 149
- Baez ACT 2019: Chat 50
- UCR ACT Seminar 4
- General 75
- Azimuth Code Project 110
- Statistical methods 4
- Drafts 10
- Math Syntax Demos 15
- Wiki - Latest Changes 3
- Strategy 113
- Azimuth Project 1.1K
- - Spam 1
- News and Information 148
- Azimuth Blog 149
- - Conventions and Policies 21
- - Questions 43
- Azimuth Wiki 718

## Comments

Thanks Dara!! I am becoming cross-eyed in my old age :)

Does MS==mean squared error? That seems strange given the scatter plots and cross correlation plots. How are you calculating these errors?

`>> Great catch Blake!! > I just wanted to note that Dara made the catch, not me! Thanks Dara!! I am becoming cross-eyed in my old age :) > I get MS = 0.09 for anom vs. link. > but MS=0.61 for anom vs. link (time shifted by -6) > In other words the time shifted link produces more error. Does MS==mean squared error? That seems strange given the scatter plots and cross correlation plots. How are you calculating these errors?`

yes.

I do linear fit in Mathematica and ask for ANOVA TABLE and I just report the MS.

`>Does MS==mean squared error yes. I do linear fit in Mathematica and ask for ANOVA TABLE and I just report the MS.`

I may be dim on these things but consider this contrived situation:

Suppose that the ENSO signal goes something like sin(2t) where t is years from zero. And then consider there is some other metric that we discover that goes like cos(2t). This definitely shows a lead of $\pi/4$ years, and perfect correlation with that lead applied. One could say that the cos(2t) signal is anticipating the sin(2t) signal and so you could use that as a predictor. So when you see a peak in the cos(2t) signal, you could imagine that the ENSO sin(2t) peak is $\pi/4$ years in the future.

However ... we also know that ENSO is somewhat erratic, so that the cos(2t) signal could potentially "hang" at the peak for an indeterminate amount of time. If it is indeed correlated to ENSO, the predictor of $\pi/4$ years in the future keeps moving to the right as we wait for the cos(2t) to release from the peak. In the end it doesn't buy us anything other than we improve our understanding of measures that are correlated.

That is the important consideration with standing wave dynamics, in that there is likely no new information to be gained. I think this also is related to the Kolmogorov complexity that Dara is getting at elsewhere in a discussion thread. The cos(2t) does not provide a lot more useful information than sin(2t), having only a phase separation.

Please correct me if I am missing something. This is a private chat discussion so feel free to beat me up :)

`I may be dim on these things but consider this contrived situation: Suppose that the ENSO signal goes something like sin(2t) where t is years from zero. And then consider there is some other metric that we discover that goes like cos(2t). This definitely shows a lead of $\pi/4$ years, and perfect correlation with that lead applied. One could say that the cos(2t) signal is anticipating the sin(2t) signal and so you could use that as a predictor. So when you see a peak in the cos(2t) signal, you could imagine that the ENSO sin(2t) peak is $\pi/4$ years in the future. However ... we also know that ENSO is somewhat erratic, so that the cos(2t) signal could potentially "hang" at the peak for an indeterminate amount of time. If it is indeed correlated to ENSO, the predictor of $\pi/4$ years in the future keeps moving to the right as we wait for the cos(2t) to release from the peak. In the end it doesn't buy us anything other than we improve our understanding of measures that are correlated. That is the important consideration with standing wave dynamics, in that there is likely no new information to be gained. I think this also is related to the Kolmogorov complexity that Dara is getting at [elsewhere in a discussion thread](http://forum.azimuthproject.org/discussion/1533/informationnon-increase-a-kolmogorov-complexity-theorem/). The cos(2t) does not provide a lot more useful information than sin(2t), having only a phase separation. Please correct me if I am missing something. This is a private chat discussion so feel free to beat me up :)`

Can you post the code? Preferably in some kind of text form, I have trouble loading CDFs.

`> I do linear fit in Mathematica and ask for ANOVA TABLE and I just report the MS. Can you post the code? Preferably in some kind of text form, I have trouble loading CDFs.`

code for ANOVA

anom and link are the two columns of interest in your csv for the links

`[code for ANOVA](http://files.lossofgenerality.com/anom_link_ml.jpg) anom and link are the two columns of interest in your csv for the links`

Those numbers are different to what I got in R in #41 which were .82 .78 and .75 respectively. They are the "Residual stadard error:" entries near the bottom of each output.

In your output the p-value and F-score indicate that significance is increasing as you go from fit1 to fit3, which seems to contradict fit1 having the lowest residual error, but would be at least qualitatively consistent with the residual errors I got in R.

Could we be misinterpreting the MS column?

`Those numbers are different to what I got in R in [#41](http://forum.azimuthproject.org/discussion/1523/crunch-time/?Focus=13487#Comment_13487) which were .82 .78 and .75 respectively. They are the "Residual stadard error:" entries near the bottom of each output. In your output the p-value and F-score indicate that significance is increasing as you go from fit1 to fit3, which seems to contradict fit1 having the lowest residual error, but would be at least qualitatively consistent with the residual errors I got in R. Could we be misinterpreting the MS column?`

I could write to Techsupport and ask them, I could not be certain.

But all that assumes that there is a relatively ok linear fit, but actually there is not, so I am wondering if we these numbers make any sense at all.

`>Could we be misinterpreting the MS column? I could write to Techsupport and ask them, I could not be certain. But all that assumes that there is a relatively ok linear fit, but actually there is not, so I am wondering if we these numbers make any sense at all.`

Paul, re #153, I get your example, but the anomaly and link strength are unlike the sin-cos example in some critical ways. The example works because both signals are completely deterministic and autonomous, neither of which is the case in for the climate signals. If you added some independent noise to both the signals, and did not know their true functional form (ie that they really are just sin and cos) then averaging the one with the time shifted version of the other would give you better predictions then either corrupted signal alone, since the noises would at least partially cancel. You would still be able tell that they were auto/cross correlated and estimate the lags from the correlation plots.

`Paul, re #153, I get your example, but the anomaly and link strength are unlike the sin-cos example in some critical ways. The example works because both signals are completely deterministic and autonomous, neither of which is the case in for the climate signals. If you added some independent noise to both the signals, and did not know their true functional form (ie that they really are just sin and cos) then averaging the one with the time shifted version of the other would give you better predictions then either corrupted signal alone, since the noises would at least partially cancel. You would still be able tell that they were auto/cross correlated and estimate the lags from the correlation plots.`

I was thinking about this 10-day sampling problem and have been looking at Graham's code which computes the link-strengths and messing around with some things. What would be most useful to people here:

1) The link strength computed everyday from whenever it first can be computed to 2013?

2) The link strength computed on the first of each month?

Let me know and I'll run one tonight. The first will probably take a few hours, the second not so much.

`I was thinking about this 10-day sampling problem and have been looking at Graham's code which computes the link-strengths and messing around with some things. What would be most useful to people here: 1) The link strength computed everyday from whenever it first can be computed to 2013? 2) The link strength computed on the first of each month? Let me know and I'll run one tonight. The first will probably take a few hours, the second not so much.`

The scatter plots and correlation seem to indicate there is at least a weak linear relationship. If there was no linear relationship, we should not be seeing p-values of $10^{-4}$, $10^{-17}$ and $10^{-30}$. Also your fit2 and fit3 models agree with the R models on coefficients and the F-statistic and are consistent on p-values, but your fit1 differs a lot from the corresponding R model in coefficients, but they approximately agree on F-statistic and p-value.

Actually I was forgetting that standard residual error is

rootmean squared error and the squares of my errors agree with your MS values for fit2 and fit3, but not fit1.`> But all that assumes that there is a relatively ok linear fit, but actually there is not, so I am wondering if we these numbers make any sense at all. The scatter plots and correlation seem to indicate there is at least a weak linear relationship. If there was no linear relationship, we should not be seeing p-values of $10^{-4}$, $10^{-17}$ and $10^{-30}$. Also your fit2 and fit3 models agree with the R models on coefficients and the F-statistic and are consistent on p-values, but your fit1 differs a lot from the corresponding R model in coefficients, but they approximately agree on F-statistic and p-value. Actually I was forgetting that standard residual error is *root* mean squared error and the squares of my errors agree with your MS values for fit2 and fit3, but not fit1.`

I see the problem: fit1 is predicting the link strength from the anomaly, the same problem problem you spotted in my models in #145. The low error is because link strength has a smaller variance, so there less error to explain away. F-statistic, p-value and R^2 (if it was computed) still agree because they are scale independent.

`I see the problem: fit1 is predicting the link strength from the anomaly, the same problem problem you spotted in my models in #145. The low error is because link strength has a smaller variance, so there less error to explain away. F-statistic, p-value and R^2 (if it was computed) still agree because they are scale independent.`

The second is more immediately useful since we could just directly plug it into the calculations we are doing now. The first is more generally useful. It is not possible to replicate Ludescher's results from 2 fror example

`> Let me know and I’ll run one tonight. The first will probably take a few hours, the second not so much. > 1) The link strength computed everyday from whenever it first can be computed to 2013? > 2) The link strength computed on the first of each month? The second is more immediately useful since we could just directly plug it into the calculations we are doing now. The first is more generally useful. It is not possible to replicate Ludescher's results from 2 fror example`

Actually I misread your definition of 2. The link strength computed just sampled for the first may be too noisy. We really want the monthly values to be averages. Probably we want 1 and then we can compute the monthly averages from that.

`Actually I misread your definition of 2. The link strength computed just sampled for the first may be too noisy. We really want the monthly values to be averages. Probably we want 1 and then we can compute the monthly averages from that.`

John wrote at

I am currently traing to understand those indices a little better. They seam to be partially based on the wind directions. I have though (again) a problem with the data in particular the entry for

850 mb Trade Wind Index (175°West-140°West) 5°North-5°South Central Pacific:

which displays the following data:

850 MB TRADE WIND INDEX(175W-140W)5N 5S CENTRAL PACIFIC

ORIGINAL DATA

YEAR JAN FEB MAR APR MAY JUN JUL AUG SEP OCT NOV DEC

1979 6.3 10.0 9.8 7.9 7.8 9.2 6.9 6.9 6.7 5.9 8.6 7.9

1980 6.9 7.5 3.8 5.7 5.2 6.8 9.0 11.1 10.4 10.4 9.2 10.6

1981 11.2 8.1 8.9 6.8 9.7 9.1 11.5 9.6 12.1 8.7 9.5 9.7

1982 10.1 8.3 7.5 6.1 5.5 6.4 6.8 6.4 2.7 0.7 -1.3 0.1

1983 1.9 2.4 -0.2 -0.2 2.7 5.8 8.3 10.2 10.3 9.8 7.8 10.8

1984 11.7 11.4 7.1 9.3 8.9 9.8 10.0 9.7 8.0 10.1 8.9 10.1

1985 8.6 13.1 10.3 9.5 10.2 9.1 8.5 9.9 9.8 7.0 6.3 7.3

1986 10.8 8.5 11.2 8.5 7.7 9.6 7.5 6.4 6.6 6.0 4.6 5.8

1987 7.6 5.1 5.0 3.3 2.2 5.2 4.2 5.5 6.8 5.7 6.9 7.5

1988 10.0 8.5 10.0 8.2 10.2 9.8 11.9 10.7 13.4 11.6 12.1 11.2

1989 12.4 10.9 10.8 10.5 10.0 8.1 9.8 8.8 7.5 9.8 7.6 6.3

1990 8.7 6.1 8.1 8.1 8.1 8.8 7.1 9.0 7.5 7.2 6.5 6.9

1991 8.8 9.7 8.0 8.1 6.4 6.4 7.6 8.3 5.6 6.0 3.5 2.0

1992 1.2 3.8 2.4 2.9 5.0 7.5 6.4 8.3 5.8 7.1 4.5 6.0

1993 6.6 6.7 7.4 5.3 5.9 6.1 8.2 6.6 6.5 6.1 6.4 9.3

1994 9.9 10.9 8.5 8.1 7.3 8.7 7.3 7.1 6.2 3.8 3.9 7.4

1995 9.7 7.5 9.3 8.0 8.7 8.7 9.1 9.9 8.8 7.9 8.4 6.8

1996 11.3 12.1 10.7 9.0 9.1 9.5 9.8 8.2 9.3 8.0 8.3 9.9

1997 8.7 11.3 6.8 5.8 5.9 1.3 6.8 3.3 0.7 -3.2 -0.7 2.6

1998 2.3 2.3 1.8 2.5 6.5 8.0 9.8 9.0 8.4 10.2 10.3 11.6

1999 12.4 11.6 10.6 9.4 9.3 8.8 10.7 9.7 7.5 9.2 10.3 11.3

2000 12.4 13.5 10.3 9.4 8.7 9.0 8.2 9.9 7.0 9.2 8.8 9.7

2001 11.4 12.6 9.4 9.2 8.4 8.5 8.0 9.1 6.4 6.7 8.9 8.4

2002 9.3 9.4 10.2 8.0 6.9 8.1 6.1 7.9 5.1 5.6 7.9 7.1

2003 9.2 8.4 8.8 8.8 8.3 6.4 9.2 7.7 7.0 6.6 7.4 10.6

2004 8.3 11.1 10.4 8.8 8.6 9.3 6.6 8.3 7.7 6.8 7.6 6.9

2005 9.9 7.5 8.2 9.2 7.5 8.3 9.4 7.6 8.1 8.5 8.5 10.3

2006 12.1 10.8 10.3 9.1 7.4 9.0 7.2 7.2 8.2 5.2 6.1 9.2

2007 10.6 11.0 10.6 9.7 8.7 10.1 8.5 9.5 11.2 8.7 11.2 11.6

2008 12.5 12.2 10.3 9.5 8.9 9.1 8.6 8.8 9.0 8.3 10.3 11.0

2009 10.7 11.7 9.3 8.8 7.7 6.6 9.3 7.0 8.7 4.3 6.8 5.7

2010 8.9 6.1 8.6 7.1 8.9 8.7 10.2 10.4 9.2 11.1 10.0 13.5

2011 10.9 11.4 11.0 9.4 9.1 9.5 9.9 9.6 10.6 7.3 10.3 11.5

2012 11.5 10.1 10.7 9.2 9.1 7.3 8.8 7.8 8.6 6.6 8.4 8.5

2013 10.0 10.6 9.5 9.1 7.5 8.0 8.4 7.4 8.6 6.6 7.3 8.9

2014 9.9 8.0 7.8 7.5 6.5 8.6 8.5 6.3 7.5 6.3-999.9-999.9

2015-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9 2016-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9 2017-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9 2018-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9 2019-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9 2020-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9

850 MB TRADE WIND INDEX(175W-140W)5N 5S CENTRAL PACIFIC

ANOMALY

together with the following graphic:

Doesn't go together in the sense that the image shows a quantity called STD. DEV. and the data contains trade wind indices.

It is also not that that STD. DEV. is just a typo because like for example in 2005 the smallest value is 7.5, which if the data would be identical would be the lowest red dot below the null line in the middle image. In 2011 however the smallest value is a 7.3 and here the red dots don't even cross the null line. STD. DEV. sounds like standard deviation, but thats a guess. Moreover since the computation of the standard deviation in climate science seems to be at least once in a while different from the usual definition. I am asking myself now: what the hell does this graphic show?!

`John <a href="http://forum.azimuthproject.org/discussion/1523/crunch-time/?Focus=13670#Comment_13670">wrote at </a> >WebHubTel wrote: >> The latest is that there will likely be an El Nino. >Thanks! >Will this arrive in time to save Ludescher et al? That’s not a very important question, but it’s interesting for this talk. >They claimed an El Niño would arrive by the end of this calendar year. The usual definition of El Niño requires that the 3-month running mean of the Niño 3.4 index exceeds 0.5 °C for at least 5 months in a row. That can’t happen this year. >But luckily for Ludescher et al, they were using a nonstandard definition of El Niño! For them, all you need is for the monthly average Niño 3.4 index to exceed 0.5 °C. >So, if the December average is above 0.5 °C, they can pat themselves on the back and claim success. I am currently traing to understand <a href="http://www.cpc.ncep.noaa.gov/products/CDB/Tropics/figt5.gif">those indices</a> a little better. They seam to be partially based on the <a href="http://www.cpc.ncep.noaa.gov/data/indices/">wind directions.</a> I have though (again) a problem with the data in particular the entry for 850 mb Trade Wind Index (175°West-140°West) 5°North-5°South Central Pacific: >> <a href="http://www.cpc.ncep.noaa.gov/data/indices/cpac850">Data</a> >> <a href="http://www.cpc.ncep.noaa.gov/products/CDB/Tropics/figt4.gif">Graphic</a> which displays the following data: 850 MB TRADE WIND INDEX(175W-140W)5N 5S CENTRAL PACIFIC ORIGINAL DATA YEAR JAN FEB MAR APR MAY JUN JUL AUG SEP OCT NOV DEC 1979 6.3 10.0 9.8 7.9 7.8 9.2 6.9 6.9 6.7 5.9 8.6 7.9 1980 6.9 7.5 3.8 5.7 5.2 6.8 9.0 11.1 10.4 10.4 9.2 10.6 1981 11.2 8.1 8.9 6.8 9.7 9.1 11.5 9.6 12.1 8.7 9.5 9.7 1982 10.1 8.3 7.5 6.1 5.5 6.4 6.8 6.4 2.7 0.7 -1.3 0.1 1983 1.9 2.4 -0.2 -0.2 2.7 5.8 8.3 10.2 10.3 9.8 7.8 10.8 1984 11.7 11.4 7.1 9.3 8.9 9.8 10.0 9.7 8.0 10.1 8.9 10.1 1985 8.6 13.1 10.3 9.5 10.2 9.1 8.5 9.9 9.8 7.0 6.3 7.3 1986 10.8 8.5 11.2 8.5 7.7 9.6 7.5 6.4 6.6 6.0 4.6 5.8 1987 7.6 5.1 5.0 3.3 2.2 5.2 4.2 5.5 6.8 5.7 6.9 7.5 1988 10.0 8.5 10.0 8.2 10.2 9.8 11.9 10.7 13.4 11.6 12.1 11.2 1989 12.4 10.9 10.8 10.5 10.0 8.1 9.8 8.8 7.5 9.8 7.6 6.3 1990 8.7 6.1 8.1 8.1 8.1 8.8 7.1 9.0 7.5 7.2 6.5 6.9 1991 8.8 9.7 8.0 8.1 6.4 6.4 7.6 8.3 5.6 6.0 3.5 2.0 1992 1.2 3.8 2.4 2.9 5.0 7.5 6.4 8.3 5.8 7.1 4.5 6.0 1993 6.6 6.7 7.4 5.3 5.9 6.1 8.2 6.6 6.5 6.1 6.4 9.3 1994 9.9 10.9 8.5 8.1 7.3 8.7 7.3 7.1 6.2 3.8 3.9 7.4 1995 9.7 7.5 9.3 8.0 8.7 8.7 9.1 9.9 8.8 7.9 8.4 6.8 1996 11.3 12.1 10.7 9.0 9.1 9.5 9.8 8.2 9.3 8.0 8.3 9.9 1997 8.7 11.3 6.8 5.8 5.9 1.3 6.8 3.3 0.7 -3.2 -0.7 2.6 1998 2.3 2.3 1.8 2.5 6.5 8.0 9.8 9.0 8.4 10.2 10.3 11.6 1999 12.4 11.6 10.6 9.4 9.3 8.8 10.7 9.7 7.5 9.2 10.3 11.3 2000 12.4 13.5 10.3 9.4 8.7 9.0 8.2 9.9 7.0 9.2 8.8 9.7 2001 11.4 12.6 9.4 9.2 8.4 8.5 8.0 9.1 6.4 6.7 8.9 8.4 2002 9.3 9.4 10.2 8.0 6.9 8.1 6.1 7.9 5.1 5.6 7.9 7.1 2003 9.2 8.4 8.8 8.8 8.3 6.4 9.2 7.7 7.0 6.6 7.4 10.6 2004 8.3 11.1 10.4 8.8 8.6 9.3 6.6 8.3 7.7 6.8 7.6 6.9 2005 9.9 7.5 8.2 9.2 7.5 8.3 9.4 7.6 8.1 8.5 8.5 10.3 2006 12.1 10.8 10.3 9.1 7.4 9.0 7.2 7.2 8.2 5.2 6.1 9.2 2007 10.6 11.0 10.6 9.7 8.7 10.1 8.5 9.5 11.2 8.7 11.2 11.6 2008 12.5 12.2 10.3 9.5 8.9 9.1 8.6 8.8 9.0 8.3 10.3 11.0 2009 10.7 11.7 9.3 8.8 7.7 6.6 9.3 7.0 8.7 4.3 6.8 5.7 2010 8.9 6.1 8.6 7.1 8.9 8.7 10.2 10.4 9.2 11.1 10.0 13.5 2011 10.9 11.4 11.0 9.4 9.1 9.5 9.9 9.6 10.6 7.3 10.3 11.5 2012 11.5 10.1 10.7 9.2 9.1 7.3 8.8 7.8 8.6 6.6 8.4 8.5 2013 10.0 10.6 9.5 9.1 7.5 8.0 8.4 7.4 8.6 6.6 7.3 8.9 2014 9.9 8.0 7.8 7.5 6.5 8.6 8.5 6.3 7.5 6.3-999.9-999.9 2015-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9 2016-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9 2017-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9 2018-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9 2019-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9 2020-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9-999.9 850 MB TRADE WIND INDEX(175W-140W)5N 5S CENTRAL PACIFIC ANOMALY together with the following graphic: ![wind anomalies](http://www.cpc.ncep.noaa.gov/products/CDB/Tropics/figt4.gif) Doesn't go together in the sense that the image shows a quantity called STD. DEV. and the data contains trade wind indices. It is also not that that STD. DEV. is just a typo because like for example in 2005 the smallest value is 7.5, which if the data would be identical would be the lowest red dot below the null line in the middle image. In 2011 however the smallest value is a 7.3 and here the red dots don't even cross the null line. STD. DEV. sounds like standard deviation, but thats a guess. Moreover since the computation of the standard deviation in climate science seems to be at least once in a while different from the usual definition. I am asking myself now: what the hell does this graphic show?!`

Blake at #159, Daniel at #162,163

The following code allegedly interpolates the ten-day numbers to get daily numbers then takes monthly averages. It would be better to do this properly as Blake's suggestion (1), but maybe this code is useful for cross-checking.

`Blake at #159, Daniel at #162,163 The following code allegedly interpolates the ten-day numbers to get daily numbers then takes monthly averages. It would be better to do this properly as Blake's suggestion (1), but maybe this code is useful for cross-checking. ~~~~ fpath10 <- "C:/Users/Work/Desktop/average-link-strength-1948-2013.txt" fpathmonths <- "C:/Users/Work/Desktop/monthly-S-1950-2013.csv" x <- read.table(fpath10)[,2] # x[] has 32*73+1 = 2337 values. # That's 64 leap-day-less years at 10 day intervals, plus 1. # Starts from day 730, and goes until day 24090, where day 1 is the first of January 1948. # x[1], ie day 730, is last day of 1949. x[2] is 10th Jan, 1950. .... x[2337] is last day of 2013 # linear interpolation to get daily figures 1950 through 2013 daily.S <- rep(0,64*365) for (d in 1:10) { daily.S[ 10*(0:2335) + d ] <- ( x[1:2336]*(10-d) + x[2:2337]*d ) / 10 } # make a 3-column matrix of monthly figures mlens <- c(31,28,31,30,31,30,31,31,30,31,30,31) csm <- c(0,cumsum(mlens)) monthly.S <- matrix(0, nrow=64*12, ncol=3) ym <- 1 for (y in 1:64) { for (m in 1:12) { monthly.S[ym, 1] <- 1950+y-1 monthly.S[ym, 2] <- m monthly.S[ym, 3] <- mean(daily.S[(y-1)*365 + csm[m]:(csm[m+1]-1)]) ym <- ym+1 } } colnames(monthly.S) <- c("Year", "Month", "S") write.csv(monthly.S, file=fpathmonths, row.names=FALSE) ~~~~`

Blake wrote:

Please do this one! As Daniel noted, we need this to compute monthly average link strengths, to compare these with Niño 3.4's, which are also monthly averages!

If you get this information to us in time for Daniel / Dara to draw scatter plots and compute correlations soon, I can use those in my talk. You can put a file of daily link strengths on github next to the other link strength files or email it to me. Do you need someone to invite you to the Azimuth Project github account?

`Blake wrote: > I was thinking about this 10-day sampling problem and have been looking at Graham’s code which computes the link-strengths and messing around with some things. What would be most useful to people here: > 1) The link strength computed every day from whenever it first can be computed to 2013? Please do this one! As Daniel noted, we need this to compute monthly average link strengths, to compare these with Niño 3.4's, which are also monthly averages! If you get this information to us in time for Daniel / Dara to draw scatter plots and compute correlations soon, I can use those in my talk. You can put a file of daily link strengths [on github next to the other link strength files](https://github.com/azimuth-project/el-nino) or email it to me. Do you need someone to invite you to the Azimuth Project github account?`

Thanks for answering all my questions so clearly and carefully, Daniel. I will not try to fix the graphs that I'd mixed up in my earlier comments. Instead, I'll try to draft the next part of my talk now... but there will be some questions for you and/or Dara.

The numbers and graphs here may change if Blake gives you daily link strengths in time and you redo your calculations using those. But the overall message probably won't change. So:

17) Do high average link strengths between points in the El Niño basin and points elsewhere in the Pacific really increase the chance that an El Niño is coming? It is hard to tell from the work of Ludescher

et al.One reason is that they treat El Niño as a binary condition, either on or off depending on whether the Niño 3.4 index for a given month exceeds 0.5 or not. This is not the usual definition of El Niño, but the real problem is that they are only making a single yes-or-no prediction each year for 65 years: does an El Niño occur during this year, or not? 31 of these years (1950-1980) are used for training their algorithm, leaving just 34 retrodictions and one actual prediction (1981-2013, and 2014).

So, there is a serious problem with small sample size.

We can learn a bit by taking a different approach, and simply running some linear regressions between the average link strength and the Niño 3.4 index for each month. There are 766 months from 1950 to 2013, so this gives us more data to look at. Of course, it's possible that the relation between average link strength and Niño is highly nonlinear, so a linear regression may not be appropriate. But it is at least worth looking at!

Daniel Mahler and Dara Shayda of the Azimuth Project did this and found the following interesting results.

`Thanks for answering all my questions so clearly and carefully, Daniel. I will not try to fix the graphs that I'd mixed up in my earlier comments. Instead, I'll try to draft the next part of my talk now... but there will be some questions for you and/or Dara. The numbers and graphs here may change if Blake gives you daily link strengths in time and you redo your calculations using those. But the overall message probably won't change. So: 17) Do high average link strengths between points in the El Niño basin and points elsewhere in the Pacific really increase the chance that an El Niño is coming? It is hard to tell from the work of Ludescher _et al_. One reason is that they treat El Niño as a binary condition, either on or off depending on whether the Niño 3.4 index for a given month exceeds 0.5 or not. This is not the usual definition of El Niño, but the real problem is that they are only making a single yes-or-no prediction each year for 65 years: does an El Niño occur during this year, or not? 31 of these years (1950-1980) are used for training their algorithm, leaving just 34 retrodictions and one actual prediction (1981-2013, and 2014). So, there is a serious problem with small sample size. We can learn a bit by taking a different approach, and simply running some linear regressions between the average link strength and the Niño 3.4 index for each month. There are 766 months from 1950 to 2013, so this gives us more data to look at. Of course, it's possible that the relation between average link strength and Niño is highly nonlinear, so a linear regression may not be appropriate. But it is at least worth looking at! Daniel Mahler and Dara Shayda of the Azimuth Project did this and found the following interesting results.`

18) Here is a scatter plot showing the Niño 3.4 index as a function of the average link strength

on the same month:The coefficient of determination, $R^2$, is 0.0175. That's quite low! But we get a $p$-value of 0.000237, so the effect is statistically significant: it's just very weak.

Of course, we are not really interested in using the average link strength to predict the Niño 3.4 index

at the same moment in time, but this is worth looking at as a kind of baseline.(The number above is the "multiple $R^2$" from Daniel's comment 41, which is apparently just the same as $R^2$ in the current example, when we're doing linear regression with just one regressor. There is also an "adjusted $R^2$", which "attempts to take account of the phenomenon of the $R^2$ automatically and spuriously increasing when extra explanatory variables are added to the model.

Here's the full story

)

`18) Here is a scatter plot showing the Niño 3.4 index as a function of the average link strength _on the same month_: <img src = "http://math.ucr.edu/home/baez/climate_networks/mahler_nino3.4_versus_link_strength_no_lag.png" alt = ""/> The [coefficient of determination](http://en.wikipedia.org/wiki/Coefficient_of_determination), $R^2$, is 0.0175. That's quite low! But we get a $p$-value of 0.000237, so the effect is statistically significant: it's just very weak. Of course, we are not really interested in using the average link strength to predict the Niño 3.4 index _at the same moment in time_, but this is worth looking at as a kind of baseline. (The number above is the "multiple $R^2$" from Daniel's [comment 41](http://forum.azimuthproject.org/discussion/1523/crunch-time/?Focus=13487#Comment_13487), which is apparently [just the same as $R^2$](http://stats.stackexchange.com/questions/90793/whats-the-difference-between-multiple-r-and-r-squared) in the current example, when we're doing linear regression with just one regressor. There is also an "adjusted $R^2$", which "[attempts to take account of the phenomenon of the $R^2$ automatically and spuriously increasing when extra explanatory variables are added to the model](https://en.wikipedia.org/wiki/Coefficient_of_determination#Adjusted_R2). Here's the [full story](http://forum.azimuthproject.org/discussion/1523/crunch-time/?Focus=13487#Comment_13487) > Residual standard error: 0.8157 on 766 degrees of freedom Multiple R-squared: 0.0175, Adjusted R-squared: 0.01621 F-statistic: 13.64 on 1 and 766 DF, p-value: 0.000237 )`

19) Here is a scatter plot showing the Niño 3.4 index as a function of the average link strength

six months earlier:Now $R^2$ is 0.088. So, the link strength explains 8.8% of the variance in the Niño 3.4 index 6 months later. This is not much - but interestingly, it's much more than when we try to relate them at the same moment in time! And the $p$-value is less than $2.2 \cdot 10^{-16}$, so the effect is statistically significant.

(The full story:

)

`19) Here is a scatter plot showing the Niño 3.4 index as a function of the average link strength _six months earlier_: <img src = "http://math.ucr.edu/home/baez/climate_networks/mahler_nino3.4_versus_link_strength_6_months_earlier.png" alt = ""/> Now $R^2$ is 0.088. So, the link strength explains 8.8% of the variance in the Niño 3.4 index 6 months later. This is not much - but interestingly, it's much more than when we try to relate them at the same moment in time! And the $p$-value is less than $2.2 \cdot 10^{-16}$, so the effect is statistically significant. (The [full story](http://forum.azimuthproject.org/discussion/1523/crunch-time/?Focus=13487#Comment_13487): > Residual standard error: 0.7822 on 760 degrees of freedom Multiple R-squared: 0.08826, Adjusted R-squared: 0.08706 F-statistic: 73.57 on 1 and 760 DF, p-value: < 2.2e-16 )`

20) Of course, we could also try to use Niño 3.4 to predict

itself. Here is a scatter plot showing the Niño 3.4 index as a function of the Niño 3.4 index six months earlier.The coefficient of determination, $R^2$, is 0.1619. So, this is better than using the average link strength.

(The full story:

)

`20) Of course, we could also try to use Niño 3.4 to predict _itself_. Here is a scatter plot showing the Niño 3.4 index as a function of the Niño 3.4 index six months earlier. <img src = "http://math.ucr.edu/home/baez/climate_networks/mahler_nino3.4_versus_nino3.4_6_months_earlier.png" alt = ""/> The [coefficient of determination](http://en.wikipedia.org/wiki/Coefficient_of_determination), $R^2$, is 0.1619. So, this is better than using the average link strength. (The [full story](http://forum.azimuthproject.org/discussion/1523/crunch-time/?Focus=13487#Comment_13487): > Residual standard error: 0.75 on 760 degrees of freedom Multiple R-squared: 0.1619, Adjusted R-squared: 0.1608 F-statistic: 146.8 on 1 and 760 DF, p-value: < 2.2e-16 )`

21) Finally, we could try to predict Niño 3.4 using

bothitselfandthe average link strength 6 months earlier. Here is a scatter plot showing that:Here the $x$ axis is an optimally chosen linear combination of average and link strength and Niño 3.4 - one that maximizes $R^2$.

In this case $R^2 = 0.22$.

(The full story:

)

`21) Finally, we could try to predict Niño 3.4 using _both_ itself _and_ the average link strength 6 months earlier. Here is a scatter plot showing that: <img src = "http://math.ucr.edu/home/baez/climate_networks/mahler_nino3.4_versus_link_strength_and_nino3.4_6_months_earlier.png" alt = ""/> Here the $x$ axis is an optimally chosen linear combination of average and link strength and Niño 3.4 - one that maximizes $R^2$. In this case $R^2 = 0.22$. (The [full story](http://forum.azimuthproject.org/discussion/1523/crunch-time/?Focus=13526#Comment_13526): > Residual standard error: 0.7231 on 759 degrees of freedom Multiple R-squared: 0.222, Adjusted R-squared: 0.22 F-statistic: 108.3 on 2 and 759 DF, p-value: < 2.2e-16 )`

$R^2$ values:

Yes, it is the prediction of the 2 variable linear regression model

`$R^2$ values: + current link -> current ANOM: 0.016 + current link -> future ANOM: 0.087 + current ANOM -> future ANOM: 0.16 + [current link, current ANOM] -> future ANOM: 0.22 --- > (What’s the x axis here??? I guess it’s an optimally chosen linear combination of average and link strength and Niño 3.4 - one that maximizes R2?) Yes, it is the prediction of the 2 variable linear regression model`

Thanks, Daniel! I am filling in the blanks now... there are still question marks in comments 170 and 171, but I will fill those in soon. I have to hit the gym now!

`Thanks, Daniel! I am filling in the blanks now... there are still question marks in comments 170 and 171, but I will fill those in soon. I have to hit the gym now!`

re #168

True, but but the fact that the $R^2$ is so low here is what makes link strength interesting at all. It shows that what little information link strength gives about the future value of the Niño 3.4 index is independent of, and thus additional to, the information contained in the current value of the index. This suggests the link strength is (a proxy for) some external driving force on the Niño 3.4 index. This is also brought out in the link-ANOM xcorr plot.

The xcorrplot also shows a noticeable negative peak at around 14 years in the future. Not sure what that means.

`re #168 > Of course, we are not really interested in using the average link strength to predict the Niño 3.4 index at the same moment in time, but this is worth looking at as a kind of baseline. True, but but the fact that the $R^2$ is so low here is what makes link strength interesting at all. It shows that what little information link strength gives about the future value of the Niño 3.4 index is independent of, and thus additional to, the information contained in the current value of the index. This suggests the link strength is (a proxy for) some external driving force on the Niño 3.4 index. This is also brought out in the link-ANOM xcorr plot. The xcorrplot also shows a noticeable negative peak at around 14 years in the future. Not sure what that means.`

The plot in #171 is the same 6 month lag ANOM-ANOM plot as in #170

`The plot in #171 is the same 6 month lag ANOM-ANOM plot as in #170`

Daniel R^2 being so low means that linear fit failed for the data, as we can see visually. R^2 closer to one means the linear fit was perfect.

Coefficient of determination

I got the same R^2 from Mathematica finally and our results agree, however R^2 being so low for

current link -> future ANOMandcurrent link -> current ANOMmeans that link did no give a good linear fit.This is congruence to my k-NN findings earlier

`>True, but but the fact that the R^2 is so low here is what makes link strength interesting at all. Daniel R^2 being so low means that linear fit failed for the data, as we can see visually. R^2 closer to one means the linear fit was perfect. [Coefficient of determination](http://en.wikipedia.org/wiki/Coefficient_of_determination#Interpretation) I got the same R^2 from Mathematica finally and our results agree, however R^2 being so low for **current link -> future ANOM** and **current link -> current ANOM** means that link did no give a good linear fit. This is congruence to my k-NN findings earlier`

True, but the point is that the $R^2$ for link -> future ANOM is more than 5 times higher than for link -> current ANOM. This means that around 80% of the current link -> future ANOM $R^2$ gets added to the current ANOM -> future ANOM $R^2$ when we build a model combining them. As a result we can account for 22% of of the variance of ANOM 6 months in the future with just two numbers (significantly better than the 16% with we get with current ANOM alone). Given the complexity of climate systems this seems interesting. We do not expect anything near perfect prediction of ANOM 6 months in the future using only past ANOM and/or link strength values. Using a full grid of pressure values only gets around 0.37.

`> I got the same R^2 from Mathematica finally and our results agree, however R^2 being so low for current link -> future ANOM and current link -> current ANOM means that link did no give a good linear fit. True, but the point is that the $R^2$ for link -> future ANOM is more than 5 times higher than for link -> current ANOM. This means that around 80% of the current link -> future ANOM $R^2$ gets added to the current ANOM -> future ANOM $R^2$ when we build a model combining them. As a result we can account for 22% of of the variance of ANOM 6 months in the future with just two numbers (significantly better than the 16% with we get with current ANOM alone). Given the complexity of climate systems this seems interesting. We do not expect anything near perfect prediction of ANOM 6 months in the future using only past ANOM and/or link strength values. Using a full grid of pressure values only gets around 0.37.`

R^2 formula is not accumulative nor linear, so the

5 timesconcept in terms of scaling might not be valid, basically you say the linear fit is very bad but much better than another much worse fit.I have not seen any place people taking two R^2 and divide them by each other to compare.

However,you should take the opinion of a statistician, my guy feel from the math for R^2 is your argument could be countered and questioned, again take what I say with caution.

`>True, but the point is that the R^2 for link -> future ANOM is more than 5 times higher than for link -> current ANOM R^2 formula is not accumulative nor linear, so the **5 times** concept in terms of scaling might not be valid, basically you say the linear fit is very bad but much better than another much worse fit. I have not seen any place people taking two R^2 and divide them by each other to compare. However,you should take the opinion of a statistician, my guy feel from the math for R^2 is your argument could be countered and questioned, again take what I say with caution.`

My interpretation of the R^2 numbers is: there is no good linear fit which we could see visually and the link strength did not aid in bettering the linear fit at all. k-NN regression shows that link strength could better the forecasts but there are other very simple patterns of numbers which are even better in forecast, in other words the link strength is not delivering better forecasts as claimed by the paper in question.

`My interpretation of the R^2 numbers is: there is no good linear fit which we could see visually and the link strength did not aid in bettering the linear fit at all. k-NN regression shows that link strength could better the forecasts but there are other very simple patterns of numbers which are even better in forecast, in other words the link strength is not delivering better forecasts as claimed by the paper in question.`

$R^2$ measures the fraction of variance of the dependent variable explained by the independent variable and variances are additive.

Writing $E[X]$ for the mean of random variable $X$ and using the linearity of $E$:

$Var(X) \triangleq E[X^2] - E[X]^2$

$Var(X+Y)=E[(X+Y)^2] - E[X+Y]^2$

and if $X$ and $Y$ are independent then $Var(X+Y)=Var(X) + Var(Y)$

For additivity of $R^2$ suppose there are 5 mutually independent variables $A, B, C, D , E$ and 3 variables $X, Y, Z$ with

If $A, B, C, D, E$ are all Gaussian then there are optimal reression models:

then

This should still approximately hold when $R^2_{X|Y}$ is small compared to $R^2_{Z|X}$ and $R^2_{Z|Y}$, as it is in this case.

`> R^2 formula is not accumulative nor linear $R^2$ measures the fraction of variance of the dependent variable explained by the independent variable and variances are additive. Writing $E[X]$ for the mean of random variable $X$ and using the linearity of $E$: + $Var(X) \triangleq E[X^2] - E[X]^2$ + $Var(X+Y)=E[(X+Y)^2] - E[X+Y]^2$ + $=(E[X]^2 + E[XY] + E[Y]^2) - (E[X]^2 + E[X]E[Y] + E[Y]^2)$ + $=(E[X]^2 - E[X]^2) + (E[Y]^2 - E[Y]^2) + (E[XY] - (E[X]E[Y])$ + $=Var(X) + Var(Y) + Cov(X,Y)$ and if $X$ and $Y$ are independent then $Var(X+Y)=Var(X) + Var(Y)$ For additivity of $R^2$ suppose there are 5 mutually independent variables $A, B, C, D , E$ and 3 variables $X, Y, Z$ with + $X = A/k + B$ + $Y = C/l + D$ + $Z = A + C + E$ If $A, B, C, D, E$ are all Gaussian then there are optimal reression models: + $\hat{Z} = k X - k E[X] + E[Z] = k A + E[C] + E[E]$ + $\hat{Z} = k Y - l E[Y] + E[Z] = l C + E[A] + E[E]$ + $\hat{Z} = k X + l Y - k E[X] - l E[Y] + E[Z] = k A + l C + E[E]$ then + $R^2_{Z|X} = Var(A)/(Var(A)+Var(C)+Var(E))$ + $R^2_{Z|Y} = Var(C)/(Var(A)+Var(C)+Var(E))$ + $R^2_{Z|XY} = (Var(A)+Var(C))/(Var(A)+Var(C)+Var(E))$ + $= R^2_{Z|X} + R^2_{Z|Y}$ This should still approximately hold when $R^2_{X|Y}$ is small compared to $R^2_{Z|X}$ and $R^2_{Z|Y}$, as it is in this case.`

Correlation between link and anom is 0.134 which is substantial and the cov between link and anom is 0.034 which is again substantial, therefore the linearity for Var fails.

anom and link are dependent vars, but that does not guarantee that one could help the other in forecast. Example: the rotary motion of the bicycle wheels are dependent on the path on the road, but they cannot (alone) help in forecasting where the bike goes.

`Correlation between link and anom is 0.134 which is substantial and the cov between link and anom is 0.034 which is again substantial, therefore the linearity for Var fails. anom and link are dependent vars, but that does not guarantee that one could help the other in forecast. Example: the rotary motion of the bicycle wheels are dependent on the path on the road, but they cannot (alone) help in forecasting where the bike goes.`

Watch how one index seems to lead another index, moving east to west. Concentrate on the middle peak around 2012.

`Watch how one index seems to lead another index, moving east to west. Concentrate on the middle peak around 2012. ![ninos](http://imageshack.com/a/img540/2896/nyoT2u.gif)`

Very interesting. That there are positive Nino2 temperature anomalies in the east about 15 months before the more westerly Nino4 seems to contradict the idea that surface warm water is sloshing from east to west as a simple flattening of the thermocline tilt.

So where does the heat in Nino2 come? Is it plausible that there's subduction of deep western warm water under a cooler central Pacific Nino4 and then upwelling in the east?

We know that there is a transfer of water from ocean to land, raining on Indo-China and Australia at a peak La Nina. In most accounts westerlies pile up water in the west and as with sand dunes, blows water westwards off the top. I'd expect sea levels to be lowered in the central and eastern Pacific. So perhaps the source of easterly winds and the shift in the Walker circulation is pulled rather that pushed?

`Very interesting. That there are positive Nino2 temperature anomalies in the east about 15 months before the more westerly Nino4 seems to contradict the idea that surface warm water is sloshing from east to west as a simple flattening of the thermocline tilt. So where does the heat in Nino2 come? Is it plausible that there's subduction of deep western warm water under a cooler central Pacific Nino4 and then upwelling in the east? We know that there is a transfer of water from ocean to land, raining on Indo-China and Australia at a peak La Nina. In most accounts westerlies pile up water in the west and as with sand dunes, blows water westwards off the top. I'd expect sea levels to be lowered in the central and eastern Pacific. So perhaps the source of easterly winds and the shift in the Walker circulation is pulled rather that pushed?`

Here is some code that simulates the idea.

The function

`testr2`

starts with uncorrelated variables A, B, C, D, E, F and uses them to construct correlated variables X, Y, Z.The variables X, Y, Z all share a dependency on F, but the weight of F in X and Y is small compared with A, B, C, D, but F has a normal weight in Z.

The function then constructs the models X->Y, X->Z, Y->Z, [X,Y]->Z and returns their respective $R^2$ values.

Now run testr2 200 times and build a linear model relating the $R^2$ values of the models:

The result should look something like

This shows that there is an almost perfect fit for the fomula

$R^2_{Z|XY} = R^2_{Z|X} + R^2_{Z|Y} - k R^2_{Y|X}$, $k ≤ 1$

This was the reasoning behind the statement in #177

`> Correlation between link and anom is 0.134 which is substantial and the cov between link and anom is 0.034 which is again substantial, therefore the linearity for Var fails. > anom and link are dependent vars, but that does not guarantee that one could help the other in forecast. Example: the rotary motion of the bicycle wheels are dependent on the path on the road, but they cannot (alone) help in forecasting where the bike goes. Here is some code that simulates the idea. The function`testr2` starts with uncorrelated variables A, B, C, D, E, F and uses them to construct correlated variables X, Y, Z. The variables X, Y, Z all share a dependency on F, but the weight of F in X and Y is small compared with A, B, C, D, but F has a normal weight in Z. The function then constructs the models X->Y, X->Z, Y->Z, [X,Y]->Z and returns their respective $R^2$ values. testr2 <- function(w) { # independent components A <- rnorm(1000) B <- rnorm(1000) C <- rnorm(1000) D <- rnorm(1000) E <- rnorm(1000) F <- rnorm(1000) # common to X, Y and Z # mixing coefficients ax <- rnorm(1) bx <- rnorm(1) fx <- w*rnorm(1)*min(abs(ax),abs(bx)) cy <- rnorm(1) dy <- rnorm(1) fy <- w*rnorm(1)*min(abs(cy),abs(dy)) az <- rnorm(1) cz <- rnorm(1) ez <- rnorm(1) fz <- rnorm(1) # correlated "observed" variables X = ax * A + bx * B + fx * F Y = cy * C + dy * D + fy * F Z = az * A + cz * C + ez * E + fz * F # linear regression between observed variables XY <- lm(Y ~ X) ZX <- lm(Z ~ X) ZY <- lm(Z ~ Y) ZXY <- lm(Z ~ X + Y) # function to return R^2 of a model r2 <- function(model) { summary(model)$r.squared } #return the list of R^2 values of models which should obey r2(ZXY) ~=~ r2(ZX) + r2(ZY) - k r2(XY), with k ~ O(1) sapply(list(XY=XY,ZX=ZX,ZY=ZY,ZXY=ZXY),r2) } Now run testr2 200 times and build a linear model relating the $R^2$ values of the models: d <- as.data.frame(t(as.matrix(replicate(200, testr2(0.2))))) summary(lm(ZXY ~ ., d)) The result should look something like Call: lm(formula = ZXY ~ ., data = d) Residuals: Min 1Q Median 3Q Max -0.031957 -0.001208 -0.000439 0.000682 0.032017 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.0007550 0.0006743 1.12 0.264 XY -0.3000197 0.2969204 -1.01 0.314 ZX 0.9991259 0.0024138 413.92 <2e-16 *** ZY 0.9988262 0.0023788 419.88 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.005993 on 196 degrees of freedom Multiple R-squared: 0.9994, Adjusted R-squared: 0.9993 F-statistic: 1.017e+05 on 3 and 196 DF, p-value: < 2.2e-16 This shows that there is an almost perfect fit for the fomula $R^2_{Z|XY} = R^2_{Z|X} + R^2_{Z|Y} - k R^2_{Y|X}$, $k ≤ 1$ This was the reasoning behind the statement in #177 > True, but the point is that the R2 for link -> future ANOM is more than 5 times higher than for link -> current ANOM. This means that around 80% of the current link -> future ANOM R2 gets added to the current ANOM -> future ANOM R2 when we build a model combining them. As a result we can account for 22% of of the variance of ANOM 6 months in the future with just two numbers (significantly better than the 16% with we get with current ANOM alone). Given the complexity of climate systems this seems interesting. We do not expect anything near perfect prediction of ANOM 6 months in the future using only past ANOM and/or link strength values. Using a full grid of pressure values only gets around 0.37.`

Hello Daniel this is the issue:

You made your vars completely random i.e. independent of each other, but anom and link are not independent they have corr and var values that are substantial.

Therefore your example is not a fair comparison.

`Hello Daniel this is the issue: >A <- rnorm(1000) You made your vars completely random i.e. independent of each other, but anom and link are not independent they have corr and var values that are substantial. Therefore your example is not a fair comparison.`

Hi Dara. The regressions are done between X, Y, Z, which are all correlated to some degree. The A..F are the independent components out of which X, Y, Z are built.

`Hi Dara. The regressions are done between X, Y, Z, which are all correlated to some degree. The A..F are the independent components out of which X, Y, Z are built.`

Dara, I have added some comments to the code in #184 to cllarify what is being done.

`Dara, I have added some comments to the code in #184 to cllarify what is being done.`

Thanx Daniel, let's end this particular discussion here and your comments and code are most valuable. Generally this is how the inferences are made by the code and output from the code and carefully postmortem study, not by loads and loads of English sentences making arguments.

I was trying to help with the stats, but I prefer it is done by you or John himself, to be on the safe side.

Let's see what else John needs

`Thanx Daniel, let's end this particular discussion here and your comments and code are most valuable. Generally this is how the inferences are made by the code and output from the code and carefully postmortem study, not by loads and loads of English sentences making arguments. I was trying to help with the stats, but I prefer it is done by you or John himself, to be on the safe side. Let's see what else John needs`

Say you had two measures: Result A was the measure of maximum El Nino and what people are interested in predicting. But then someone discovered a separate Result B, which is another measure that always precedes Result A by a few months.

But eventually someone found out that B was simply the time derivative of A, or the rate of change of A. Would we consider that a good forecasting measure?

What do you say in response to someone that brings up such a pedantic question? I am always interested to know if I am just fooling myself with apparently interesting findings.

Incidentally, this comes up in electronics where a measure of current can lead a measure of voltage (and vice versa). But yet again, is this really forecasting anything? If the idea was to forecast a dangerous electrical shock potential, then perhaps one monitors the current rushing in to a discharge capacitor. One can think of placing a peak detector in that case, warning the person of imminent shock potential. But what happens if the voltage keeps rising and the peak is delayed? That peak detector does us no good, because the shock is made worse. The early warning sign with shock hazards is feeling your hair standing on its end, which is just a sensitive detector of increasing electrostatic voltage building up.

`Say you had two measures: Result A was the measure of maximum El Nino and what people are interested in predicting. But then someone discovered a separate Result B, which is another measure that always precedes Result A by a few months. ![measure](http://imageshack.com/a/img537/4002/XFQW8l.gif) But eventually someone found out that B was simply the time derivative of A, or the rate of change of A. Would we consider that a good forecasting measure? What do you say in response to someone that brings up such a pedantic question? I am always interested to know if I am just fooling myself with apparently interesting findings. Incidentally, this comes up in electronics where a measure of current can lead a measure of voltage (and vice versa). But yet again, is this really forecasting anything? If the idea was to forecast a dangerous electrical shock potential, then perhaps one monitors the current rushing in to a discharge capacitor. One can think of placing a peak detector in that case, warning the person of imminent shock potential. But what happens if the voltage keeps rising and the peak is delayed? That peak detector does us no good, because the shock is made worse. The early warning sign with shock hazards is feeling your hair standing on its end, which is just a sensitive detector of increasing electrostatic voltage building up.`

Way back on November 13th, in comment #24 of this thread, David Tanzer wrote:

That's right. For my talk I'm sticking with their original method - or more precisely, Graham's interpretation of what they must have meant when they explained their original method.

I think this is fine for now, since

1) we know for sure that the formula for running means described in their paper could not have been the formula they actually used,

2) Graham's way of fixing the formula makes lots of sense and could easily be what they meant,

3) the results he got with his fixed version rather closely match theirs (though not perfectly).

I think it's a pretty minor issue compared to many others we're dealing with. There are some really huge questions about their approach I will try to (gently) raise in my talk.

Okay, I'll do this sometime. For some reason I'm dreading it, I guess because introducing myself while simultaneously raising touchy points isn't something I enjoy. But it would be good to resolve these issues before we publish a paper on this stuff, if we decide to to that.

`Way back on November 13th, in comment [#24](http://forum.azimuthproject.org/discussion/1523/crunch-time/?Focus=13443#Comment_13443) of this thread, David Tanzer wrote: > There was something wrong with their notation for the running means, which Nad picked up on. She wrote to the authors to ask about this, but they didn’t respond. Graham ended up making an assumption about what they actually meant. On that basis, he closely replicated their results. Then he simplified their measures, and still obtained similar results. That's right. For my talk I'm sticking with their original method - or more precisely, Graham's interpretation of what they must have meant when they explained their original method. I think this is fine for now, since 1) we know for sure that the formula for running means described in their paper could not have been the formula they actually used, 2) Graham's way of fixing the formula makes lots of sense and could easily be what they meant, 3) the results he got with his fixed version rather closely match theirs (though not perfectly). > Graham and John, can you say how much of these issues have a significant bearing on the comparison between our data and theirs? And how much of it limits our ability to make a deeper critique of their methodology. I think it's a pretty minor issue compared to many others we're dealing with. There are some really huge questions about their approach I will try to (gently) raise in my talk. > I had posted a draft of a letter to Ludescher et al to address some of these questions, and discussed it on the forum here. John if you haven’t done so already, I suggest that you write this letter, because you are going to know how to put the question most clearly and in a professional academic way. Okay, I'll do this sometime. For some reason I'm dreading it, I guess because introducing myself while simultaneously raising touchy points isn't something I enjoy. But it would be good to resolve these issues before we publish a paper on this stuff, if we decide to to that.`

Daniel wrote:

Damn! Thanks! I fixed it: it was a cut-and-paste error: I forgot to change the image link in #170 when writing #171.

(I'm always trying to analyze the causes of these stupid mistakes I make. In this particular case I was in a rush, trying to finish some work before meeting my wife to go to the gym.)

`Daniel wrote: > The plot in #171 is the same 6 month lag ANOM-ANOM plot as in #170. Damn! Thanks! I fixed it: it was a cut-and-paste error: I forgot to change the image link in #170 when writing #171. (I'm always trying to analyze the causes of these stupid mistakes I make. In this particular case I was in a rush, trying to finish some work before meeting my wife to go to the gym.)`

WebHubTel wrote:

It would be a good forecasting measure if we could reliably use it to make good forecasts! For forecasting, that's all that counts.

If the bumps in the Niño 3.4 index always looked sort of like the blue curve here:

perhaps we could reliably predict that when its first derivative hits a maximum

now, the index itself would hit a maximum in 6 months (or something like that: adjust the details to taste).But in fact bumps in the Niño 3.4 index are far more irregular in nature, so I doubt that simple methods like this work very well... at least, not for forecasting 6 months in advance.

For the very short term, probably linear extrapolation works reasonably well.

The good thing is that we can try lots of methods for predicting the Niño 3.4 index - silly-sounding methods and profound-sounding methods - and use statistics to quantify how well they work.

As a professional pedant, I make my living answering pedantic questions! The question you raise here is related to the question of whether certain supposed faster-than-light communication schemes actually work, in situations where not only the phase velocity but the group velocity of light exceeds $c$.

(I'm too lazy to explain exactly how these questions are related, but maybe you can guess.)

`WebHubTel wrote: > But eventually someone found out that B was simply the time derivative of A, or the rate of change of A. Would we consider that a good forecasting measure? It would be a good forecasting measure if we could reliably use it to make good forecasts! For forecasting, that's all that counts. If the bumps in the Niño 3.4 index always looked sort of like the blue curve here: <img src = "http://imageshack.com/a/img537/4002/XFQW8l.gif" alt = ""/> perhaps we could reliably predict that when its first derivative hits a maximum _now_, the index itself would hit a maximum in 6 months (or something like that: adjust the details to taste). But in fact bumps in the Niño 3.4 index are far more irregular in nature, so I doubt that simple methods like this work very well... at least, not for forecasting 6 months in advance. For the very short term, probably linear extrapolation works reasonably well. The good thing is that we can try lots of methods for predicting the Niño 3.4 index - silly-sounding methods and profound-sounding methods - and use statistics to quantify how well they work. > What do you say in response to someone that brings up such a pedantic question? As a professional pedant, I make my living answering pedantic questions! The question you raise here is related to the question of whether certain supposed faster-than-light communication schemes actually work, in situations where not only the phase velocity but the group velocity of light exceeds $c$. (I'm too lazy to explain exactly how these questions are related, but maybe you can guess.)`

By the way, Blake is computing daily average link strengths. His laptop may be slower than mine; it took him 12 hours to compute 10,000 of 24,000 pts.

`By the way, Blake is computing daily average link strengths. His laptop may be slower than mine; it took him 12 hours to compute 10,000 of 24,000 pts.`

Continuing my talk draft:

22) What can we conclude from this?

Using a linear model, the average link strength on a given month accounts for only 8% of the variance of Niño 3.4 index 6 months in the future. This doesn't sound very good, and indeed it's not.

However, there are more interesting things to say than this!

Both the Niño 3.4 index and the average link strength can be computed from the surface air temperature of the Pacific. The Niño 3.4 index explains 16% of its own variance 6 months into the future; the average link strength explains 8%, and taken together they explain 22%. So, they contain a fair amount of

independentinformation about the future Niño 3.4 index - in fact, a surprisingly large amount for just 2 variables, given how difficult the Niño 3.4 index is to predict.For comparison, Mahler used a random forest variant called ExtraTreesRegressor to predict the Niño 3.4 index 6 months into the future from much larger collections of data. Out of the 778 months available he trained the algorithm on the first 400 and tested it on the remaining 378.

The result: using a full world-wide grid of surface air temperature values explains only 23% of the Niño 3.4 index 6 months into the future. A full grid of surface air pressure values does considerably better, but explains only 34% of the variance. Using

twelve monthsof the full grid of pressure values only gets around 37%.(Are these last three percentages $R^2$'s???)This makes it seem more respectable that taken together, the Niño 3.4 index and average link strength at a given month explain 22% of the variance of Niño 3.4 six months in the future.

`Continuing my talk draft: 22) What can we conclude from this? Using a linear model, the average link strength on a given month accounts for only 8% of the variance of Niño 3.4 index 6 months in the future. This doesn't sound very good, and indeed it's not. However, there are more interesting things to say than this! Both the Niño 3.4 index and the average link strength can be computed from the surface air temperature of the Pacific. The Niño 3.4 index explains 16% of its own variance 6 months into the future; the average link strength explains 8%, and taken together they explain 22%. So, they contain a fair amount of _independent_ information about the future Niño 3.4 index - in fact, a surprisingly large amount for just 2 variables, given how difficult the Niño 3.4 index is to predict. For comparison, Mahler used a random forest variant called [ExtraTreesRegressor](http://scikit-learn.org/dev/modules/generated/sklearn.ensemble.ExtraTreesRegressor.html) to predict the Niño 3.4 index 6 months into the future from much larger collections of data. Out of the 778 months available he trained the algorithm on the first 400 and tested it on the remaining 378. The result: using a [full world-wide grid of surface air temperature values](https://9d8c7e6f260445992e26252c839fa61414632ec6-www.googledrive.com/host/0B4cyIPgV_Vxrb2wxUnFteXVwWHM) explains only 23% of the Niño 3.4 index 6 months into the future. A [full grid of surface air pressure values](https://2bd877097dcd3604697f17379c3f2232c9730610-www.googledrive.com/host/0B4cyIPgV_VxrbDFJS1dTVWxhV28) does considerably better, but explains only 34% of the variance. Using _[twelve months](http://www.googledrive.com/host/0B4cyIPgV_VxrX0lxSUxHU2VLN28/window-pressure-anom-predict.html)_ of the full grid of pressure values only gets around 37%. **(Are these last three percentages $R^2$'s???)** This makes it seem more respectable that taken together, the Niño 3.4 index and average link strength at a given month explain 22% of the variance of Niño 3.4 six months in the future.`

Daniel: there is a question marked by ??? in my previous comment (#151) which I hope you can help me answer.

I'm wondering if the "scores" for ExtraTreeRegressor can be seen as $R^2$ values. For example, I see a "score" in line 281 of this notebook:

Does this mean you get $R^2 = 0.23$ for the 378 test months? It would be nice to know the $R^2$ for the predicted El Niño 3.4 versus the actual El Niño 3.4.

What does the figure of 0.928 mean?

`Daniel: there is a question marked by ??? in my previous comment (#151) which I hope you can help me answer. I'm wondering if the "scores" for ExtraTreeRegressor can be seen as $R^2$ values. For example, I see a "score" in line 281 of [this notebook](https://9d8c7e6f260445992e26252c839fa61414632ec6-www.googledrive.com/host/0B4cyIPgV_Vxrb2wxUnFteXVwWHM): > print model.score(X=train, y=nino34.iloc[0:400,-1]), model.score(air[400:788,:], nino34.iloc[400:788, -1]) > 0.928413218247 0.235866144462 Does this mean you get $R^2 = 0.23$ for the 378 test months? It would be nice to know the $R^2$ for the predicted El Niño 3.4 versus the actual El Niño 3.4. What does the figure of 0.928 mean?`

The ExtraTreesRegressor.score method returns the $R^2$.

The line

prints the scores for both the train sets and the test set. Knowing the train set is useful for diagnostic purposes. Tree ensemble models like ExtraTreesRegressor are expected to get very high score on the train set since the effectively memorize the train set similary to knn models. With linear models on the other hand you do not want to see a big difference between train and test score, since that indicates overfitting. If the train scores are not what you expect that is a sign to take a closer look at the model.

`The [ExtraTreesRegressor.score](http://scikitlearn.org/stable/modules/generated/sklearn.ensemble.ExtraTreesRegressor.html#sklearn.ensemble.ExtraTreesRegressor.score) method returns the $R^2$. The line print model.score(X=train, y=nino34.iloc[0:400,-1]), model.score(air[400:788,:], nino34.iloc[400:788, -1]) 0.928413218247 0.235866144462 prints the scores for both the train sets and the test set. Knowing the train set is useful for diagnostic purposes. Tree ensemble models like ExtraTreesRegressor are expected to get very high score on the train set since the effectively memorize the train set similary to knn models. With linear models on the other hand you do not want to see a big difference between train and test score, since that indicates overfitting. If the train scores are not what you expect that is a sign to take a closer look at the model.`

Also I had a closer look at the

`summary.lm`

method in R and it seems we should quote the Multiple R-squared values.`Also I had a closer look at the `summary.lm` method in R and it seems we should quote the Multiple R-squared values.`

I am only remotely looking at what you do here, that is I am scanning your comments here merely for seeing whether there is a comment to my comment at #163 about the strange temperatures graphics. so my question might eventually sound displaced, but I found this sentence strange:

Are the SSTs the only parameters which enter Nino3.4? How is Nino3.4 computed from the temperatures?

`I am only remotely looking at what you do here, that is I am scanning your comments here merely for seeing whether there is a comment to my comment at <a href="http://forum.azimuthproject.org/discussion/1523/crunch-time/?Focus=13710#Comment_13710">#163</a> about the strange temperatures graphics. so my question might eventually sound displaced, but I found this sentence strange: >Both the Niño 3.4 index and the average link strength can be computed from the surface air temperature of the Pacific. Are the SSTs the only parameters which enter Nino3.4? How is Nino3.4 computed from the temperatures?`

As an average.

`> " How is Nino3.4 computed from the temperatures? " As an average.`

Thanks for clarifying the 'score', Daniel.

Good - I read about $R^2$ on Wikipedia and elsewhere and came to the same conclusion myself, mainly because the "multiple $R^2$" will be easier for me to explain if needed. In the talk draft I've already changed all the $R^2$ values from "corrected $R^2$" to "multiple $R^2$".

Luckily this changes the numbers very little, which suggests that nothing very sneaky is going on in our examples!

`Thanks for clarifying the 'score', Daniel. > Also I had a closer look at the summary.lm method in R and it seems we should quote the Multiple R-squared values. Good - I read about $R^2$ on Wikipedia and elsewhere and came to the same conclusion myself, mainly because the "multiple $R^2$" will be easier for me to explain if needed. In the talk draft I've already changed all the $R^2$ values from "corrected $R^2$" to "multiple $R^2$". Luckily this changes the numbers very little, which suggests that nothing very sneaky is going on in our examples!`