Filling in Punta Laguna

Yesterday Mann archived a new version of his data, a version supposedly without any infilling. Gavin Schmidt is pretending that this was there all along and demanding an apology from one of his readers, who had the temerity to question Schmidt.

I’ve taken a quick look at the new (Sept 4, 2008) version of Mann’s Punta Laguna and plotted the 4 series below. It sure looks like infilling has taken place for this supposedly “original” data (though this would not be the case for annually resolved data.) The resolution of this data is a little sparser than many other series. This data also has a pronounced MWP depending on whether it’s flipped or not. (Correlation between O18 and temperature for one mollusk species is negative and one is positive, so I presume that Mann flips one of them over, rather than investigating whether the proxy has any merit if different mollusks yield different correlations.

Or perhaps individual mollusks have their own methods of teleconnecting to world temperature – different channels on the telecommunication band. FM88 – the CLAM. Or maybe P coronatus is sort of a guy’s channel, while C iloswegi is a comedy network. Or maybe the opposite. You need special crystals to tune your teleconnections.

Each of the 4 plots show the HadCRU3 gridcell temperature (annual) with the C13 and O18 values for the two mollusk species plotted using the right-hand scale. The solid dots show the actual values from WDCP archive – in each case, there are only 4 values for the entire 20th century (not decadal resolution). I have plotted the “infilled” Mann version in red, with the “original” Mann version in blue. Values from 1993 to 1998 have been removed in the “original” version; otherwise they are identical. [Update: thanks to an observation by Hu McCulloch, see comments below, the graphics below replace an earlier version in which 3 input values went to the wrong column in my read of WDCP data].

   
   

Figure 1. The above figures are matched to Mann proxies as follows. 381 – C13.C_ilosvayi; 382- C13.P_coronatus; 383 – O18.C_ilosvayi; 384 – O18.P_coronatus. This has been done based on the following nomenclature in Mann’s script call ids: curtis_1996_d13c, curtis_1996_d13cpyro, curtis_1996_d18o, curtis_1996_d18opyro. “pyro” is believed to be P. coronatus, leaving the ones with no appelation to be allocated to C. ilosvayi.

Update – The above tables hown the reported correlations, which are said to be “significant” even though there aren’t very many values. Here’s what Mann says about this:

we assumed n = 98 nominal degrees of freedom for annually resolved records, and n  8° 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 l p=0.10 value. 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.

These sentences bear some reflection in the context of this and other series.

17 Comments

  1. Sam Urbinto
    Posted Sep 5, 2008 at 9:25 AM | Permalink

    Ah, pretty pictures.

  2. Clark
    Posted Sep 5, 2008 at 9:36 AM | Permalink

    So, let me get this straight:

    1. There are series of unproven proxy value to temperature.

    2. The “original” data is actually a sparse list of real-live data with other “data” points added in by Mann, without explanation or rationale.

    3. The 20th century data is almost completely absent from some series.

    4. Mann simply creates data for these series into the modern period where they could be tested against detailed temperature records.

    5. For yet other records, the most recent real-live data points are deleted by Mann. These values are then replaced by others that Mann creates.

    Did I misstate anything?

  3. Richard deSousa
    Posted Sep 5, 2008 at 9:45 AM | Permalink

    As usual Gavin is playing the role of bully.

  4. deadwood
    Posted Sep 5, 2008 at 9:47 AM | Permalink

    The stick handling techniques of the Team are easier to follow now than they were 10 years ago.

    Its rather sad though that the National Academy of Science can be so easily fooled. Or are they complicit?

  5. Clark
    Posted Sep 5, 2008 at 10:03 AM | Permalink

    Deadwood: The PNAS is another odd journal because its referring is not run by a professional staff, but by members of the National Academy (which is basically a collection of high-profile scientists). If one of those Academy members is sympathetic to you or your point of view, you can get pretty much anything published. The journal has tried to add a layer of professional review on top of Academy member review, but it’s nothing like peer review at another paper.

  6. Jeff Alberts
    Posted Sep 5, 2008 at 10:08 AM | Permalink

    snip = I know Steve is always careful not to say that (and this post will probably get axed), but you have to call a spade a spade…

  7. Craig Loehle
    Posted Sep 5, 2008 at 10:11 AM | Permalink

    When one calculates R2 values with sparse data it is very important to calculate the adjusted-R2 instead, which subtracts for degrees of freedom. If you only have 3 data points and calculate an intercept and slope, you only have 1 df left and your adj R2 goes down precipitously. Thus in the above cases, the adj R2 would be much less than Mann reports, and in a case with only 1 data point would show that you have no predictive value (most software won’t even compute a regression using 1 point!).

  8. Not sure
    Posted Sep 5, 2008 at 10:53 AM | Permalink

    Gavin Schmidt is pretending that this was there all along and demanding an apology from one of his readers, who had the temerity to question Schmidt’s untrue assertion that the SI data at the time of his original comment was not infilled.

    Your tax dollars at work.

  9. Craig Loehle
    Posted Sep 5, 2008 at 11:47 AM | Permalink

    In the left hand panels, the fit would have been quite good if all the red dots were connected/used in the regression. This leads to two questions:
    1) why are the red dots from about 1900 to 1940 left out and replaced with an artificial line?
    2) Did Mann use the red line (or even blue line) interpolated annual “data” to calibrate/test the relationship instead of the raw data? That would be very naughty indeed…

  10. Posted Sep 5, 2008 at 12:33 PM | Permalink

    The use of isotopes as paleo thermometer is basically flawed as a affirming the consequent fallacy:

    If it’s raining then the streets are wet. The streets are wet.
    Therefore, it’s raining.

    Consequently: if it is warm then the isotopes are heavy; the isotopes are heavy but is warm?

    You can see here that higher humidity levels gives identical results to warm.

  11. Luis Dias
    Posted Sep 5, 2008 at 2:51 PM | Permalink

    Quite interesting. Is this phenomenon bounded to a single example, or the problem is in many other areas? Is this still unknown?

    Also, what do you mean by this?

    “It sure looks like infilling has taken place for this supposedly “original” data (though this would not be the case for annually resolved data.)”

    The infilled data is then what, monthly resolved data? Decadal?

    Thanks. Great work.

  12. Posted Sep 5, 2008 at 4:00 PM | Permalink

    Steve —
    The Punta Laguna series are admittedly goofy, but are not as bad as you have read them to be. Part of the problem is that you have misread the WDCP data file.

    This file has 7 columns, 1-3 for depth, age and date, 4-5 for C. ilosvayi, and 6-7 for P.coronatus. For 1985.3, 1972.1 and 1939.0, columns 6-7 are blank, while for 1923.6, 1915.8 and 1900.4, columns 4-5 are blank. Somehow your code has attributed the non-missing values to the wrong two columns, with the result that some values seem to appear out of nowhere, while others just disappear.

    Only series 383 (column 4, C.i d18O) and 382 (column 7, P.c d13C) passed the “screening”, and they did that for all 3 periods. These are graphed correctly below as red X’s, with Mann’s annual version as a blue line:

    Up to the last real observation (1993), Mann & Co have simply linearly interpolated the raw data to an annual frequency (fair enough), with no dropped or added values.

    However, their nonlinear extrapolations to 1998 are truly thin-air fabrications. I doubt that the brief extrapolations affect the results much for these two series, but if the same procedure was followed for other series, there could be a big problem.

    A further minor confusion is that in your Sept 3 post, the series are numbered 381-384 in the order of the columns in the WDCP file, when in fact Mann et al have for some reason listed these in their XLS file as numbers 383, 381, 384, and 382, in that order. The proxy numbers in the CA folder /data/images/mann.2008/proxy0381.gif, etc, do correspond the XLS sequence, however. (I was also thrown for a while by the fact that proxies 381-384 appear on lines 382-5 of the XLS file, but this is only because line 1 contains header labels so that proxy 1 appears on line 2 etc.)

    Steve: My bad. I had an error in my read-in code. Unusually for WDCP, the table is tab-separated (which I prefer and use myself, but missed). I should have used:
    yuc=read.table(url,fill=TRUE,skip=1,sep=”\t”)
    Instead I used:
    yuc=read.table(url,fill=TRUE,skip=1)
    I’ll re-do and re-issue. I allocated series 381-384 to WDCP columns (which I added to the caption and shold have done earlier) as follows. 381 – C13.C_ilosvayi; 382- C13.P_coronatus; 383 – O18.C_ilosvayi; 384 – O18.P_coronatus. This has been done based on the following nomenclature in Mann’s script call ids: curtis_1996_d13c, curtis_1996_d13cpyro, curtis_1996_d18o, curtis_1996_d18opyro. “pyro” is believed to be P. coronatus, leaving the ones with no appelation to be allocated to C. ilosvayi.

    • Craig Loehle
      Posted Sep 5, 2008 at 6:16 PM | Permalink

      Re: Hu McCulloch (#11), Hu: nicely done. Now, can anyone tell if the actual data or the interpolated data was used to regress against instrumental temperatures?

  13. Steve McIntyre
    Posted Sep 5, 2008 at 7:12 PM | Permalink

    Here is script for the above calculation (which requires that you have already downloaded the Hadcru.nc file. I’m experimenting a bit with temperature correlations and will report further on this.

    ##TEMPERATURE CORRELATION
    #spot check Curtis
    # id start end r1850_1995 r1850_1949 r1896_1995
    #381 curtis_1996_d13c -1597 1993 NA NA NA
    #382 curtis_1996_d13cpyro -1597 1993 0.3971 0.5158 0.5987
    #383 curtis_1996_d18o -1597 1993 0.6273 0.7467 0.8693
    #384 curtis_1996_d18opyro -1597 1993 NA NA NA

    #Mann PRoxy Data
    #download.file(“http://data.climateaudit.org/data/mann.2008/mann.tab”,”temp.dat”,mode=”wb”)
    #load(“temp.dat”); length(mann) #[1] 1209
    #returns list of 1209 tables with columns headed year, proxy, count
    #might change this to time series format at a later stage
    load(“d:/climate/data/mann.2008/mann.tab”)

    #MAnn Proxy Info #list of 1209 series
    load(“d:/climate/data/mann.2008/details.tab”)

    #HAdcru3 Temperature Data
    library(ncdf)
    v< -open.ncdf("d:/climate/data/gridcell/hadcrut3/HadCRUT3.nc")
    instr <- get.var.ncdf( v, v$var[[1]]) # 1850 2006
    #dim(instr)# [1] 72 36 1883
    #this is organized in 72 longitudes from -177.5 to 177.5 and 36 latitudes from -87.5 to 87.5

    read.hadcru3<-function (lat,long){
    j<- 19+ floor(lat/5 +.01);
    i<- 37+ floor(long/5 +.01);
    gridcell3<-instr[i,j,]
    index<-length(gridcell3)%%12
    gridcell3<-c(gridcell3,rep(NA,12-index) )
    close.ncdf(v)
    read.hadcru31850)
    yuc[temp,]
    # depth BP year O18.C_ilosvayi C13.C_ilosvayi O18.P_coronatus C13.P_coronatus
    #1 0 -43.000 1993.0 -0.56 -4.93 -2.36 -4.97
    #2 1 -34.940 1985.3 -0.54 -4.64 NA NA
    #3 4 -10.754 1962.1 -0.27 -5.05 NA NA
    #4 7 13.432 1939.0 -0.29 -5.06 NA NA
    #5 9 29.556 1923.6 NA NA -1.54 -5.07
    #6 10 37.618 1915.8 NA NA -2.98 -6.53
    #7 12 53.741 1900.4 NA NA 0.13 -4.34
    #8 13 61.802 1892.7 -0.98 -4.83 -2.21 -6.45

    #spot check Curtis
    # id start end r1850_1995 r1850_1949 r1896_1995
    #381 curtis_1996_d13c -1597 1993 NA NA NA
    #382 curtis_1996_d13cpyro -1597 1993 0.3971 0.5158 0.5987
    #383 curtis_1996_d18o -1597 1993 0.6273 0.7467 0.8693
    #384 curtis_1996_d18opyro -1597 1993 NA NA NA

    ####
    #PLOT COMPARISON
    gridcell=ts.annavg( read.hadcru3(lat=details$lat[k],long=details$long[k]) )
    k=381; curtis=”C13.C_ilosvayi”
    #k=382 ;curtis=”C13.P_coronatus”#Curtis ID #this is C13.P_coronatus
    #k=383;curtis=”O18.C_ilosvayi”
    #k=384; curtis=”O18.P_coronatus”

    mean.obs=mean(gridcell,na.rm=T);sd.obs=sd(gridcell,na.rm=T);
    c(mean.obs,sd.obs) # -0.03530296 0.38290608

    mean.est=mean(yuc[temp,curtis],na.rm=T); sd.est=sd(yuc[temp,curtis],na.rm=T);
    c(mean.est,sd.est) # -4.6383333 0.9693176

    y= (yuc[temp,curtis]-mean.est) *sd.obs/sd.est +mean.obs; y
    round(y,3)
    #[1] -0.166 NA NA NA NA NA NA -0.751 0.272 0.233 0.083 NA 0.118

    z=(mann[[k]]$proxy-mean.est)*sd.obs/sd.est+mean.obs;#z

    par(mar=c(3,4,2,4))
    plot(c(time(gridcell)),gridcell,xlab=””,ylab=”deg C”,axes=FALSE)
    axis(side=2,las=1);axis(side=1);box()
    points(yuc$year[temp],y,col=2,pch=19,cex=1.2)
    lines(mann[[k]]$year,z,lty=1,col=2,lwd=2)
    title(paste(“Punta Laguna “,curtis))
    at0= -9: -2; at1= at= (at0-mean.est)*sd.obs/sd.est+mean.obs
    axis(side=4,at=at1,labels=paste(at0),col=2,las=1)
    text(1930,-.7,paste( “r 1896-1995: “,round(details$rtable.r1896_1995lf[k],2)),font=2,pos=4,cex=.8)
    text(1930,-.9,paste( “r 1850-1995: “,round(details$rtable.r1850_1995lf[k],2)),font=2,pos=4,cex=.8)
    Data=readmann(k)
    z=(Data$proxy-mean.est)*sd.obs/sd.est+mean.obs;#z
    points(Data$year,z,col=4,pch=19,cex=.1)

  14. Steve McIntyre
    Posted Sep 5, 2008 at 7:36 PM | Permalink

    Here is a first calculation of correlations to gridcell annual temperature, making no attenpt right now to assess authocorrelation. Archived data is shown separately from the tables SD1.xls and r1209.xls.

    stat=as.matrix( details[381:384, c(22:24,32:35)]);stat=round(t(stat),3) ;stat
    stat=rbind(stat,array(NA,dim=c(3,4) ) )
    dimnames(stat)[[2]]=as.character(details$id)[381:384]
    row.names(stat)[8:10]=c(“emu1850_1996″,”emu1896_1995″,”emu1850_1949”)

    for (k in 1:4) { year=mann[[380+k]]$year
    temp1=(year>=1850)&(year=1896)&(year=1850)&(year< =1949);
    fred=ts.union(gridcell,ts(mann[[k+380]]$proxy,start=year[1]))
    stat[8:10,k]= round( c( cor(fred[temp1,],use=use0)[1,2], cor(fred[temp2,],use=use0)[1,2],cor(fred[temp3,],use=use0)[1,2] ),3)
    }
    stat
    curtis_1996_d13c curtis_1996_d13cpyro curtis_1996_d18o curtis_1996_d18opyro
    r1850_1995 NA 0.397 0.627 NA
    r1896_1995 NA 0.599 0.869 NA
    r1850_1949 NA 0.516 0.747 NA
    rtable.r1850_1995 0.132 0.329 0.432 -0.179
    rtable.r1850_1995lf 0.261 0.397 0.627 -0.270
    rtable.r1896_1995 -0.213 0.363 0.581 -0.188
    rtable.r1896_1995lf -0.293 0.599 0.869 -0.281
    emu1850_1996 0.135 0.419 0.622 -0.004
    emu1896_1995 -0.206 0.354 0.565 -0.182
    emu1850_1949 0.097 0.422 0.645 0.199

    `

  15. Posted Sep 8, 2008 at 11:07 AM | Permalink

    Re Craig Loehle, #11,

    Hu: nicely done. Now, can anyone tell if the actual data or the interpolated data was used to regress against instrumental temperatures?

    Thanks, Craig! As for your question, see my Comment #23 on the Proxy Screening by Correlation thread. In their SI, MZHBMRN08 indicate that they first interpolated the actual multi-year data to annual frequency, then ran an unspecified filter (perhaps Butterworth) on it to iron out fluctuations less than 20 years or so, and then computed their correlations. They used 15-2 = 13 DOF to interpret their full sample correlations, on the assumption that they had about 15 independent observations on each series after this filter. But when, (as with Steve’s example of Punta Laguna) there were in fact only 9 or 10 observations to start with, this would greatly understate the critical value, and so greatly overstate the significance of the correlation.

    In Loehle and McCulloch (2008), the low-frequency data were likewise first interpolated and then smoothed (with a 29-year centered mean), but since each series had already been calibrated to temperature by their authors, this procedure did not affect the validity of these calibrations. My 95% CI’s for the mean of the up to 18 series did use t critical values that assumed that the variance of each series about the global mean was computed from about 60 independent tridecadal observations, but the effect of DOF on t critical values is much weaker than the effect of DOF on R2 critical values.

    As I mention in the comment linked above, in order to make correct inference about any calibration of the MZHBMRN08 series, the handful of actual observations for each series should have been regressed on the corresponding (local or global) temperatures. Far fewer than 484 would have appeared significant, and this is particularly true of the non-dendro subset.

  16. Posted Sep 27, 2008 at 7:23 PM | Permalink

    Thanks Hu,

    I stared at these graphs for a long time trying to understand them before. I could see the shape wasn’t right like the wrong curve was on the wrong graph, now it makes much more sense.

    I came back today after my m08 reconstruction of the reconstruction showed how important these series are.