### The Battle of the Models - Keenan vs. Met

The latest contestants to enter the climate debate stage are mathematician Douglas J. Keenan and the UK's Meteorological Office or Met for short. On November 8, 2012, Lord Donoughue, a member of the UK Parliament posed the question:

*To
ask Her Majesty's Government ... whether they consider a rise in global
temperature of 0.8 degrees Celsius since 1880 to be significant. [HL3050].*

The Met Office answered, "the temperature rise since about 1880 is statistically significant," meaning that the temperature rise since 1880 is attributable to human caused global warming.

Keenan challenged the Met Office's reply to Parliament on the basis that the model chosen by Met's leading statistician, Dr. Julia Slingo, was not suited to answer the question. Keenan offered up his own choice of model that he claims was much better suited. As you might imagine, a series back and forth opinion editorials and commentaries ensued which sparked a rather heated debate.

*Douglas Keenan – Editor of
InforMath.org *

Much has already been written about the tone of the debate, Met's failure to provide specifics on their claims of statistical significance, and of Keenan's brazen indictment of Met's Chief Executive Officer, John Hirst and leading scientist Dr. Slingo. It is the ruckus between Keenan and the Met that has captured most people's interest and is the basis of most articles on the matter.

Being a numbers guy, I'm less interested in the tone of the debate and more interested in the models that the Met and Keenan espouse. I'll lay out the strengths and weaknesses of the two models and offer a closing assessment of which model is more suited to determining the statistical significance of the temperature rise of 0.8 degrees C since 1880.

Statistical significance is a somewhat nebulous term when it
comes to non-stationary time series data such as the temperature record. **Most measurements of processes in nature
will show statistically significant results if the time series is long enough
simply because Mother Nature does many cyclical or curvilinear things in big
ways.** The following graphic illustrates this.

*The Pacific Decadal Oscillation
is a natural cycle with significant trends*

In medical research, statistical significance is usually determined by comparing a control group to an experimental group. However, in climate research there is no control group, making the testing for statistical significance a more mathematically tangled affair, which sparks much disagreement among scientists on how to best do it.

The central debate between Keenan and the Met Office revolves around which model provides a better "goodness of fit" (GOF) to the temperature record and has the best predictive skills from which to infer statistical significance. Rather than focus on the question of statistical significance in post-industrial temperatures, which I'll leave to Met and Keenan to debate, this article takes an objective look at the two models they espouse and why one model is better suited to determine statistical significance.

The Met uses the AR(1) model. They claim their model shows statistical significance in current temperature trends, inferring recent temperature trends are being caused by human activities. Keenan's charge is that AR(1) is not robust enough to make such determination. He espouses the ARIMA(3,1,0) model which he claims is more robust. He claims his model shows today's temperature trends are not statistically significant and more likely caused by natural variability.

**Decomposing The Models**

Both models begin with the letters AR, meaning that both
models employ **A**uto-**R**egression, a stochastic process that
weights the sum of previous measurements in the time series and white noise
error to produce a predictive output for the next measurement. In the case of
Met's AR(1) model, the one (1) means it is a "first-order autoregression" model.
In practical terms, first order means the temperature from January of last year
has a direct impact on January's temperature for this year but the temperature
from January two years ago or earlier has no direct relationship to this year.

The use of the AR(1) model is premised on the assumption that changes in temperatures have only first order relationships as described above. That's a rather bold and unjustified assumption in my opinion. There are known natural cycles and events, lasting from a few years to many decades, which affect the earth's temperature. The notion that temperatures more than one year back can have some effect on today's temperature seems intuitive to me. This is one of Keenan's criticisms of the AR(1) model – a criticism that I agree with.

The ARIMA(3,1,0) model proposed by Keenan is parameterized to use third-order autoregression as indicated by the value three (3). Keenan's model assumes the three previous temperature measurements and the white noise term contribute to the current temperature measurement. If the sum of the three previous temperature measurements is positive, the output will resemble a low pass filter with the high frequency part of the white noise decreased in amplitude. We'll see that low pass filtering effect when we start looking at comparative graphs of AR(1) vs. ARIMA(3,1,0) model output.

