Tuesday, April 13, 2010

The bunnies have been asking what's wrong with the rather strange roughly cylindrical brown object that Craig Loehle left lying on the sidewalk at Atmospheric Environment. The constant nattering has woken Eli from his retirement nap so he might as well have a go.

This is a case, where when you see the trick you want to go punch the referees, decapitate the editor and shoot the author into space with a one way ticket. It starts with a paper by Hofmann, Butler and Tans in the same journal which looked at how fast the CO2 mixing ratio has been growing since ~1800. They find a close correlation between population growth and the mixing ratio, indeed, in 1987 (note added in proof) Newell and Marcus had suggested that monitoring CO2 mixing ratio growth was a reasonable way of monitoring global population.

HBT show that the CO2 mixing ratio has been growing very close to exponentially since ~ 1960 when monitoring at Mauna Loa started, although going back further than that on the basis of ice cores shows slower growth in earlier times. They then try to tie everything together to suggest

Besides showing the insight gained by removing pre-industrial CO2 in explaining the curvature in the Mauna Loa CO2 record, and organizing the confusion on growth rates, does this new analysis of the CO2 record have any further use for the future? One possibility is to use it to watch for the expected and necessary break of the exponential nature of CO2 and its growth rate. Fig. 4 shows an exponential fit to the past 14 years of Mauna Loa anthropogenic CO2 data and the residual between the data and the exponential function. The residual has varied from zero less than 1% over this time period. A sustained downward trend of the residual by more than 1% (in the absence of any known major volcanic event) would be a clear sign of a change in the nature of anthropogenic atmospheric carbon dioxide, and a possible gauge of progress in the inevitable need to limit atmospheric CO2. Similarly, a break in the close relation between anthropogenic CO2 and population would signal that a change in ‘‘business as usual’’ had occurred.

For anyone who needs a shorter HBT, they said watch for deviations from exponential growth on short to medium term intervals to spot changes in emissions. They did not predict future emissions but they did discuss the entire record. Let Eli see what Loehle did according to Loehle

A paper by Hofmann et al. (2009, this journal) is critiqued. It is shown that their exponential model for characterizing CO2 trajectories for historical data is not estimated properly. An exponential model is properly estimated and is shown to fit over the entire 51 year period of available data. Further, the entire problem of estimating models for the CO2 historical data is shown to be ill-posed because alternate model forms fit the data equally well. To illustrate this point the past 51 years of CO2 data were analyzed using three different time-dependent models that capture the historical pattern of CO2 increase. All three fit with R2 > 0.98, are visually indistinguishable when overlaid, and match each other during the calibration period with R2 > 0.999. Projecting the models forward to 2100, the exponential model comes quite close to the Intergovernmental Panel on Climate Change (IPCC) best estimate of 836 ppmv. The other two models project values far below the IPCC low estimates. The problem of characterizing historical CO2 levels is thus indeterminate, because multiple models fit the data equally well but forecastvery different future trajectories.

Watch the pea. The IPCC estimate is for A PARTICULAR EMISSION SECENARIO, it is not a projection of the last fifty years of measuring CO2 mixing ratios, nor is it the only emission scenario. Fast work there folks. Still, Eli is a trusting bunny so let us go on. The three fitting forms and Loehle's values for the parameters are

and this is what Craig shows,but Eli made a somewhat fuller picture, using the Loehle's values (click on the graphs for better views)

The exponential model, although not perfect, does match the CO2 mixing ratio much better at earlier times than 1959, so, contrary to Loehle, the three models are not equally plausible. But there's more. Eli also did his own fit.

The saturating curve here is much closer to the data and the exponential curve. Still, young bunnies, there is still something like junior high school wrong with what Craigie poo did. Flip to the continuation for the answer

------------------------------------------------

Loehle extrapolated all of his fits from the year dot. Little ones can spot this in his value for b in the exponential fit and for a in the quadratic and d in the saturated one. If the time is replaced the three fitting forms by [t-1985] you get much much improved saturated fits. The quadratic fit still sucks at early times (and late ones) but even it comes closer to the other two..

The extrapolations are shown below.

Model

Year 0 Loehle

Year 0 Eli

Year 1985 Eli

2000

Exp

370

369

369

Quad

367

367

370

Saturated

368

370

370

2050

Exp

513

513

510

Quad

479

479

490

Saturated

467

500

512

2100

Exp

846

846

819

Quad

651

651

668

Saturated

567

729

834

The referees and the editor missed the slight of hand, they also did not read Hoffman, et al. The graphs demonstrate that quadratic fits are not industrial grade for this problem, and the table and graphs show that Loehle's fits are fishy. Doesn't anyone check the programming? Where are the auditors when you need them.

44 comments:

Up late working and saw your post when taking a break. Hope my post at Tamino' didn't come across as pushy, it just occurred me at the time that you'd mentioned taking a look on an earlier thread.

Regarding Craig's analysis. Wow....words fail me.

If I were a brave bunny I'd hope over the CA and ask Stephen McIntyre to do an official audit on Loehle's paper. But I'm a meek and scared lagomorph these days, and the acolytes at CA make my ears twitch, it is really annoying.

