Supplementary Information, McIntyre and McKitrick (2003)

Proxy Data

For the emulation of MBH98, the data set pcproxy used is here. This data set was originally obtained from Mann's ftp site here according to directions from Scott Rutherford (at Mann's direction), but this file was deleted after the publication of MM03. This file was only used for our emulation of MBH98. 

Recollated Tree Ring Data

We noticed problems with the principal component series in pcproxy.txt, most notably that the principal component series in this file were not correctly calculated. Accordingly, we carried out a complete re-collation of the tree ring proxy series used in MBH98 for 5 of 6 regions, using site identification information at the Nature SI and archived series at WDCP. 

In each of the 5 regions below, the "sitelist" is the list of WDCP identifications used in the collations and the "dataset" is the freshly collated tree ring information used in principal component calculations. For these calculations, the portions of the series prior to 1400 were not included in the collation. 

NOAMER (North America region): dataset sitelist  The Nature SI provided a list of ITRDB identifications. In nearly all cases, the SI identifications directly matched a WDCP site. A few aliases were noticed in the MBH98 listing and the following WDCP identifications substituted in our collation: ak001 for AK001C; wy005 for WY005B; the site AR045 could not be identified at WDCP.  

SOAMER (South America): sitelist dataset. Again, the ITRDB identifications are obtained from the Nature SI. The following aliases were identified: chil002 for CHIL002S; chil007 for CHIL007S;  

AUSTRAL (Australia-New Zealand) sitelist dataset. The ITRDB identifications were obtained from the Nature SI and no aliases were noticed. 

STAHLE/OK (Texas-Oklahoma) sitelist dataset. In this case, the Nature SI listed the names for 16 sites, but did not provide ITRDB identifications. The names for the 16 sites were looked up at the WDCP database and the identifications shown in the sitelist assigned and then used for re-collation.

STAHLE/SWM (Southwestern US- Mexico) sitelist dataset. The same procedure was followed as for the Stahle/OK series.  Subsequent to MM03, we noticed that a couple of sites had more than one version archived at WDCP and that the start dates of the version here did not match the start date of the version listed in the Nature SI. Subsequent work has used an updated version of this dataset.