There are three parts to the model: "AR" (**A**uto**R**egression), which has been discussed above. The "I" (**I**ntegrative) term indicates that the
model applies an integrative function to change the data from non-stationary to
stationary structure to make it more suitable for autoregression. Keenan's
model applies one (1) order for differencing. The "MA" (**M**oving **A**verage) term
indicates the model applies a moving average to its output. Because Keenan's
value is zero (0), no moving average is applied to the model's output.

Keenan's model is flexible and more adaptive but considerably more complex. This is the chief criticism of Met regarding Keenan's choice of model. I'll be exploring that criticism shortly.

The Battle of the Met's
AR(1) and Keenan's ARIMA(3,1,0) Models

To determine "the better model," there are a number of objective measures we can look at to assess the model's parameterization, goodness of fit to observed temperatures, and predictive skills. In Keenan's supplementary material, he explains that he tested the Met's AR(1) and his preferred ARIMA(3,1,0) models against the GISS temperature series. In Met's response to Parliament of Keenan's criticisms, they used the HadCRUT4 temperature series.

I personally don't think which temperature series is used to test the models matters. If it does, then we have other systematic issues to deal with that go beyond the choice of models. For the purpose of this article, I'll use the latest GISS GHCNv3 temperature series to compare the two models. All of my statistical analysis is done using the statistical package R with the "Forecast" library – a package designed for AR and ARIMA modeling.

**Let the Analysis Begin!**

**Test for Statistical Significance
–** I performed appropriate tests to determine if the results I was
getting were, themselves, statistically significant and could be relied upon. Again,
**I'm not testing the statistical
significance of the temperature series but rather I'm testing the reliability
of the tests I'm performing on the models.** This is standard practice in
mathematics. In most cases an ANOVA test was performed. In all cases, I tested
at the 95% Confidence Interval (*alpha*
= 0.05, *p* < 0.05). The tests I
performed on both models were as follows:

**Sigma^2** - When we
compare the fit of the models to the data, we look for the smaller standard deviation
of the residual (error) of the model. The residual represents the portion of
the data the model can't explain. By "can't explain," I mean the model fails to
account for. The lower the residual, the better the model fits to the data and explains
it. The standard deviation of the
residual is called the Sigma^2.

**Akaike Information
Criterion test with correction (AICc)** – AICc is a method of testing a
given model from a set of models. The best model is one that minimizes the
Kullback-Leibler distance between the model and the data. Kullback-Leibler
distance measures the difference between the probability distributions of the
model as compared to the actual temperature series distribution. The model that
has the best GOF (or least divergence between the two distributions) is the better
model.

**Bayesian Information Criterion
(BIC)** – The chief criticism Met levies against Keenan's choice of
model is the greater number of parameters used in his chosen model makes it unnecessarily
complicated and prone to error. The
BIC test is similar to the AICc test but penalizes the model for over-parameterization
(being unnecessarily complicated for the data under analysis). It is a better
test when comparing models with differing orders of parameterization as in our
testing.

**Coefficient of
Determination (R ^{2}) and Residual (1-R^{2})** - The R

^{2}measure determines what percentage of the variance in the trend produced by the model output can be explained by the model. It is a measure of how well the model is designed and parameterized to account for the variance and trend in the data. The residual of the coefficient, 1-R

^{2}, tells us what percentage of the variance in the model trend is due to residual (error) and is not accounted for by the model. The higher the R

^{2}and the lower the 1-R

^{2}, the better the model is designed and parameterized.

** F Ratio** - The

*F*ratio is part of the ANOVA test of the model's performance. It tells us the ratio of the variability of the trend explained by the model to the unexplained variability of the trend, each divided by the degrees of freedom (

**) of autoregression used by the model. The higher the F ratio, the more useful the model is to analyze the temperature data.**

*df***Lag Correlations**
– The lag correlations tell us how well the parameterizations of the
models perform their autocorrelation to the time series data. There are two
objective measures of lag correlation: ACF (residual **A**uto**C**orrelation **F**unction) and PACF (residual **P**artial **A**uto**C**orrelation **F**unction). Collectively, these two measures
determine how well the parameters to be used for the "AR", "I", and "MA"
functions perform for the respective models. They can also be used to test for
the best parameters to be used. ACF and PACF are best analyzed visually in
graph form.