It's still early here, but I don't understand. There are already t-basing constants (d) in the saturated equation, and in the quadratic one (a and b), and in the exponential one (b). In other words, if you can get a better fit by replacing t with t-1985, then you can get exactly the same better fit by choosing different values of these constants. Can't you?

Eli, forget my question. I got this working in matlab using lsqcurvefit, and replicated your findings... except for the 'saturated' equation.

I find that its solution is extremely poorly conditioned and sensitively dependent on getting the numerics right. For me, it converges -- very very slowly -- to the quadratic solution, not, as Eli found, to the exponential one.

I suspect there may be more than one local optimum for this function, all very very flat, and which one you get depends on your starting point and whether you get there all the way, on not prematurely stopping the iteration.

Nick, I think it is only the saturated equation that has a t shift, and your question is legitimate there. I suspect that the reason for Eli getting a different result here is more related to the above mentioned fucked-up-ness of this function.

The sensitivity of the saturated model was my initial concern, if it closely matches an exponential over the calibration period that would suggest very little constraint on the saturation. It is possible that Loehle's model fit was heavily determined by the initial values of the parameters, which just by chance gave a low projection. My though was to perform a Markov Chain Monte Carlo simulation and see what the posterior distributions looked like. I suspect for the saturated model they are rather broad.

that's precisely it. You need three free parameters to fit the calibration period, and three parameters is sufficient for a smooth curve like this. The fourth parameter is a giveaway that there is a near-singularity somewhere.

To Nick Barnes: yes, you can choose the "reference point" for time (t0) to be anything you want, and adjusting the constants still gives the same fit.

BUT: when you use just "t" (i.e., t0=0) with all t values very near, say, 2000, then terms like t^2 are all near 4000000, and small differences between t^2 values which are important, get lost in the trailing significant digits. So, unless your computer tracks 30 or so significant digits, your calculation of the parameters can be highly inaccurate.

It's often a good idea, when curve-fitting, to transform "t" to "t - mean(t)" (in the lingo of R) so that the mean value of the predictor variable is zero.

OK, so we're looking at the rounding effects. But the biggest difference is in the saturated case, where it's completely unclear to me either what the equation is supposed to represent (what is "saturated" about it) or how Loehle chose his parameters. I note that Loehle's curve is the other way around from Eli's (b values have opposite signs). If we constrain ourselves to negative b, what is the best fit?

Nick Barnes, I think the point is that there are many different values for the parameters of the saturated model that all give almost equally fits to the data, but which give very different extrapolations. This happens because there is a flat minimum (or perhaps several flat minima) to the cost function used to fit the model.

I was also originally suspicious of using the data with the annual cycle left in, as the end-points may have an substantial effect on the extrapolation. Ideally if you are going to write a comment paper, you need to justify any changes you make in the method used in the original paper (and show that the results are robust to your choices).

