Climate Audit

pete

that’s not what I said … SO please take a little trouble before throwing stones.

It’s a direct quote of what you said. If you get in a snit whenever you’re corrected you’ll just turn this place into an echo chamber.

Steve: you also stated:

It terms of statistics, you’re wrong here. For a given level of correlation, significance will decrease with sample size.

as though that was news to me and as tho it corrected the prior point.,Your point is correct but not at issue in the quote. As I observed above, the benchmarks are different in the comparison depending on the length so your criticism is invalid.

pete

As I observed above, the benchmarks are different in the comparison depending on the length so your criticism is invalid.

The benchmarks over a shorter time period have to be harder (for a given significance), because the smaller sample size makes false positives more likely.

This means that a proxy that passes the (easier) benchmark over the longer time period might still fail the (harder) benchmark over the shorter time period.

So while 13% would be expected to pass the easier screening by chance, some smaller number (say x%) would be expected pass the harder screening you recommend. Therefore it makes sense to compare 40% (484/1209) to 13%, or to compare 28% (342/1209) to x%. Comparing 28% to 13% is a mistake.

bernie

Pete:
If I understand the point, I am not sure that increasing the difference in the sample sizes of 1850 to 1949 compared to 1850 to 1995 would be particularly material to the point that Steve was making, even if technically you are correct.