I could write a book about how to interpret ACF and PACF.
So, to keep it simple we'll look at one important aspect of the above graphs.
Met's AR(1) model is on the left side with ACF at the top and PACF at the
bottom. Keenan's ARIMA(3,1,0) model is on the right. Notice that each graph has
dotted blue lines. This is referred to as the cutoff band. When a model's
parameterization is optimum, its lag values will remain inside the cutoff. **A special note is in order here. In R, a Lag
0 line is drawn for the ACF graph that extends well outside of the cutoff
lines. This is a well-known defect in the ACF graphics function so it can be
ignored.**

Note that Met's AR(1) model makes excursions outside of the
cutoff lines. Keenan's ARIMA(3,1,0) model's lag values remain inside the cutoff
lines, indicating that **Keenan's model autocorrelates
better **to the time series data.

**Diebold-Mariano Test**
– This tests the null hypothesis that the two models perform equally as
well against the same time series data. The test compares the residuals (error)
of each model. The results of the Diebold-Marino test are as follows:

DM = -2.4648

Forecast horizon = 1

Loss function power = 2

p-value = 0.01499

alternative hypothesis: two.sided

The Diebold-Marino test is statistically significant (*p* = 0.015), meaning we reject the null
hypothesis that both models perform equally. The DM value is negative,
indicating that **Keenan's ARIMA(3,1,0)
model performed better than Met's AR(1)**.

**The Box-Ljung Test
for Residual Dependency** – The Box-Ljung test is a diagnostic tool
used to test the null hypothesis that the model does not exhibit a lack of fit
to the time series. The test is applied to the residuals of a time series after
fitting an ARIMA(p, d, q) model to the data. The test examines autocorrelations
of the residuals for dependency. If the autocorrelations are very small, we
conclude that the model does not exhibit significant lack of fit.

This is somewhat inverse logic and a little bit difficult to follow. Ideally, you want the model to fail the Box-Ljung test (have a p-value > 0.05) supporting the null hypothesis that the observed fit in the model is due to the model's predictive skills (I'm oversimplifying here but hopefully you mathy people get the point).

**The Proof Is
In the Tasting of the Pudding**

So far, I've thrown out a bunch of model tests and made my best attempt to explain what they mean. If you've reached this place in the article, congratulations, you managed to stay awake, LOL! For those who found the above boring, stick around. It's going to get more fun from here.

The proof is in the tasting of the pudding. Both the Met's AR(1) and Keenan's ARIMA(3,1,0) models are time series predictive models, meaning they are mathematical crystal balls that can tell the future, given some hint of the past.

This is what the models render when we plug all of the temperature data into them:

The actual GISS GHCNv3 temperature trend is the black line in the graphic. Met's AR(1) model output is in red. Keenan's ARIMA(3,1,0) model output is in green. If you study this graph closely, you might notice something odd towards the current timeline in the graph. Here, let me show you:

One of the criticisms of the AR(1) model is it's ability to follow a non-stationary trend falls off as the time series gets longer. Above, you can see that up until 1960, the AR(1) model tracks the actual temperature trend fairly well. However, as we approach the year 1990 and forward, it starts to drift negative in relation to the actual trend. It is a "drifting" model whereas Keenan's ARIMA(3,1,0) is driftless.

Earlier in this article, I mentioned that the ARIMA(3,1,0) model uses third order auto-correlation. In mostly positive or mostly negative trends it tends to act somewhat as a low pass filter, meaning it favors the lower frequencies of the trend and tends to be less reactive to the high frequency variability components of the trend. You can see that the green line of the ARIMA(3,1,0) model is less variable than the AR(1) line. The key performance indicator of these models is how well they predict the trend and not so much how well they account for high frequency components in the trend.

Now, we're going to see how well the two models serve as crystal balls. In this next test, I removed the years 2002 through 2012 from the GISS GHCNv3 temperature series. I intentionally selected 2002 as my truncation date because the years from 1997 through the end of 2001 provide a "hint" that the rate of global warming was beginning to stagnate. I wanted both models to have equal opportunity to get the hint.

The following two graphics overlay the model's predictions over the actual temperature trend which includes the 2002 through 2012 actual temperature anomalies.

**First up, Keenan's ARIMA(3,1,0) model prediction for the
years 2002 through 2012:**