VAGANOV: There was no site information for the Vaganov region. In this instance, we were unable to re-collate and used the 3 Vaganov series (#81-83) as supplied by Mann.

Updated Data Versions

For non-PC series, we noticed that archived versions at WDCP differed from MBH98 versions. MBH98 do not provide data citations (see AGU for a data citation policy) and the basis for data citation applied here is noted in the comments.

Retained Principal Component Series

For each region, we used the number of principal components listed in the weights file, which coincided with information at the Nature SI: NOAMER -9; SOAMER - 3; AUSTRAL - 4; OK- 3; SWM -9; VAGANOV -3. For each region, as stated in MM03, we used a conventional principal components calculation over the period for which all sites have values. Conventional principal component calculations, which MBH98 claim to use, require that there be no missing data.  In the calculations of MM03, we calculated principal components for each site over the maximum available period. The script is here. The relationship between the periods in which MBH principal components were applied (according to the pcproxy.txt file) and the period during which all selected sites in the region are available as shown here. (Postscript: subsequent to MM03, Mann et al. reported that they dealt with missing data in tree ring calculations by "stepwise" calculation of principal components in which the rosters were periodically changed. The difference in procedure resulted in the NOAMER PC1 and the STAHLE/SWM PC1 not being available to the AD1400 roster. This difference in procedure has been characterized by Mann et al. in very invidious terms as "selective censoring", but was actually a reasonable attempt to interpret MBH98 procedures from the public record, since Mann had refused to clarify procedures upon request. We have subsequently examined the principal component calculations of MBH98 in considerable detail and have carried out a reconstruction with the inclusion of the NOAMER PC1 and STAHLE/SWM PC1 in the AD1400 roster, showing the same results as MM03. This has been submitted to a senior journal and is under review as of June 2004.) 

Corrected Dataset

The corrected dataset incorporating updated versions of series as annotated above and freshly calculated principal component series as used in MM03 is here.

Calculations

To say that the description within MBH98 of their methodology is terse is an understatement. The algorithms described here were successful in replicating the MBH98 reconstruction in the 20th century to a high degree of precision; the correlations declined in earlier centuries, but still captured most MBH98 features. I (SM) notified Prof. Mann of this in a private email and requested clarification, perhaps in the form of a less terse description of the MBH98 methodology. Prof. Mann stated that he was unable to deal with this request because of numerous other time commitments. The algorithms described here are accordingly based upon a careful consideration of the public data and considerable experimentation.  (Postscript June 2004 - this represents the position as at October 2003 and is provided here to represent our approach as at MM03.)

(a) Temperature Principal Components

MBH98 purports to establish relationships between the proxies and 16 temperature principal components calculated from the Climate Research Unit (CRU) instrumental temperature database, using a subset of 1,082 out of 2,592 cells and the 79-year period from 1902-1980 as a calibration period.  MBH98 reported that they carried out a principal component analysis of the above temperature data, using "conventional" principal component analysis.  There is considerable missing data. Since "conventional" methods fail in the presence of missing data, the precise methodology used by MBH98 remains undisclosed. 

MBH have archived the gridcell locations, temperature principal components, EOFs and eigenvalues used in MBH98 at two locations: ftp://eclogite.geo.umass.edu/pub/mann/MANNETAL98 also ftp://ftp.ngdc.noaa.gov/paleo/paleocean/by_contributor/mann1998. These were downloaded and re-collated.  A script for downloading, collating and saving these and other MBH98 datasets is here.

MBH did not archive the gridcell standard deviations used to normalize the gridcell temperature series prior to principal component calculations. This information is required to calculate the NH average temperature after reconstruction of the temperature principal components using proxies. The current CRU dataset differs from the version used in MBH98. Because the now-obsolete version was unavailable, the gridcell standard deviations were estimated from the current CRU dataset. We do not expect the calculations to be especially sensitive to this approximation, or, alternatively, if they are, then this would raise questions about the robustness of the procedure. The CRU dataset was downloaded from CRU in July 2003 and collated into an R-table of dimension 1769 ( months from 1856 Jan to 2003 May) x 72 (longitude groups W to E) x 36 (latitude groups N to S).  The Jones marker of -9999 was changed to NA. The dataset was truncated to the 1082 MBH cells, the location of which was obtained from the MBH98 FTP site. Following MBH98, the dataset was then truncated to 1902-95 and then a Z-transformation (subtracting mean and dividing by standard deviation) was carried out cellwise. The scaling factor for each cell was saved as "mann.scale.tab".  Four cells were found to have no observations. These locations were saved as a tracking vector (of length 1082) is saved as "nil.tab". For each cell, the cosine of the latitude is used as an area-weighting factor. The scaled dataset obtained above is then multiplied by this area-weighting factor to obtain a scaled latitude-weighted temperature dataset, which is saved for further use.  This dataset is not used further in this paper, but was used in calculations regarding the MBH98 calculations of temperature principal components, which will be discussed in a later paper. The script for the above is here. 

(c) Reconstruction

The reconstruction of the NH temperature average is done through 4 functions.  

RPC  for a dataset of proxies proxy calculates the reconstructed TPCs.  The periodization defined in MBH98 is parameterized in the R-object period.m and the selections of eigenvectors in each period is parameterized in the R list select of length12. The proxies are Z-transformed basis 1902-1980. For all periods prior to 1970, the roster of proxies used in the period is defined as the set of proxies available in the first year of the period.  After 1970, the roster is the year-by-year available proxies.  For each period, using the roster of proxies and selection of eigenvectors, a matrix of coefficients G is calculated using Gcalc.  Then for each year in the period, the proxies are regressed against the selection of coefficients (using the proxy weighting factors as sent by  Prof. Mann's associate) thereby yielding the "reconstructed" TPCs for the selected eigenvectors.

Gcalc   (see after function RPC )  calculates a matrix G of regression coefficients given a roster of proxies (as a logical vector of length 112) and a selection of eigenvectors (as a logical vector of length 16) for the relationship between the proxies and TPCs for the calibration period 1902-1980. (All calculations are done after Z-transforming both proxies and TPCs to period 1902-1980, but this is not done in this function.)  This calculation is formally equivalent to MBH98's "least-squares solution to the overdetermined matrix equation.

 The script for generating Figures 7 and 8 is here. The MBH98 NH reconstruction (nhmann) is loaded, as are the MBH98 TPCs (pc) , eigenvectors (eof) , eigenvalues (lambda), weights (weights1), proxies (proxy), together with the above functions.  Figure 7a is column 2 of nhmann; Figure 7b is column 1 of rpc;  by applying the function RPC to proxy, a matrix MBH of reconstructed TPCs is obtained; Figure 7c is column 1 of MBH;  Then the corrected and updated proxy table proxy4 is loaded; a matrix new of reconstructed TPCs is obtained by applying the function RPC to the new table proxy. Figure 7d is column 1 of new. Figure 7 is then plotted.

The NH average using the reconstructed TPCs is now calculated. The reconstructed TPCs are expanded to 1082 cells using the eigenvectors eof, eigenvalues lambda, and the selection of eigenvectors for each period in list select. The MBH98 periodization is included as data within the function.  The location of the Mann cells is loaded within the function as is the vector nil of cells with no data.  This dataset is re-scaled using the scaling factors calculated in the original Z-transform of the CRU temperature data. The average is then taken of the dataset truncated to NH cells (using information on location of Mann cells). Areal weighting was already allowed for in the transformation described above. This gives Figure 8b; Figure 8a is the same as Figure 7a above.

Oct. 2003.

Feb. 4, 2005 - 2 hyperlinks edited.