• pete
  Posted Jan 21, 2009 at 8:22 AM | Permalink

  Re: bernie (#6), if you decrease the sample size from 1850–1995 to 1850–1949 you decrease the number of significant proxies from 484 to 466.

  If you increase the number of sceening tests from one (1850–1995) to three (1850–1995, 1850–1949, and 1896–1995) you decrease the number of significant proxies from 484 to 342. (At one point Steve also adds a fourth test, screening for consistent sign).

  So you’re right, the reduced sample size isn’t the most important factor here.

  There are really two points to be made:
  (1) additional screening reduces false positives but increases false negatives
  (2) reducing the sample size increases false negatives without decreasing false positives

Fred

The Church of Climate Scientology marches on.

Obviously Statistics is not in the curriculum.

Peter Thompson

I apologize for asking a probably stupid question, but if a proxy’s correlation to temperature is deemed significant in one screening period, but not another, then what exactly can it tell you about past temperatures? To a layman, the obvious answer is “not much” insofar as has been proven to be a proxy for temperature some of the time. When in the past? When not? Shouldn’t a proxy correlate with the instrumental record at any point for which you have data, all other things being equal?

• pete
  Posted Jan 21, 2009 at 8:48 AM | Permalink

  Re: Peter Thompson (#10)

  if a proxy’s correlation to temperature is deemed significant in one screening period, but not another, then what exactly can it tell you about past temperatures?

  The proxy should be correlated in both periods, but this does not guarantee it will be significant in both periods.

  If a proxy is significant, then we know it is correlated (to some specified level of confidence). If a proxy is not significant, then we don’t know whether it’s correlated or not.

  Significant is a tricky piece of stats jargon; it has more to do with our knowledge than about the underlying process. And the less data we use, the less we know — which is why I recommend taking the correlation over the full period.

pete

Only 143 pass the above simple test ( 15.4%). This is only slightly above Mann’s 13%. Pick two-daily keno raises the bar a bit as well. One is forced to conclude that these dendro proxies as a group are not significant as well.

This is particularly wrong: you’re comparing the triple-screening result to the single-screening threshold — apples and oranges.

Mark T.

Statistically speaking, that’s an issue of stationarity, bernie, and it is assumed from the beginning in all of the reconstruction work, without ever being demonstrated to be true.

Mark

John Archer

OT (with apologies):

I’m trying to locate a graph I’ve seen (or rather a link to it) of delta T versus CO2 concentration. It’s the essentially logarithmic one.

I’d be very grateful if you could help. (My email is johnwdarcher=at=hotmail=dot=co=dot=uk.) Thanks.

Jeff Id

I do enjoy the Mann posts. Let’s see. First I see comments about tighter screening. This sounds good but it is IMHO not correct, the spurious correlations are too high? With all the infilling and autocorrelation a claim of X percent pass so this must be temp is not a fair test. Removing the infilled data I was able to get 39% to pass a negative linear slope screen. Had to mention it again. Still, it is kind of the crux of the paper so I understand why people keep going back to that.

The Briffa proxies are interesting because they pass to such a high percentage and yet have substantial infilling. It might be interesting to look at the degrees of freedom of the infilled portion of the data after RegEM compared to the normal set. Oops, I have a meeting and have to go.

Clark

So am I correct in reading what you wrote, for the MXD series they:

1. Remove the last 50-odd years from the series because the proxy gives the ‘wrong’ result (um, I mean there is ‘divergence’).

2. They replace this deleted data with made-up data taken from the proxies that pass the correlation-with-temperature screening.

3. After this, then they test the newly ‘improved’ series in a correlation test.

And that finally, the MXD series and the Luterbacher temperature-as-a-proxy-for-temperature series are the only group to pass better than random correlation?

Steve McIntyre

#1. the salient point in my inline comment could have been expressed without the extra editorializing. It was very late when I wrote this and there was no need to be Gavinesque. Sorry bout that.

#5,8. pete and DC, you’re both fixing on a point that was passim to the post. Mann asserted that a “passing” yield of 484 out of 1209 proxies showed something as only ~13% would pass by chance. The observation that I was attempting to highlight was the stratification of the yield: 100% of the Luterbacher instrumental series “passed”. No wonder – they are instrumental series. How can these be included in a statistical analysis of “proxies” in a serious article?

The second point was the very high yield (88%) in the M08 infilled version of Rutherford Mann ete al ReGEM gridded versions of Briffa MXD series, where the yield is far higher than run-of-mill dendro proxies.

The third point is the very low yield of low-freq proxies and dendro proxies.

I think that I may have misinterpreted one of Pete’s points above and my reply, while valid in one perspective, may have missed a more salient issue. The yield of red noise series that are “significant” in all 3 comparisons will necessarily be smaller than the yield for any one of the comparisons. What is the “right” benchmark allowing for pick two daily keno and autocorrelation? Dunno. Mann’s “test” is typically home-made, in the sense that there is no statistical reference for his test. So it’s very hard to tie his stuff to literature known off the Island. Is 13% the “right” number allowing both for 3 comparisons and pick two -daily keno? Dunno. Might be higher, might be lower.

What is the “right” sort of red noise from which to construct a test? Autocorrelation matters a lot – see Santer et al 2008.

And if (say) 19% of the dendro proxies passed the Mann test, what exactly does that mean? Does it mean that the series are a bit more autocorrelated than the red noise template? Or does it mean that they have a decodable signal?

At this point, I really don’t want to place undue weight on benchmarks, as these take time to analyse. My issue is primarily the simple one – Luterbacher and Briffa-Rutherford-Mann RegEM MXD.

In light of these comments (which I take as review comments), I’ve edited the post slightly to remove some passim references to the questionable Mannian benchmarks (which I will return to on another occasion).

• Mark T.
  Posted Jan 21, 2009 at 2:35 PM | Permalink

  Re: Steve McIntyre (#24),

  Or does it mean that they have a decodable signal?

  A colleague sent me a list of signal processing lemmas today, one of which seems entirely relevant.

  “In the limit, given enough matched filtering, enough cross-correlation, and enough coherent integration, there is no need for a signal.”

  There are two others that don’t really apply.

  Mark

• pete
  Posted Jan 21, 2009 at 6:22 PM | Permalink

  Re: Steve McIntyre (#24)

  #5,8. pete and DC, you’re both fixing on a point that was passim to the post.

  The point I’m fixing on is the very first substantive paragraph in this post; the passim points that you have removed are simply a follow-on from this original error. It’s hard to get to the intended point when you begin by confusing statistical and physical significance.

  In the spirit of offering a “review” comment:

  If your main point is the stratification, you would be better off ignoring the sub-period tests altogether and simply compare the strata using the original test. The way the post is laid out at present gives the impression that the stratification only shows up when using the triple-screening process, which would make the statistical mistakes in the triple-screening process relevant. But Luterbacher, for example, is still going to give 71/71 with the original test.

  • Deep Climate
    Posted Jan 21, 2009 at 11:33 PM | Permalink

    Re: pete (#49),

    Steve said:

    Code 9000 dendro proxies make up over 927 of 1209 M08 proxies. Only 143 pass the above simple test ( 15.4%).

    Then Pete said:

    If your main point is the stratification, you would be better off ignoring the sub-period tests altogether and simply compare the strata using the original test.

    For what it’s worth, 243 of the code 9000 candidate proxies pass under the original test (25%).

    I’m not sure which codes or other criteria correspond to “low frequency” (said to pass 7/51 under the more restrictive statistical screening), so I can’t review that particular statement at present.

    And, by the way, it’s a bit of a stretch to state that I “remarked in effect” that “30% of the proxies failed” what you term an “elementary precaution.” Rather, that is a conclusion you have drawn from data I drew to your attention. I’ve made it very clear that I don’t consider the sub-period significance test appropriate, although a much less stringent test for correlation fluctuation might be.

    Indeed my original point has been lost, namely, that only one of the 342 proxies that passed all three tests had an opposing sign in the sub-period correlations, which ran counter to your speculation:

    Out of 1209 proxies, 308 have opposite “low frequency” orientations between early-miss and late-miss. Presumably many are “significant” in both directions.

Steve McIntyre

#15. I totally reject the idea of ex post picking – which all too often are connected to things like “false positives”.

As I’ve said on many occasions, if, say, white spruce treeline ring width chronologies are believed to be a temperature proxy, then climate scientists have take them all. They can’t say – after the fact – that ones passing correlation tests are “proxies” and ones failing aren’t. If the first study missed a salient factor, then they need to do a fresh out-of-sample study using new ex ante selection to avoid data snooping bias.

EVen worse is mining the entire ITRDB data set for correlations with complete disregard for whether the series meet any ex ante criteria of being a temperature proxy.

You can get high RE statistics easily – Yule’s classic spurious regression between alcoholism and C of E marriages has a high RE. If you have “Method Wrong” (in Wegman’s phrase), the work isn’t saved by a high RE.

DC seems to take a different position.

JaneHM

Steve

More voodoo to be published in Nature tomorrow with new adjustments for Antarctica 1957 – 2006:

“The majority of weather stations on Antarctica sit around the coast, with only two providing an unbroken record from the continent’s interior. Steig and colleagues overcame this lack of data by using satellite data and statistical techniques to fill in the gaps…
“Eric has done a very clever analysis with extremely sparse data,””

http://www.newscientist.com/article/dn16460-even-antarctica-is-now-feeling-the-heat-of-climate-change.html

http://www.nature.com/nature/journal/v457/n7228/abs/nature07669.html

hmm

STEVE-
I’m confused on one point: Why doesn’t the “divergence” issue itself dominate the debate? Doesn’t divergence directly prove that these “proxies” are not as good as claimed? How can proxies not be calibrated over the entire instrumental period; what justification is there for this? Did they cherry pick the best correlated timeframe so that the computated errors would appear lower than we are observing in real life?

Shouldn’t the highest level of observed divergence in the instrumental period be added into the error computation in the reconstruction period (and if so, was it)? If it can diverge this much now, it could in the past too IMO.

I would suggest that instead of deconstructing their studies, you create your own study starting with the raw data available. You use your own methods and get your own results and compare to theirs.

Edward

#32 The Ford Prefect
You state: “Proxies are aligned with temperature records”. Is that an assumption? Based on my reading of this blog that statement may or may not be true and certainly depends on the particular proxy you would like to discuss.

Mark T.

Yes, and what we were trying to understand was a question of relevance to the issue at hand, which I think has been detailed sufficiently since: not so much.

Mark

Sam Urbinto

Hey, where’s the infrared photos?

The ones in places with higher concentrations of carbon dioxide are especially cute.

William

#32 Ford Prefect
When I veiwed your image it teleconnected for me and brought up an image of an old I love Lucy re-run. Would that be considered a high correlation to temperature?

Alan Wilkinson

Using proxies that are not understood to the ridiculous extent that their relevance is not predictable prior to use is never going to yield reliable estimates of temperature.

When they are anyway non-linear to the point of sign reversal they become farcical.

Yet again this is not science, it is politics.

Steve McIntyre

For those following along at home: the other 307 with opposite orientations were non-significant proxies.

Don’t think so. Of these 308 proxies, 151 are considered “significant” in at least one of the 3 versions.

temp1= !(sign(details$rtable.r1850_1949lf)== sign(details$rtable.r1896_1995lf)) ;sum(temp1) #308
temp2=temp1& (passing$whole|passing$latem|passing$earlym);sum(temp2) #151

The usual interpretation of such failures in reconstruction literature is that the posited relationship was not “verified”. I agree that the Socotra case seems to be especially bad by yielding both a “significant” positive and “significant” negative relationship, whereas the others appear to have (say) a “significant” positive relationship in one period and an “insignificant” negative relationship in another period – not particularly edifying.

Statistically, it’s useful to ponder how you can have both a”Significant” positive and “significant” relationship. The explanation is almost certainly that the data is highly autocorrelated, that the degrees of freedom are much fewer than assumed by Mann (the type of argument used by Santer against Douglass) and thus neither relationship is actually “significant” allowing for a very modest number of degrees of freedom. I’ll do these calcs tomorow. The knock-on impact of applying Santer-style methods will probably be to eliminate quite a few “significant” relationships, but that’s just a guess.

I remind readers that Mann’s benchmarking lacks any statistical reference and cannot be relied on. I haven’t parsed whether this should be 13% or some other number or considered what an 18% yield actually means (given that M08 is very obscure on this.) As I noted above, the yield represented by “484” is highly inflated by inclusion of Luterbacher and perhaps by RegEM’ed MXD data. The percentage of “significant” dendro proxies is remarkably low even with Mannian under-allowance for autocorrelation.

• pete
  Posted Jan 22, 2009 at 4:53 AM | Permalink

  Re: Steve McIntyre (#62)

  I’d prefer to do it this way:

  temp1 <- sign(details$rtable.r1850_1949lf) != sign(details$rtable.r1896_1995lf) ; sum(temp1)
  # 308
  temp2 <- temp1 & (passing$whole & passing$latem & passing$earlym) ; sum(temp2)
  # 40

  You can’t determine if the correlations have opposite sign unless you can determine the sign, which you can’t do unless your estimated correllation is significant.

  Statistically, it’s useful to ponder how you can have both a”Significant” positive and “significant” relationship. The explanation is almost certainly that the data is highly autocorrelated,

  Actually, this is about as many as I’d expect by chance:

  Of the 484, say 376 are genuine and 108 spurious [108 / (1209 – 376) ~ 13%]
  Of those 108, 76 survive the two sub-period tests [108 * (342/484) ~ 76]
  Of those 76, 38 should have opposite signs [76/2 = 38; it’s a coin flip since they could have the same sign by chance]

  Which is pretty close to the 40 I got earlier.

  The percentage of “significant” dendro proxies is remarkably low even with Mannian under-allowance for autocorrelation.

  28% (260/927) compared to a threshold of 13%. That’s not “remarkably low”.

  Steve: As a start, 484 is a fake number because it includes 71 out of 71 Luterbachers. Plus the 88% yield of M08 massaged Rutherford-RegEMd Briffa MXD data is suspect. Whatever you do needs to be done without these datasets [codes 2000, 7500]

  • RomanM
    Posted Jan 22, 2009 at 2:59 PM | Permalink

    Re: pete (#65),

    Actually, this is about as many as I’d expect by chance:

    Of the 484, say 376 are genuine and 108 spurious [108 / (1209 – 376) ~ 13%]
    Of those 108, 76 survive the two sub-period tests [108 * (342/484) ~ 76]
    Of those 76, 38 should have opposite signs [76/2 = 38; it’s a coin flip since they could have the same sign by chance]

    Which is pretty close to the 40 I got earlier.

    Actually there are several things wrong with these calculations. As I pointed out in Comment 54 of the Mann Correlation Mystery thread, 13% is not correct. If you read the screening procedure on page two of the SI,

    The corresponding one-sided p=0.10 significance thresholds are |r|=0.11 and |r|=0.34 respectively. For the shorter (100 year) calibration intervals, we assumed n=98 nominal degrees of freedom for annually resolved records, and n=8 degrees of freedom for decadal resolution records. The corresponding one-sided p=0.10 significance thresholds are |r|=0.13 and |r|=0.42 respectively. Owing to reduced degrees of freedom arising from modest temporal autocorrelation, the effective p value for annual screening is slightly higher (p~0.128) than the nominal (p=0.10) value.

    you will notice that they themselves point out that the critical value is for a one-sided test at the .10 level (which I verified in my linked comment). However, as soon as they take the absolute value of the correlation coefficient before making the comparison, they convert it into a two-sided test with a corresponding doubling of the significance level. The nominal level is p = 0.20 and with their “slightly higher” correction gives an error rate of 25.6% not 13%.

    I would also take exception with your statement that “it’s a coin flip since they could have the same sign by chance”. In fact, because the two correlations overlap for 54 years of each validation period, I would expect a pretty strong positive correlation between the two correlations thereby increasing the proportion of pairs with the same signs.

    I wrote a little script to simulate correlated two independent normal series in the same fashion (since I didn’t have access to “climate model data” 😉 ):

    #n = number of repetitions
    samp = function(n) {
    result = matrix(NA,nrow=n, ncol=2)
    for (i in 1:n) {
    rand1 = rnorm(146)
    rand2 = rnorm(146)
    result[i,1] = cor(rand1[1:100],rand2[1:100])
    result[i,2] = cor(rand1[47:146],rand2[47:146])}
    result}

    #answers will vary
    test=samp(10000)
    cor(test[,1],test[,2]) # 0.5307599

    #number with same signs
    sum(test[,1]*test[,2] >0) # 6789

    In several runs, the correlation of the two calculated correlations was between .50 and .55 and the percentage of cases where both of the signs were the same was about 68%.

    • pete
      Posted Jan 22, 2009 at 6:10 PM | Permalink

      Re: RomanM (#69)

      you will notice that they themselves point out that the critical value is for a one-sided test at the .10 level (which I verified in my linked comment). However, as soon as they take the absolute value of the correlation coefficient before making the comparison, they convert it into a two-sided test with a corresponding doubling of the significance level. The nominal level is p = 0.20 and with their “slightly higher” correction gives an error rate of 25.6% not 13%.

      I read this as converting a two-sided p=0.10 test with thresholds +/-0.13 into a one-sided p=0.10 test with threshold +0.13.

      Thanks for the correction re “coin flip”, I was a little careless there.

      Um, Steve? Could you perhaps delete my post above? That huge chunk of whitespace was meant to be a +/- sign; might have to brush up on my TeX.

    • RomanM
      Posted Jan 22, 2009 at 6:53 PM | Permalink

      Re: pete (#78),
      All you have to do is calculate the probability that that a correlation (without absolute value exceeds the critical value. here is an R program:

      #Calculate probability correlation greater than r

      probcor = function(r,df) {
      tval = r/sqrt((1-r^2)/df)
      prob = pt(tval, df,lower.tail = F)
      list(t=tval, p =prob)}

      df = c(8, 13, 98, 144)
      r =c(.42, .34, .13, .11)

      probcor(r,df)$p # 0.11344236 0.10750085 0.09867576 0.09312986

      Since the distribution is symmetric around 0, the absolute value doubles the probability. The differences from .10 come from rounding the critical value to two digits (why would anyone round the values when computers are being used???). My suspicion is that someone got the idea that since the comparison on the negative side was unnecessary due to the absolute value, the test was therefore “one-sided” – an error that I’ve seen before in elementary courses…

      Steve: Roman, there’s probably a native function that does this even more directly. I can directly reverse this with qnorm:

      qnorm( 1-c (0.11344236, 0.10750085, 0.09867576, 0.09312986 ),sd=1/sqrt(c(8,13,98,144))
      #0.4272419 0.3438944 0.1302222 0.1101438

      A question: I’ve been taking the tanh of this expression for the Fisher transformation i.e.

      tanh(qnorm( 1-c (0.11344236, 0.10750085, 0.09867576, 0.09312986 ),sd=1/sqrt(c(8,13,98,144)))
      #0.4030138 0.3309497 0.1294911 0.1097005

      IT doesn’t make much difference but isn’t it required in some of these calcs?

    • pete
      Posted Jan 22, 2009 at 8:16 PM | Permalink

      Re: RomanM (#80), thanks got it now.

      (getting some weird errors, sorry if this double-posts)

    • Deep Climate
      Posted Jan 23, 2009 at 11:33 AM | Permalink

      Re: RomanM (#69),

      The corresponding one-sided p=0.10 significance thresholds are |r|=0.11 and |r|=0.34 respectively.

      You will notice that they themselves point out that the critical value is for a one-sided test at the .10 level (which I verified in my linked comment). However, as soon as they take the absolute value of the correlation coefficient before making the comparison, they convert it into a two-sided test with a corresponding doubling of the significance level.

      This appears to be a misinterpretation of the text. Mann et al are simply using a short hand notation to describe the two possible one-sided tests. But in the actual significance test for the vast majority of individual proxies, only one or the other is applied (otherwise it would be a two-sided test).

      At least, that’s how I understood it. But if you’re not convinced, check the data:

      temp1= (abs(details$rtable.r1850_1995) > .1) ;sum(temp1) # 930
      temp2 = temp1 & passing$whole ;sum(temp2) # 484

      It looks like there were quite a few proxies that would have passed a two-sided test, if it had been applied. But in the vast majority of cases it wasn’t.

    • RomanM
      Posted Jan 23, 2009 at 12:19 PM | Permalink

      Re: Deep Climate (#82),

      You are totally and completely confused. Have you not taken an elementary stat course?

      What we are talking about here is a SINGLE comparison involving the correlation of the proxy to a target sequence to try to determine whether the calculated correlation could be as far from zero as it is if there actually was no relationship between the two series.

      In this case, we have what is called a two-sided test because the decision is made to say that there is a relationship if the correlation is either large enough on the positive OR large enough on the negative side. If we were looking for positive correlation only, we would ignore the negative values (no matter how far below zero) and only do our check for large enoughness on the positive side. That is called a one-sided test.

      Now, the value that has been selected by the Mann will accept a non-existent, spurious correlation about 10% of the time when a one-sided test is done . That is what the SI correctly states (which my calculation #79 verifies – if you don’t agree with the calculations, by all means, point out the errors). Now in an elementary stat course, you would learn that using that same cutoff value on both sides of zero (for a symmetrically distributed statistic) will double the error rate. By taking the absolute value of that statistic, this is EXACTLY what you are doing even though it might look like you are only comparing on one-side, a mistake made by the folks who wrote this paper and selected the WRONG cutoff value for their stated error rate. Got it?

      There will be a test on this on Monday and you really need to pull up your grades if you expect to pass. 😉

    • Deep Climate
      Posted Jan 23, 2009 at 2:09 PM | Permalink

      Re: RomanM (#83),

      From Mann et al 2008:

      Where the sign of the correlation could a priori be specified (positive for tree-ring data, ice-core oxygen isotopes, lake sediments, and historical documents, and negative for coral oxygen-isotope records), a one-sided significance criterion was used. Otherwise, a two-sided significance criterion [e.g. spelo] was used.

      [Emphasis added]

      Here is the result of applying the two-sided test to every proxy: 902 pass.

      temp1= (abs(details$rtable.r1850_1995) > .106) ;sum(temp1) # 902

      Here’s the number that actually passed: 484

      temp2 <- temp1 & passing$whole; sum(temp2) # 484

      Honestly, I don’t know what else to say.

    • RomanM
      Posted Jan 23, 2009 at 2:29 PM | Permalink

      Re: Deep Climate (#85),

      So show me where in Mann’s program where he actually did a one-sided test.

    • Deep Climate
      Posted Jan 23, 2009 at 6:00 PM | Permalink

      Re: RomanM (#86),

      So show me where in Mann’s program where he actually did a one-sided test.

      Here are some code snippets from gridproxy.m

      One sided positive significance test for most of the proxies (tree-rings etc.):

      if (z(3,i)==9000 | z(3,i)==8000 | z(3,i)==7500 | z(3,i)==4000 | z(3,i)==3000 | z(3,i)==2000) &…
      x(kk,i+1)>-99999 & x(kkk,i+1)>-99999 &…
      z(1,i)>=ilon1 & z(1,i)=ilat1 & z(2,i)<=ilat2 & z(ia,i)>=corra

      One-sided negative significance test for code 7000 proxies:

      if ( z(3,i)==7000) &…
      x(kk,i+1)>-99999 & x(kkk,i+1)>-99999 &…
      z(1,i)>=ilon1 & z(1,i)=ilat1 & z(2,i)<=ilat2 & z(ia,i)<=corra2

      Two-sided test for code 5000/6000 proxies only (relatively small number of proxies)

      if (z(3,i)==6000 | z(3,i)==5000) &…
      x(kk,i+1)>-99999 & x(kkk,i+1)>-99999 &…
      z(1,i)>=ilon1 & z(1,i)=ilat1 & z(2,i)<=ilat2 & abs(z(ia,i))>=corra3

      These are for the annually resolved proxies, but I see similar code for the decadally resolved proxies too. It appears that Mann did what he said he did.

    • RomanM
      Posted Jan 24, 2009 at 10:14 AM | Permalink

      Re: Deep Climate (#88),

      Score one for DC and the hockey team. You are corect, they do one-sided tests for the proxies in several of the code snippets that you list. My bad! I had not gone through all of the Matlab programs in detail and overlooked the pieces that you listed. I agree that Prof. Mann et al. are aware of the difference between one-sided and two-sided tests and apologize to them for having unfairly implied the opposite. 😦

      Now having eaten some crow, I have to ask why in a serious scientific study, anybody would do a “two-sided test for code 5000/6000 proxies only (relatively small number of proxies)”. If you don’t know in advance whether the correlation is positive or negative (indicating that the physical processes behind that link are poorly and/or incompletely identified), how could you assume that that relationship is specifically to temperature in a linear fashion? Black box teleconnection at its finest!

• Deep Climate
  Posted Jan 23, 2009 at 12:53 PM | Permalink

  Re: Steve McIntyre (#62),

  Of these 308 proxies, 151 are considered “significant” in at least one of the 3 versions.

  I think it’s worth having a little more precision on this issue, focusing only on the “significant” sub-period proxies.

  First of all in the screening process (using HF):

  temp1 <- sign(details$rtable.r1850_1949) != sign(details$rtable.r1896_1995) ; sum(temp1) # 229
  temp2 <- temp1 & (passing$latem & passing$earlym) ; sum(temp2) # 1
  temp3 <- temp1 & (passing$latem | passing$earlym) ; sum(temp3) # 102

  So out of 229 opposite signed proxies, only one was “significant” and passed the screen. Among proxies that passed in one sub-period only, 101 were opposite-signed in the other subperiod – all rejected.

  Only 1/229 opposite-signed proxies was significant.

  Now let’s look after conversion to low frequency for calibration. Here there is no further test for significance. So we would expect some “passing” proxies to change, even to the point of changing signs.

  temp1 <- sign(details$rtable.r1850_1949lf) != sign(details$rtable.r1896_1995lf) ; sum(temp1) # 308
  temp2 <- temp1 & (passing$latem & passing$earlym) ; sum(temp2) # 40
  temp3 <- temp1 & (passing$latem | passing$earlym) ; sum(temp3) # 136

  So the correct number, as pete first showed, is 40/308, not 151.

  Based on a complete understanding of the process, I don’t find any of these numbers particularly surprising or problematic.

Steve McIntyre

Mann stuff is always needles in eyes. Consider the following:

We assumed n = 144 nominal degrees of freedom over the 1850–
1995 (146-year) interval for correlations between annually resolved
records (although see discussion below about influence of
temporal autocorrelation at the annual time scales), and n = 13 degrees
of freedom for decadal resolution records. The corresponding
one-sided P = 0.10 significance thresholds are |r| = 0.11 and
|r|~ 0.34, respectively.

They have 8 annual series with stated correlations less than r=0.11. It looks like they used something lower – about .106.

More annoying is this. They have 3 series that don’t meet the low-frequency guidlines:

id code r1850_1995 rtable.r1850_1995lf
483 li_1998_d13c 6001 NA -0.3907627
1043 suk_1987_koreadrought 5001 NA -0.3540128
1205 zhang_1980_region3 5001 NA -0.3541390

Here’s it’s not that they used a different number, as there are a few series in the same range that “pass”. I guess maybe they forgot to flip in time.

# id code r1850_1995 rtable.r1850_1995lf
#316 chesapeak_paleosalinity 3001 0.3886 0.3886366
#382 curtis_1996_d13cpyro 4001 0.3971 0.3970749
#1055 thompson_2003_dasuopu 8001 0.3423 0.3422881

• Deep Climate
  Posted Jan 22, 2009 at 5:53 PM | Permalink

  Re: Steve McIntyre (#72),

  Hmmm …

  These two end in 1960 – could that reduce d.f. just enough to fail?

  1043 suk_1987_koreadrought 5001 NA -0.3540128
  1205 zhang_1980_region3 5001 NA -0.3541390

  Can’t explain the other one though.

Steve McIntyre

#74. Doubt it. Remember that Mann infills everything for these calcs, so these series have data from 1960-90 (hence the correlations) even if the data’s fictitious.

Peet

This is particularly wrong: you’re comparing the triple-screening result to the single-screening threshold — apples and oranges.

pete

So when Mann says …

The corresponding one-sided p=0.10 significance thresholds are |r|=0.11 and |r|=0.34 respectively.

… he’s talking about two different one-tailed thresholds: r > 0.11 for a priori positive correlations and r < -0.11 for negative correlations?

That would explain #72

• Deep Climate
  Posted Jan 24, 2009 at 10:07 AM | Permalink

  Re: pete (#89),

  So when Mann says …

  The corresponding one-sided p=0.10 significance thresholds are |r|=0.11 and |r|=0.34 respectively.

  … he’s talking about two different one-tailed thresholds: r > 0.11 for a priori positive correlations and r < -0.11 for negative correlations?

  That’s exactly the point I’ve been trying patently to explain to the CA crowd; I think they get it now.
  [Steve: while I’ve taken exception to many aspects of Mann’s methods, this was not an issue that I had raised or taken exception to. Roman had – and has agreed that he misunderstood Mann’s methods (which is easy enough to do), but I’m not sure that anyone else had weighed in on the topic, so it’s a bit grandiose to say that you were dealing with a “crowd” on this topic.
  

  Of course, for low frequency one-tailed thresholds (|r|=0.34) the two tests would be r > 0.34 for a priori positive correlations and r < -0.34 for negative correlations.

  That would explain #72

  Sort of. I went back to that comment and noticed that the failing low-frequency proxies were code 5001 and 6001. So these would have been subjected to the two-sided test (as described in the paper and implemented in the code). I’d have to check the documentation and code, but obviously you would need a higher absolute threshold for the two-sided test to achieve the same level of significance. So it seems very plausible that that’s the explanation.

  As for the high frequency one-sided proxies that passed but were under .11, it’s clear that the actual calculated threshold was .106, which was rounded up in the SI documentation.

Hu McCulloch

Re #1, 2, Steve originally wrote,

If a “proxy” is a proxy, then it is a proxy regardless of the subperiod.

True enough — if a proxy’s relation to temperature isn’t constant across subperiods of the calibration period, there is no reason to expect it to have remained constant back in the reconstruction period, either.

However, as Pete points out (with unnecessary rudeness), Steve’s next sentence is too strong. Steve said,

It is not enough to have a “significant” relationship in the 1850-1995 period, it should also have “significant” relationship in the 1850-1949 and 1896-1995 periods (Mann’s late-miss and early-miss periods.)

Requiring the proxy to pass three hurdles in this manner — a “Leap 3” criterion, as it were — is just as size-distorting, but in the opposite direction, as Mann’s “Pick 2 Keno” criterion. In fact, as the subperiods get shorter, the power of a test at the level in question gets smaller, and so one is almost guaranteed not to find a significant correlation in subperiods if they are sufficiently short.

Although it is too much to ask that the data show a significantly nonzero relationship in each subperiod as well as the whole period, a very different but valid double check is whether the relationship is significantly non-constant between subperiods. This can be done with a standard Chow test, breaking the whole period in half.

The standard two-regime Chow test could be easily modified in the spirit of the Mann 08 paper to split the whole period into three approximately equal periods, and then do a (single) test to see if the coefficients are equal in all 3 periods. Howevever, the power of such a test would be reduced. My inclination would be to stick with 2 periods, but if someone wants to do 3 that’s fine too.

The difference in the two approaches is that Steve would require the series to reject-reject-reject non-zero slope in the full and 2 subsamples, whereas I would have it reject nonzero slope in the full sample, and then not reject unequal coefficients between periods. My suggestion is still a double test of sorts, but it is more of the nature of leap one hurdle and then don’t fall on your face, rather than leap two hurdles.

On the opposite distortion generated by Mann’s Pick-2 Keno procedure, newbies are strongly urged to go back and read the 9/20 thread The Mann Correlation Mystery, and then to browse through the remainder of the Mann 08 Category on CA.

• Deep Climate
  Posted Jan 24, 2009 at 10:37 AM | Permalink

  Re: Hu McCulloch (#90),

  This can be done with a standard Chow test, breaking the whole period in half.

  The standard two-regime Chow test could be easily modified in the spirit of the Mann 08 paper to split the whole period into three approximately equal periods, and then do a (single) test to see if the coefficients are equal in all 3 periods.

  Your idea is a better version of the idea I raised back in Deep Climate (#15),

  I do think some sort of check on correlation “fluctuation” might be appropriate, although requiring two sub-period correlations to also be significant would be too restrictive. One possible simple test would be a sign test, requring the sub-periods to be of the same sign as the overall correlation, but not necessarily to the level of significance. (The sub-periods should also not overlap for this particular purpose, of course – half periods would seem appropriate).

  One implication of any sub-period test, of course, is that the same test should be applied in the validation exercise to the calibration subperiod. So you divide the whole period in halves for the “whole” reconstruction Chow test. And you would do the same thing within the two sub-period validation reconstructions.

Hu McCulloch

Re #82, 83, 89, I think Pete and DC are right here — Mann just means that the 1-tailed critical value for r is +.11 or -.11, as the case may be, but then confusingly merged this into the single statement that the absolute value of the critical value is +.11, which could easily be taken as saying that the critical value for the absolute value of r is .11. As RomanM correctly objects, the latter would be a 2-tailed test, but Mann seems to be giving it the former interpretation. (A more exact value is .1067.)

Three valid problems with Mann’s critical value are:

a) Mann’s Pick2 procedure, which greatly disorts its validity, as Steve showed earlier on the 9/20 thread The Mann correlation Mystery,

b) The “one size fits all” adjustment for serial correlation claiming that the SC-adjusted test size is about 0.128 for the annually resolved proxies and 0.1 for the decadal proxies. I can see applying a single adjustment to each category of proxies, but surely tree rings widths, grape harvests, and ice cores have completely different serial correlation properties. An average serial correlation should have been calculated for each category, reported, and used to compute SC-adjusted critical values at a common p-value.

c) The “Voodoo” data mining procedure that applies a critical value for a single test to 1209 tests, and that was the original topic of this thread and its predecessor. More on that later.

Hu McCulloch

Back, now to Voodoo Correlations and Correlation Picking, which was the original topic here. (See also Industrial Strength Voodoo Correlations).

As I discussed at Comment #19 of the original Voodoo Correlations and Correlation Picking thread, the standard Bonferroni size adjustment states that the individual test size α* may have to be set as low as α/n in order for the probability of a false reject to be no more than α. This is just an inequality that covers the worst possible case, but is not much different, for small test sizes, than the Šidàk adjustment, which is exact in the case of independent tests.

With n = 1209 and a desired collective test size of α = .10, the adjusted individual test size could have to be as low as α* = .1/1209 = .0000827. The Šidàk value of 1-.91/n = .0000871 isn’t different enough to be worth calculating. For the high-frequency proxies, the Bonferroni-adjusted 1-tailed critical value becomes 0.3068 (or -.3068, as the case may be) rather than 0.1067 (or -.1067) as for a single test. The low-frequency adjusted 1-tailed hurdle is 0.8228 (or -.8228).

Scanning the SI ProxyNames spread sheet, I see at least a few TRW’s that are greater than .3068, plus a number of other proxies with even higher value. I take this to mean that only these proxies are clearly significant at the .1 level, and that only these should be included in any reconstruction period where all 1209 are available, after eliminating the Pick-2 Keno trick and doing a more detailed adjustment for serial correlation.

If only a handful of the 1209 are active for a given period (say the MWP), the adjusted critical values should be recomputed using the actual number of active proxies and the selection procedure redone with these adjusted critical values.

Since the tests are in fact positively correlated (in the sense of being positively dependent), the Bonferroni and even Šidàk criteria are admittedly probably too severe. A more powerful test that all the proxies are zero (and therefore of determining that at least 1 is non-zero) would be based on an F statistic that took the correlations of the proxy residuals on temperature (as well as their individual serial correlation) into account.

Unfortunately, the matrix of the unrestricted correlations of these residuals is singular, since there are far more proxies than observations. However some sort of test like this might still be computable under strong restrictions on the covariances, such as that the members of each class of proxies have a common correlation coefficient with temperature, and that different classes (TRW vs ice core, eg) have uncorrelated errors.

A further problem with an F test is that it is intrinsically based on a 2-tailed alternative hypothesis for the coefficients. Perhaps this can be gotten around by Monte Carlo methods, but it’s a big complication if the one-tailed priors are to be maintained.

UC

propick.m is very difficult to follow, modified the code to show some values,

if sum(d3)>=igr % at least ‘igr’ correlation of ‘igrid’ pass the CL
n4=n4+1;
aa
bb=abs(aa)
ra=max(bb)
[iii,jjj]=find(bb==ra)
if sum(iii)>2
rra(n4)=aa(1,1); %07/20/06
else
rra(n4)=aa(iii,jjj);
end

rra(n4), pause
d4(n4)=i+1;
end

%% two nearest correlations

aa =

-0.1852
0.1093

%% both significant ? 🙂

%% absolute value
bb =

0.1852
0.1093

%% max abs correlation
ra =

0.1852

%% selected rr(n4)
ans =

-0.1852

Where is rtable1209late generated, btw?

Steve McIntyre

As you observe, the code numbers are slightly different tho very close to this derivation.

However you can’t get .337 because the steps for N are too big and rounding won’t do it.

I presume that this must have done the calculation some Mannian way with Mannian rounding errors accumulating along the way.

Regardless, my derivation reconciles all the retentions and permits further analysis, which I’m doing.

• Deep Climate
  Posted Jan 24, 2009 at 3:18 PM | Permalink

  Re: Steve McIntyre (#102),
  Sure, it’s close enough for your purposes – no problem with that.

RomanM

Re:Steve’s inline comment on RomanM (#79)

The tanh transform that you are using is an approximation and is less accurate for small degrees of freedom. It is used mainly when you wish to calculate confidence intervals or test whether the correlation might be equal to some non-zero value.

For the test of zero correlation, the statistic

t = r / sqrt((1-r^2)/df)

where df = sample size – 2 is used. This statistic is distributed exactly as a t random variable with df degrees of freedom. To calculate the critical value, you use the inversion formula which I gave in the post. A simple R script for that calculation is

#calculate upper quantiles for testing corr = 0

quant = function(alpha, df) {
tval = qt(alpha,df, lower.tail=F)
(sqrt( tval^2 /(tval^2 +df)))}

df = c(8,13,98,144)

quant(.10,df)
#[1] 0.4427959 0.3506884 0.1292418 0.1066760

quant(.05,df)
#[1] 0.5493568 0.4408608 0.1654298 0.1366643

The corresponding values from gridproxy appear to be:

alpha = .10: .425 .337 .136 .106

alpha = .05: .521 .418 .164 .128

These differences could not be due to rounding as we know it.

• Deep Climate
  Posted Jan 24, 2009 at 3:55 PM | Permalink

  Re: RomanM (#104),

  df = c(8,13,98,144)

  quant(.10,df)
  #[1] 0.4427959 0.3506884 0.1292418 0.1066760

  quant(.05,df)
  #[1] 0.5493568 0.4408608 0.1654298 0.1366643

  The corresponding values from gridproxy appear to be:

  alpha = .10: .425 .337 .136 .106

  alpha = .05: .521 .418 .164 .128

  These differences could not be due to rounding as we know it.

  I think some of the values got transposed:

  alpha = .10: .418 .337 .128 .106

  alpha = .05: .521 .425 .164 .136

  For example, .136 is the threshold used for the “whole” period two-sided test (i.e. df=144 with .05 two sided).

  See Deep Climate (#101)

  for the four “whole” recon parameters (df=13, df=144).

  • RomanM
    Posted Jan 24, 2009 at 4:24 PM | Permalink

    Re: Deep Climate (#109),

    Makes sense and the smaller values are close, but it still isn’t rounding so who knows where the values came from.

    • Deep Climate
      Posted Jan 24, 2009 at 4:42 PM | Permalink

      Re: RomanM (#112),

      Steve has come very close for df=144, 13, so my hunch is that his method would come closer for the sub-period cases too.

    • RomanM
      Posted Jan 24, 2009 at 5:55 PM | Permalink

      Re: Deep Climate (#113),

      What I have calculated is correct – this material has been my bread and butter for 40 years. If the results differ by as much as they do, particularly for the lower degrees of freedom, then what was done by Mann et al is not the same process. What did they do?

    • Deep Climate
      Posted Jan 24, 2009 at 6:54 PM | Permalink

      Re: RomanM (#114),

      I’m not disputing the correctness of your calculations. I’m just saying Mann et al values are quite close to the approximation that Steve used (but still not exact).

      tanh( qnorm(.90, sd=1/sqrt(98))) # 0.1287379
      tanh( qnorm(.95, sd=1/sqrt(98))) # 0.1646430
      tanh( qnorm(.90, sd=1/sqrt(8))) # 0.4244413
      tanh( qnorm(.95, sd=1/sqrt(8))) # 0.5237864

      alpha = .10: .418 .337 .128 .106
      alpha = .05: .521 .425 .164 .136

      On the other hand, if you add 1 to each df, you get this from your program:

      df = c(9,14,99,145)
      quant(.1,df) # 0.4186622 0.3382816 0.1285896 0.1063083
      quant(.05,df) # 0.5214044 0.4259020 0.1645995 0.1361950

      That’s even closer – so who knows?

    • RomanM
      Posted Jan 24, 2009 at 8:11 PM | Permalink

      Re: Deep Climate (#115),

      That’s even closer – so who knows?

      So who knows what? There is a right way and there are wrong ways. You don’t just “add 1 to each df” and you don’t invent new ways to calculate “principal components” because there is no valid reason to do so. Are you implying that what they have calculated these values correctly? Just a little wrong? Or did they (as you may have discovered) use the incorrect value for the degrees of freedom?

    • Deep Climate
      Posted Jan 24, 2009 at 9:16 PM | Permalink

      Re: RomanM (#116),
      The question I asked was the same as yours: “who knows where the values came from?” It’s strictly a “reverse engineering” puzzle for me and I’m not implying that the method was correct, whatever it was.

      As I said before, as far as I know your threshold values are correct. And so they probably either used an approximation or made some other error that gave low threshold results, especially for the low d.f. cases. So I don’t really see that we have a disagreement here.

      Steve: I think that you’re being careful here, but just in case: Roman is simply trying to do things in an intelligible way; don’t assume that Mannian “benchmarks” calculated intelligibly are “correct”. Roman doesn’t, nor do I. There’s much still to disentangle.

    • RomanM
      Posted Jan 25, 2009 at 8:10 AM | Permalink

      Re: Deep Climate (#117),

      Given your previous history as a poster, I am leery of your motives (unquestioning vindication of the team and their methods – they simply can’t be wrong in any way) and interpret everything you say in that context. If I see any admissions from you that some of the publications under consideration might have inappropriate and/or incorrectly done analyses or unjustifiable conclusions, I might see your posts in a different light.

      Following up on your observation that the degrees of freedom might not be correctly used, here is some R output

      df = c(8,13,98,144)
      floor(1000*quant(.10,1+df))/1000 #[1] 0.418 0.338 0.128 0.106
      floor(1000*quant(.05,1+df))/1000 #[1] 0.521 0.425 0.164 0.136

      Only one of the values is still different (.338 is out by .001). To achieve this, all of the degrees of freedom are (incorrectly) increased by one and all of the results are rounded down to three decimal places (presumably as a “conservative” approach).

      Let’s call .338 a “typing error” by whomever typed up the matlab script. 🙂

    • RomanM
      Posted Jan 25, 2009 at 9:29 AM | Permalink

      Re: RomanM (#120),

      … wait a minute! Rounding down isn’t “conserative”. It allows for slightly easier acceptance as a proxy, not the opposite.

Steve McIntyre

#105. You’ve changed topics here. The estimation of uncertainties (which is a separate can of worms) has nothing to do whether there is a “significant” yield of correlated proxies and whether Mann has misrepresented the situation here.

Even if the uncertainties were done correctly (about which I’m dubious), it doesn’t justify the inclusion of Luterbacher instrumental data into the reporting of proxy yields – different topic.

• Deep Climate
  Posted Jan 24, 2009 at 4:12 PM | Permalink

  Re: Steve McIntyre (#106),

  OK, point taken. Another way to get at this issue would be to analyze the “pass” rate by proxy length. For example, what is the “pass” rate for proxies going back to 1000 AD? 800 AD? That could be interesting, right?

Steve McIntyre

#104. Roman, you’re undoubtedly right in how they should be calculated, but you just don’t see t-stats in the Mann corpus, so I’m guessing that the values are derived from the normal approximation. But we need to keep the correct t-formula at hand in the event that we get to trying to do these things based on methods known to the ROM.

Deep Climate

Off topic, but interesting from a blog management perspective – did anyone wonder about #81? The link goes to some Russian dining website. Looks like there’s a bot that took an earlier comment of pete’s, mangled the name slightly and added a link. What will those spammers think of next!

• pete
  Posted Jan 25, 2009 at 3:49 AM | Permalink

  Re: Deep Climate (#118), I certainly wondered!

Hu McCulloch

RomanM asks in #93,

I have to ask why in a serious scientific study, anybody would do a “two-sided test for code 5000/6000 proxies only (relatively small number of proxies)”. If you don’t know in advance whether the correlation is positive or negative (indicating that the physical processes behind that link are poorly and/or incompletely identified), how could you assume that that relationship is specifically to temperature in a linear fashion? Black box teleconnection at its finest!

Good question. However, its equally interesting converse is that if one is confident that a type of proxy (eg treerings) should have say a positive but noisy relationship to temperature, then the distribution of coefficients over a large number of such proxies should not be symmetrical about zero, but should at least have a clear majority of positives.

Mann’s 1209proxyames SI file does not give all the correlations, before or after his size-distorting Pick2 procedure, but only the ones that passed the usually 1-tailed significance threshold. However, Steve evidently found a version of the file that did have all the Pick2 correlations, and was able to compute the raw correlations himself. He reported the results graphically in his 9/20/08 thread on the The Mann Correlation Mystery.

The following scatter plot from that thread shows the raw correlations on the horizontal axis and Mann’s Pick2 correlations on the vertical axis.

The red symbols are dendro, which are supposed to have an a priori postive relationship to temperature. In fact, they look about equally positive and negative, with perhaps a majority of negatives!

He didn’t provide a separate histogram of the raw dendro correlations, but the right side of the following diagram is a histogram of all the raw correlations.

The little bump at the right end are the indisputably positive Luterbacher instrumental series. Aside from them, there is no evidence of a predominance of positive correlations, despite the fact that the supposedly positive dendros (TRW + MXD) make up the great majority of the proxies (927/1138 non-Luterbachers).

Mann’s own data therefore now places very heavy burden of proof on him and his coauthors to demonstrate that treerings or densities have a generally positive relationship to temperature at all.

Incidentally, the upper graph and the lower left graph illustrate that Mann’s Pick2 cherry-picking procedure, although highly size distorting, did not maximize his reported correlations as much as it might have: Instead of selecting the more positive of the two correlations, he selected the correlation with the greater absolute value, and then discarded that result if it was negative. As a result, the histogram on the left is hollowed out, but still fairly symmetrical after eliminating the Luterbacher bump on the right. See also my comment #12 on that thread.

• UC
  Posted Jan 25, 2009 at 11:09 AM | Permalink

  Re: Hu McCulloch (#121),

  before or after his size-distorting Pick2 procedure

  In propick.m there is option to select pick-n ( igrid , %numbers of gridpoints near proxy site used to cal. r )
  so I can try how would rtable look with igrid=1. But first I did
  compare my Matlab (igrid=2 ) result with the one Mann has in his
  webpage. Only difference is this:

  rtable_myrun(:,899)

  ans =

  -5
  37.5
  5000
  -0.072659745
  -0.27887841
  -0.12528468
  -0.50476717

  >> rtable_website(:,899)

  ans =

  37.5
  -5
  5000
  -0.149443
  -0.51617378
  -0.13864945
  -0.55486167

  is this result of some lat / lon issue discussed here ?

  Steve: Yep. That’s the rain in Spain in Kenya.

UC

Nah, not that easy when there are clear all commands everywhere and separate files for high – pick and low-pick.. Just a second..

Steve McIntyre

I’d got within .001 of UC’s values for 1193 out of 1207 proxies as follows. I’ve benchmarked my calcs against rtable results and they work in the cases that I’ve studies – I’ll need to check out the other 14. This calc uses the next1, next2 allocations that I calculated a while ago and there seems to be slight differences in a few cases with UC’s Matlab. As I recall, there were some seriously weird fencepost issues in Mann’s allocations.

url=”http://data.climateaudit.org/data/mann.2008″
source(“http://data.climateaudit.org/scripts/mann.2008/utilities.txt”)
download.file(file.path(url,”mann.tab”),”temp.dat”,mode=”wb”);load(“temp.dat”)
download.file(file.path(url,”Basics.tab”),”temp.dat”,mode=”wb”);load(“temp.dat”)
instrumental.info=Basics$instrumental.info; details=Basics$details

require(“R.matlab”);require(“Rcompression”)
download.file(“http://www.meteo.psu.edu/~mann/supplements/MultiproxyMeans07/data/instrument/infilled2poles_1850.mat”,”temp.dat”,mode=”wb”)
inst=readMat(“temp.dat”)
infilled2poles=ts(inst$anngridinst,start=1850); dim(inst) #157 2592

g=function(X) ts(X[,2],start=X[1,1])
library(signal)

details$pickone.1850_1995=NA
details$pickone.1850_1995lf=NA
for (k in 1:1209) {
(index= (1:2592)[!is.na(match(instrumental.info$jones,unlist(details[k,”next1″]))) ] )
#locates mannid for next1, next2 for proxy #1486 1487
x=g(mann[[k]]);y= infilled2poles[,index]
test=ts.union(x,y )
temp=(time(test)>=1850)&(time(test)=1850) &(time(x)< =1995)
temp1=(time(x)=1850)&(time(test)< =1995)
details$pickone.1850_1995lf[k]= cor(test[temp,],use=use0)[1,2]
}

test=read.table("http://signals.auditblogs.com/files/2009/01/rtable_pick1.txt&quot😉
test=t(test)
cor(test[,4],details$pickone.1850_1995) # 0.9985831
plot(test[,4],details$pickone.1850_1995)
temp= abs(test[,4]-details$pickone.1850_1995)<.001
sum(temp)
#[1] 1193
data.frame( details[!temp,c("id","pickone.1850_1995")],test[!temp,4])
id pickone.1850_1995 test..temp..4.
308 cana156 -0.05984178 0.074809533
395 fisher_1994_agassiz 0.16454857 0.258134410
396 fisher_1995_agassiz 0.04978613 -0.065326351
405 fisher_2002_plateau -0.20314180 -0.002479786
576 meese_1994_gisp2 0.01949094 0.061961720
641 mongolia-darrigo -0.03435543 0.091591579
799 nv037 -0.07755463 -0.040603843
897 qian_2003_yriver -0.18205392 -0.154840220
899 rodrigo_1999_andaluciarainual -0.12687449 -0.048791467
944 schweingruber_mxdabd_grid28 0.06002927 0.253837030
966 schweingruber_mxdabd_grid5 0.10559929 0.075073236
977 schweingruber_mxdabd_grid6 0.09591880 0.257575630
980 schweingruber_mxdabd_grid63 0.32504908 0.167564870
987 schweingruber_mxdabd_grid7 0.21058567 0.156203550
1017 smith_2004_thicknessmm 0.08221941 0.137785490
1055 thompson_2003_dasuopu 0.19079018 0.032087830

Hu McCulloch

Re my #97 and Steve’s inline reply,

I’ve now taken a look at Luterbacher’s paper (Science 2004 vol. 303 pp. 1499-1503, with Dietrich, Xoplaki, Grosjean and Wanner), and now concur with Steve that these series should not have been used in Mann 08.

Can anyone recommend a good recipe for roast crow? 😉

This is not to say that Luterbacher’s Figure 1C is not necessarily a valid multiproxy reconstruction of European temperatures back to 1500, to within its large CI. Nor that some or even all of his raw data might not be valid inputs to a study like Mann’s. The problem is that the entire period 1900-present is not a reconstruction at all, but merely the CRU and/or GISS data with which the earlier period was calibrated. Furthermore, the period 1750-1900 has a much smaller standard error than the period 1500-1750, because fewer good proxies are available in the earlier period. The proxies themselves might be valid inputs to another multiproxy study, and the overall reconstruction itself might validly be averaged into other, independent reconstructions, but the 71 Luterbacher local reconstructions cannot be treated as if they were proxies in their own right.

Luterbacher mentions that there is one instrumental record going back to 1659, and several back to the mid-18th century. Since these are not as comprehensive as the post-1850 or post-1900 data, they cannot just be averaged to create a precise instrumental index for all of Europe. However, there is no reason these shouldn’t just be treated as proxies in their own right, i.e. noisy observations on the average temperature. In fact, their calibration might not even require estimation of a slope term, since their slope relative to the comprehensive instrumental data in the calibration period is presumably unity.

Other Luterbacher input series, like historical records of cherry blooms and canal freezes, go back further, and may be useful proxies in their own right. Luterbacher’s SI unfortunately doesn’t tabulate the raw data, or even his local or continental reconstructions, but it does give a long list of reference sources that may have the actual data, if only in hard copy form.

Steve’s diagrams copied in #121 above illustrate the anomalous behavior of the Luterbacher series — The black symbols in the top graph are Luterbacher, and naturally have a much higher correlation with instrumental data in the 1850-present period, since they are the instrumental data in the post-1900 period, and presumably draw heavily on it in the 1850-1900 period. They account for the anomalous bumps on the right of the histograms in the lower figure.

Their high correlation with the instrumental data definitely cannot be taken as representative of the “skill” of the remainder of the Mann08 1209 “proxy” data set.

It is most unfortunate that a journal named “Science” would publish an article like this without at least requiring the authors to provide their numerical results! An electronic tabulation of the input data would even further advance the cause of true science (with a small s), even if the numbers can all be found in the 51 references cited in the SI.

Incidentally, the exponent in the formula I gave for the Šidàk size correction that I gave back in comment #96 didn’t come through correctly in the post, even though it looked OK on the preview window. Using ^ for exponentiation, the formula I gave that came out
1-.91/n = .0000871
should instead have been
1-.9^(1/n) = .0000871

UC

Team’s choice for H0 seems to be AR(1) noise , but I tried how 1/f – like noise ( as discussed in here and here ) would do, correlations for original proxy data and 1/f-noise :

Quite interesting, I’ll try to clean and put this somehow online, when I have more time. ( Note that tree rings are one example in here )

Steve McIntyre

#133. I’m not sure what the question is. In the figure in #130, the plateau of high correlations are the loaded dice from Luterbacher. (In the online R data frame details.tab, I’ve collated odds and ends of proxy info strewn in various places in the Mann SI. details$id[500:600] shows that these series are Luter series.

UC

Hu,

Is the big blip around observations 500-570 for real, or a typo?

These are sample correlations of 1/f noise vs. grid-box temperatures.

Steve,

In the figure in #130, the plateau of high correlations are the loaded dice from Luterbacher.

OK. So real proxies are not that significant, if you take 1/f as benchmark. Like I said, I’ll need some time to write this out.

Trying to go through relevant journals ( Technometrics, Biometrika, I don’t like Science or Nature anymore 😉 ) on this topic,

Slowly decaying serial correlations have, if ignored, disastrous effects on classical tests and confidence intervals.

..

Ad hoc methods like smoothing and detrending usually do not lead to satisfactory results.

(from A test of location for data with slowly decaying serial correlations, J. Beran, Biometrika (1989), 76, 2, pp. 261-9 )

Compare with Mann’s:

For the decadally resolved proxies, the effect is negligible because the decadal time scale of the smoothing is long compared with the intrinsic autocorrelation time scales of the data.

So the question is, are there slowly decaying correlations in proxy series ?

Steve McIntyre

#136. How retro. Of course, you’ll have to develop some kind of digital device to figure out what it means.

Steve McIntyre

Quite so – we do use “proxies” all the time. Brown and Sundberg describe a calibration theory. The problem with Mannian stuff is not the idea of a proxy, but that that he just calls things proxies. You don’t use the weight of the spectrometer as a “proxy” for iron concentration.

Sam Urbinto

Yes, proxies can be useful if they’ve been established to actually infer the variable of interest. My point is that if the timed vibrations of a piezoelectric or the effects of gravity on a fixed weight are the variables of interest themselves, they aren’t proxies. But in the other respect, everything’s a “proxy”, even the well known relationship between applied voltage and mechanical deformation.

Seems a rather pointless discussion. 🙂