**And now, Met's AR(1) model's prediction of temperatures from
2002 to 2012:**

In both graphics, the dark blue lines are the model's temperature prediction. The light blue bands around the prediction lines are the 95% Confidence Interval. The still wider gray bands are the 90% confidence Interval.

It is obvious that Keenan's ARIMA(3,1,0) model took the hint I provided and modeled the temperatures from 2002 to 2012 with an accuracy of close to 0.05 deg. C. Met's AR(1) model completely missed the hint and assumed that temperatures would follow the momentary down trend that was occurring at the truncation of the temperature series.

An interesting surprise to me was that Keenan's ARIMA(3,1,0) model actually predicted the two peaks and one trough in temperature trends during the prediction period, albeit somewhat muted as expected. You can see the two small peaks and single trough in the dark blue line in the ARIMA(3,1,0) model prediction. Well done!

**Is There a Better Model?**

A model can always be improved upon. We hear that from climatologists all the time, right? Well, actually it's true. Referring back to my discussion of ACF and PACF, I played around with the ARIMA model class and found a set of parameters that performed slightly better than Keenan's ARIMA(3,1,0) model. My model was the ARMIA(0,1,1).

When I say slightly better, I mean my ARIMA(0,1,1) model failed several tests for significance in performance against Keenan's ARIMA(3,1,0) model, meaning that while some of my key performance metrics were ever so slightly better, they weren't better enough to be of any interest outside of the fun I had playing with it. Here's a graph of my model compared to Met's and Keenan's.

The black line is the actual temperature trend. Red is Met's AR(1) model. Purple is Keenan's ARIMA(3,1,0) model, and green is my ARIMA(0,1,1) model. As you can see, my model tracked Keenan's almost identically and was not significantly different in performance. I ran the 2001 to 2012 predictive test on my model and, while I got a linear trend similar to Keenan's, mine failed to predict the peaks and troughs because it lacked autoregression to learn patterns from the past.

**Which Model Won?**

Forgetting about my uninteresting model we're going to come full circle in the central debate – the question of which model is the better for determining the statistical significance in answer to Parliament's question?

Rather than putting my dog in the fight, I'll let Met answer the question. The Met opens their response to Her Majesty's Government with the following:

*This
briefing paper has been produced to provide background information relating to
analyses undertaken in response to a series of Parliamentary Questions on the
use of statistical models to assess the global temperature record, and to address misleading ideas currently
appearing in various online discussions.*

The emboldened text above obviously refers to Keenan (he who shall not be named). But in a breathtakingly contradictory statement, Met concludes Keenan's critique of their choice of the AR(1) model holds merit:

*The
results show that the linear trend model with first-order autoregressive noise
[Met's model] is less likely to emulate the global surface temperature
timeseries than the driftless third-order autoregressive integrated model
[Keenan's model]. ... This provides some
evidence against the use of a linear trend model with first-order
autoregressive noise [Met's model] for
the purpose of emulating the statistical properties of instrumental records of
global average temperatures, as would be expected from physical
understanding of the climate system.*

Oops! Talk about an awkward moment. Then in a face saving move, the Met closes with this amazing statement:

*These
results have no bearing on our understanding of the climate system or of its
response to human influences such as greenhouse gas emissions and so the Met Office does not base its assessment of
climate change over the instrumental record on the use of these statistical
models.*

Translation: *"We don't
need no stinking statistical models to determine statistical significance. We
can determine it by eyeballing the data and you'll just have to take our word
at it*."

What Met has done was to state to Parliament that the rise in global temperature of 0.8 degrees Celsius since 1880 is statistically significant but Keenan is right in his criticism that the model they use is useless in making that determination. But no worries, they don't need statistical models to determine statistical significance.

Do you have any idea how incredibly wrong that sounds to a statistician? Statistical significance can only be determined by using statistical tests. That's why the term "statistically" is in front of "significant." That was the whole point of Parliament's question.

*[Edit June 14, 2013]: My work was brought to the attention of Lord Bernard Donoughue, UK Parliament Labor Peer, who spearheaded the investigation into
the MET Office inquiry. He offered a kind response to my article.
*

*This article is an updated reprint of my original work of the same title published at Suyts Space, an on-line science and politics venue.*