My paper was about forecasting, not back casting. In econometrics a model forecasting cell phone sales need not back cast to be useful. Not enough digits were carried in the publication to replicate the graphs. It is suggested that a high precision platform such as Mathematica be used with the following parameters: exponential {a=259.3945817870858, b=2.9782774195066923^-13, c=0.016769296042440738}, quadratic {a=45318.539465496644, b=-46.7684328499448,c=0.0121468549016735}, and saturating {a=616.9635436975683, b=0.006760893860992654, c=310.9394987790752, d= 1945.6473702994858`}.Craig Loehle

Craig, the necessity for so many significant digits (in some case more than are supported by double precision format) seems to me an indication that the problem as posed was ill-conditioned. The fact that making a better conditioned problem by changing the reference time gives a different result suggests that it is your model fit that is unreliable, rather than Eli's?

I disagree about the back casting comment, as a statistician what matters is decided by the needs of the application, not the statistician. Just because it doesn't matter for cell phone, doesn't mean it doesn't matter for carbon dioxide forecasts. If one model can give back-casts, but the others can't that immediately means that there is a reason to chose one model over the others, and the paper ought to have pointed out that deficiency of the saturating and polynomial models.

Having now written some MATLAB code, I am not so sure the reference time is the issue. I used the yearly mean Mauna Loa data (as I don't think the annual cycle should be allowed to influence the outcome), and using a reference time of 1985 I get a plot that looks a lot like Craigs.

I suspect the problem is flat minimum, so maybe tomorrow I'll find time to try the MCMC approach and look at the uncertainty in the forecast (which I would imagine are important in econometrics as well).

Having said which, I have realised that the back-casting is a red-herring; the reason the data before 1958 are important because if you built a model using the available relevant knowledge (which includes the Law Dome data for instance) the quadratic and saturating models would no longer give the same forecasts, whereas the exponential one would be largely unchanged. I have more trust in a forecast from a model the more observational data it explains, especially if it explains it without re-calibration.

The key point is "would the comment get past the reviewers had the back-casts been plotted?". I suspect the answer is probably "no".

Eli agrees with Dikran at least in the important part The basic point is that forecasts are based either on historical data or models. In this case the forecast was based on historical data. The strong deviation from from the data of two of the fitting forms is equivalent to the 1998 cherry pick we are all familiar with.

On the centering of the fit, not so much. The farther the extrapolation, the more sensitive the fit is to noise in the data, lack of infinite precision, etc. So, saying this awkwardly, a small movement in one parameter can lead to a large on in another if you don't choose to wisely.

"259.3945817870858" taken from Loehle's paper. When a person quote many (many) more digits than are significant, it's usually a good indication that he does not know what he's doing re math, curve-fitting, or statistics.

Here is my code, which reproduces something like Craigs results, but with annual mean data, and using a reference time of 1985 to reduce numerical problems. It is all self contained, so just copy all into a single file and run. I get the following 0.5 times sum of squared errors:

I agree. I find it utterly bizarre that Craig's formulae apparently require astronomical numbers of sig fig to produce results that need only have only 4/5 sig fig resolution. Accuracy is another matter, and seems to be poor in any case.

All this curve fitting reminds Horatio more than a little of running on a hamster wheel.

Consider what has happened to the growth in (yearly) emissions from the 1990's to 2000's decribed here

"Worldwide growth in carbon dioxide emissions has doubled since the close of the 1990s, reports a study published in the early on-line edition of the Proceedings of the National Academy of Sciences."

"A team of researchers led by Michael R. Raupach of CSIRO-Australia found that global growth rate of CO2 increased from 1.1 % per year during the 1990s to 3.1% per year in the early 2000s. The growth is being fueled primarily by developing countries, especially China and India, where economies are fast-expanding and population continues to increase at a significantly higher rate than in industrialized nations."

There are many many many things that can affect the actual emissions path that is taken, so many, in fact, that it is virtually impossible to say with any confidence what the total emissions -- and resulting increase in atmospheric CO2 concentration -- will be by 2100.

That is precisely why the IPCC gives many different scenarios. As Eli indicates with his "watch the pea" comment, they are simply not in the business of "predicting" future emissions, nor should they be.

For a particular scenario, they make an assumption about what will happen over the next century to yearly emissions (eg "They will grow by 2% per year" or whatever) and then calculate, based on the details of the climate system (eg how much CO2 is absorbed by oceans and other natural sinks, etc) how much that will increase the atmospheric CO2 concentration.

Attempting to actually "forecast" future emissions (especially based on simple "curve fitting") is basically a hamster's game. It tires you out but really gets you nowhere (Horatio knows about the latter)

Having said that, Horatio has run the numbers and thinks it should be "a=259.3945817870856"

Now if I were Jonathan Leake, then tomorrow's headline might read something like this:

"Craig Loehle throws McIntyre under the bus"

In a startling revelation yesterday, Craig Loehle said that Stephen McIntyre uses inferior statistical software which is not suitable for conducting precise statistical analyses. This now calls into question the validity of all of Mr. McIntyre's audits. When approached for comment, Mr. McIntyre muttered something about whitewash and the Hockey Stick being a fraud before breaking down in tears. After composing himself, Mr. McIntyre said that Loehle has now been banned form his blog."

Being extremely suspicious of Craig's huge numbers of sig fig, I thought I'd try a simple test using the exponential formula.

I couldn't use the the full number of sig fig in my spreadsheet. I lost 1 from a, 2 from b, and 2 from c. But anyway, the fit 1958-2008 was not particularly good. I got Total squares error 59.931958 value 313.56 (actual 315.34)2008 value 384.67 (actual 385.59)2100 value 845.36

I then removed a further 3 sig fig from each parameter. The results were Total squares error 59.931958 value 313.56 (actual 315.34)2008 value 384.67 (actual 385.59)2100 value 845.36

So, only now do we see any difference. In fact, if I add 2 sig back in, the result is Total squares error 59.871958 value 313.56 (actual 315.34)2008 value 384.67 (actual 385.59)2100 value 845.36At this point, a=259.39458b=2.9782774E-13c=1.676930E-02

IOW, you only see any difference when using 8 instead of 9 digits. This suggests that Craig's use of such "precision" makes no sense and that it's a smokescreen for a poor analysis.

Anyrabett who teaches an introductory science class sees this frequently. It is a sure sign of someone who is unclear on the concept.

OTOH, as Tamino pointed out above, many fitting procedures (OLS) require subtracting numbers from each other where there is a small difference between them, and the procedure requires many significant figures. Experience is knowing when to cut and when to paste.

Rabett Run

Subscribe Rabett Run

The Bunny Trail By Email

Contributors

Eli Rabett

Eli Rabett, a not quite failed professorial techno-bunny who finally handed in the keys and retired from his wanna be research university. The students continue to be naive but great people and the administrators continue to vary day-to-day between homicidal and delusional without Eli's help. Eli notices from recent political developments that this behavior is not limited to administrators. His colleagues retain their curious inability to see the holes that they dig for themselves. Prof. Rabett is thankful that they, or at least some of them occasionally heeded his pointing out the implications of the various enthusiasms that rattle around the department and school. Ms. Rabett is thankful that Prof. Rabett occasionally heeds her pointing out that he is nuts.