# Global warming since 1995 ‘now significant’

**From the Bottom of the Heap - R**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

**Updated: 17 July 2011**

Yesterday (June 11, 2011) the BBC reported on comments by Prof. Phil Jones, of the Climatic Research Unit (CRU), University of East Anglia (UEA), that the warming experienced by the planet since 1995 was statistically significant. That the trend in the recent data was now significant when last year it was not, was attributed to the fact that one more data point (year) was now available for the analysis and Jones highlighted the utility of long-term monitoring data as more data is collected in detecting pattern in noisy data. In this post I wanted to take a closer look at these data and illustrate how we can use R to fit regression models to time series data.

We will look at the HadCRUT3v data set, using the global means version of the data. I don’t know specifically that Jones was commenting on these particular data as I have been unable to identify a source for his comments other than the information in the BBC article, but it seems plausible that Jones was talking about these data. To load these data into R we use

```
URL <span class="o"><-</span> url<span class="p">(</span><span class="s">"http://www.cru.uea.ac.uk/cru/data/temperature/hadcrut3vgl.txt"</span><span class="p">)</span>
gtemp <span class="o"><-</span> read.table<span class="p">(</span>URL<span class="p">,</span> fill <span class="o">=</span> <span class="kc">TRUE</span><span class="p">)</span>
```

The data need a little post processing to clean them up as the provided format also has, in every other row, the number of stations used to form each temperature data point. We won’t consider those data in this post, but as they provide information as to the likely uncertainty in each observation, we could use those data as case weights in the analysis. The data come as monthly means, with an annual mean column at the end. The data are anomalies from the 1961-1990 mean.

```
<span class="c1">## Don't need the even rows --- perhaps do as case weights?</span>
gtemp <span class="o"><-</span> gtemp<span class="p">[</span><span class="o">-</span>seq<span class="p">(</span><span class="m">2</span><span class="p">,</span> nrow<span class="p">(</span>gtemp<span class="p">),</span> by <span class="o">=</span> <span class="m">2</span><span class="p">),</span> <span class="p">]</span>
<span class="c1">## set the Year as rownames</span>
rownames<span class="p">(</span>gtemp<span class="p">)</span> <span class="o"><-</span> gtemp<span class="p">[,</span><span class="m">1</span><span class="p">]</span>
<span class="c1">## Add colnames</span>
colnames<span class="p">(</span>gtemp<span class="p">)</span> <span class="o"><-</span> c<span class="p">(</span><span class="s">"Year"</span><span class="p">,</span> month.abb<span class="p">,</span> <span class="s">"Annual"</span><span class="p">)</span>
<span class="c1">## Data for 2011 incomplete so work only with 1850-2010 data series</span>
gtemp <span class="o"><-</span> gtemp<span class="p">[</span><span class="o">-</span>nrow<span class="p">(</span>gtemp<span class="p">),</span> <span class="p">]</span>
```

The first step in any analysis should be to look at the data. Here we produce a simple time series plot of the annual mean temperature using base graphics commands

```
ylab <span class="o"><-</span> expression<span class="p">(</span>Temperature<span class="o">~</span>Anomaly<span class="o">~</span><span class="p">(</span><span class="m">1961-1990</span><span class="p">)</span><span class="o">~</span>degree<span class="o">*</span>C<span class="p">)</span>
plot<span class="p">(</span>Annual <span class="o">~</span> Year<span class="p">,</span> data <span class="o">=</span> gtemp<span class="p">,</span> type <span class="o">=</span> <span class="s">"o"</span><span class="p">,</span> ylab <span class="o">=</span> ylab<span class="p">,</span>
main <span class="o">=</span> <span class="s">"Global mean temperature anomaly 1850-2010"</span><span class="p">)</span>
```

to produce this figure:

The news report was about trends in temperature from 1995 onwards. We subset the data, selecting only the observations from 1995–2010 (and only variables `Year`

and `Annual`

) and plot those data

```
grecent <span class="o"><-</span> subset<span class="p">(</span>gtemp<span class="p">,</span> subset <span class="o">=</span> Year <span class="o">>=</span> <span class="m">1995</span><span class="p">,</span>
select <span class="o">=</span> c<span class="p">(</span>Year<span class="p">,</span> Annual<span class="p">))</span>
<span class="c1">## plot</span>
plot<span class="p">(</span>Annual <span class="o">~</span> Year<span class="p">,</span> data <span class="o">=</span> grecent<span class="p">,</span> type <span class="o">=</span> <span class="s">"o"</span><span class="p">,</span> ylab <span class="o">=</span> ylab<span class="p">,</span>
main <span class="o">=</span> <span class="s">"Global mean temperature anomaly 1995-2010"</span><span class="p">)</span>
```

In R, linear models fitted by least squares are fitted using the `lm()`

function. In the code chunk below, I fit a linear trend to the 1995–2010 data (`gm1`

) and a null model of no change (`gm0`

), and compare the two models with a F-ratio test

```
<span class="c1">## fit a linear model through these recent data</span>
gm1 <span class="o"><-</span> lm<span class="p">(</span>Annual <span class="o">~</span> Year<span class="p">,</span> data <span class="o">=</span> grecent<span class="p">)</span>
gm0 <span class="o"><-</span> lm<span class="p">(</span>Annual <span class="o">~</span> <span class="m">1</span><span class="p">,</span> data <span class="o">=</span> grecent<span class="p">)</span>
anova<span class="p">(</span>gm0<span class="p">,</span> gm1<span class="p">)</span>
coef<span class="p">(</span>gm1<span class="p">)</span>
```

The linear trend fits the data significantly better than the no-trend model. The rate of increase in temperatures over the period is ~0.01 (standard error = 0.005). As these are time series data and probably violate the independence assumption of the model, the standard error is potentially too small. As the data are regularly spaced in time, we can easily employ the autocorrelation function to investigate residuals correlations in the model errors. The `acf()`

function can be used for that, which produces a plot of the correlogram

`acf<span class="p">(</span>resid<span class="p">(</span>gm1<span class="p">),</span> main <span class="o">=</span> <span class="s">"Residuals of linear model"</span><span class="p">)</span>`

Although within the approximate 95% confidence intervals, there are some large correlations in the residuals, especially at lag 2. This is something we must try to account for in our model and see how that affects the estimates of the model coefficients and ultimately the significance or otherwise of the fitted trend. Other model diagnostics can be produced via the `plot()`

methods for `“lm”`

objects

```
layout<span class="p">(</span>matrix<span class="p">(</span><span class="m">1</span><span class="o">:</span><span class="m">4</span><span class="p">,</span> ncol <span class="o">=</span> <span class="m">2</span><span class="p">))</span>
plot<span class="p">(</span>gm1<span class="p">)</span>
layout<span class="p">(</span><span class="m">1</span><span class="p">)</span>
```

but given the small sample size there is not much to be worried about here, though note the observations with large residuals and their leverage, which could indicate an undue influence on the model parameters.

We can account for the autocorrelation in the data whilst estimating the linear trend by switching to fitting the model using generalized least squares (GLS) via the `gls()`

function available in the **nlme** package. GLS allows a correlation structure for the model residuals to be estimated using a simple time series model such as an first-order auto-regressive process or AR(1). We need to decide what type of process to use for the correlation structure; the corellogram shown above suggests that AR terms, possibly up to order 2, might be appropriate. So we fit a series of models using `gls()`

that refits the null and linear trend models from before (just so we are confident we are comparing like with like), plus models using AR(1) and AR(2) processes for the residuals

```
require<span class="p">(</span>nlme<span class="p">)</span>
gg0 <span class="o"><-</span> gls<span class="p">(</span>Annual <span class="o">~</span> <span class="m">1</span><span class="p">,</span> data <span class="o">=</span> grecent<span class="p">,</span> method <span class="o">=</span> <span class="s">"ML"</span><span class="p">)</span>
gg1 <span class="o"><-</span> gls<span class="p">(</span>Annual <span class="o">~</span> Year<span class="p">,</span> data <span class="o">=</span> grecent<span class="p">,</span> method <span class="o">=</span> <span class="s">"ML"</span><span class="p">)</span>
gg2 <span class="o"><-</span> gls<span class="p">(</span>Annual <span class="o">~</span> Year<span class="p">,</span> data <span class="o">=</span> grecent<span class="p">,</span>
correlation <span class="o">=</span> corARMA<span class="p">(</span>form <span class="o">=</span> <span class="o">~</span> Year<span class="p">,</span> p <span class="o">=</span> <span class="m">1</span><span class="p">),</span> method <span class="o">=</span> <span class="s">"ML"</span><span class="p">)</span>
gg3 <span class="o"><-</span> gls<span class="p">(</span>Annual <span class="o">~</span> Year<span class="p">,</span> data <span class="o">=</span> grecent<span class="p">,</span>
correlation <span class="o">=</span> corARMA<span class="p">(</span>form <span class="o">=</span> <span class="o">~</span> Year<span class="p">,</span> p <span class="o">=</span> <span class="m">2</span><span class="p">),</span> method <span class="o">=</span> <span class="s">"ML"</span><span class="p">)</span>
```

Correlation structures are specified using one of the `“corStruct”`

classes; here I used the `corARMA()`

function to fit an ARMA process, but will only use AR terms. If all you want to fit is an AR(1) process in the residuals, then `corAR1(form ~ Year)`

can be used. All the models were fitted using maximum likelihood estimation. The `anova()`

method can be used to compare the sequence of nested models. By nested, I mean that you can go from the most complex model (`gg3`

) to the simplest model (`gg0`

) by setting successive parameters to zero. The output from this process is shown below

```
<span class="go">R> anova(gg0, gg1, gg2, gg3)</span>
<span class="go"> Model df AIC BIC logLik Test L.Ratio p-value</span>
<span class="go">gg0 1 2 -25.10920 -23.56403 14.55460 </span>
<span class="go">gg1 2 3 -28.01874 -25.70097 17.00937 1 vs 2 4.909536 0.0267</span>
<span class="go">gg2 3 4 -26.49875 -23.40839 17.24937 2 vs 3 0.480010 0.4884</span>
<span class="go">gg3 4 5 -30.77369 -26.91075 20.38685 3 vs 4 6.274945 0.0122</span>
```

As we already saw, the linear trend model fits the data significantly (*p*-value = 0.0267) better than the null model of no trend. The model with the AR(1) process does not significantly improve the fit over the linear model, whilst the model with the AR(2) processes provides a significantly better fit than the linear model. A direct comparison between the linear trend model and the linear trend plus AR(2) model indicates the substantially better fit of the latter, despite requiring the estimation of two additional parameters

```
<span class="go">R> anova(gg1, gg3)</span>
<span class="go"> Model df AIC BIC logLik Test L.Ratio p-value</span>
<span class="go">gg1 1 3 -28.01874 -25.70097 17.00937 </span>
<span class="go">gg3 2 5 -30.77369 -26.91075 20.38685 1 vs 2 6.754954 0.0341</span>
```

The residuals from the best model (`gg3`

) have much-reduced auto-correlations

`acf<span class="p">(</span>resid<span class="p">(</span>gg3<span class="p">,</span> type <span class="o">=</span> <span class="s">"normalized"</span><span class="p">))</span>`

Note that we need to use the normalized residuals to get residuals that take account of the fitted correlation structure. We are now in a position to assess the significance of the fitted trend within this model and investigate the slope of trend, which we can do using the `summary()`

method:

```
<span class="go">R> summary(gg3)</span>
<span class="go">Generalized least squares fit by maximum likelihood</span>
<span class="go"> Model: Annual ~ Year </span>
<span class="go"> Data: grecent </span>
<span class="go"> AIC BIC logLik</span>
<span class="go"> -30.77369 -26.91075 20.38685</span>
<span class="go">Correlation Structure: ARMA(2,0)</span>
<span class="go"> Formula: ~Year </span>
<span class="go"> Parameter estimate(s):</span>
<span class="go"> Phi1 Phi2 </span>
<span class="go"> 0.2412298 -0.6527874 </span>
<span class="go">Coefficients:</span>
<span class="go"> Value Std.Error t-value p-value</span>
<span class="go">(Intercept) -16.163962 6.232562 -2.593470 0.0212</span>
<span class="go">Year 0.008268 0.003112 2.656464 0.0188</span>
<span class="go"> Correlation: </span>
<span class="go"> (Intr)</span>
<span class="go">Year -1 </span>
<span class="go">Standardized residuals:</span>
<span class="go"> Min Q1 Med Q3 Max </span>
<span class="go">-2.29284316 -0.68980863 0.03087847 0.51005562 1.99216289 </span>
<span class="go">Residual standard error: 0.08715603 </span>
<span class="go">Degrees of freedom: 16 total; 14 residual</span>
```

The linear trend is still significant at the 95% level, although the estimated rate of change in temperature is a little lower then in the least squares model (0.008°C year^{-1}; standard error = 0.003). The estimates of the AR(2) parameters are also shown (0.24 and -0.65). Approximate confidence intervals on the estimated parameters can be produced using the `intervals()`

function:

```
<span class="go">R> intervals(gg3)</span>
<span class="go">Approximate 95% confidence intervals</span>
<span class="go"> Coefficients:</span>
<span class="go"> lower est. upper</span>
<span class="go">(Intercept) -29.531478957 -16.163962275 -2.79644559</span>
<span class="go">Year 0.001592535 0.008267935 0.01494333</span>
<span class="go">attr(,"label")</span>
<span class="go">[1] "Coefficients:"</span>
<span class="go"> Correlation structure:</span>
<span class="go"> lower est. upper</span>
<span class="go">Phi1 -0.1675377 0.2412298 0.38958957</span>
<span class="go">Phi2 -0.9036540 -0.6527874 -0.06838202</span>
<span class="go">attr(,"label")</span>
<span class="go">[1] "Correlation structure:"</span>
<span class="go"> Residual standard error:</span>
<span class="go"> lower est. upper </span>
<span class="go">0.05109686 0.08715603 0.14866224</span>
```

At this point, it would be useful to visualise the fitted trend on the observed data. To do this, we need to predict values from the fitted model for the trend over the period of interest. The `predict()`

method is used to derive predictions for new observations. First we build a data frame containing the values of `Year`

we want to predict at (1995–2010), and then we add a `yhat`

variable to the data frame, which contains the predicted values. The predicted values are obtained using `predict(gg3, newdata = pred)`

. The `transform()`

function is used to add the `yhat`

component to our data frame `pred`

:

```
<span class="go">R> pred <- data.frame(Year = 1995:2010)</span>
<span class="go">R> pred <- transform(pred, yhat = predict(gg3, newdata = pred))</span>
<span class="go">R> with(pred, yhat)</span>
<span class="go"> [1] 0.3305672 0.3388351 0.3471030 0.3553710 0.3636389 0.3719069 0.3801748</span>
<span class="go"> [8] 0.3884427 0.3967107 0.4049786 0.4132465 0.4215145 0.4297824 0.4380503</span>
<span class="go">[15] 0.4463183 0.4545862</span>
<span class="go">attr(,"label")</span>
<span class="go">[1] "Predicted values"</span>
```

A final step, having produced the predicted values is to plot the trend on the original data series, which we do using the code below, first with the full data series and then with the 1995–2010 subset

```
layout<span class="p">(</span>matrix<span class="p">(</span><span class="m">1</span><span class="o">:</span><span class="m">2</span><span class="p">,</span> ncol <span class="o">=</span> <span class="m">2</span><span class="p">))</span>
<span class="c1">## plot full data</span>
plot<span class="p">(</span>Annual <span class="o">~</span> Year<span class="p">,</span> data <span class="o">=</span> gtemp<span class="p">,</span> type <span class="o">=</span> <span class="s">"o"</span><span class="p">,</span> ylab <span class="o">=</span> ylab<span class="p">,</span>
main <span class="o">=</span> <span class="s">"Global mean temperature anomaly 1850-2010"</span><span class="p">)</span>
lines<span class="p">(</span>yhat <span class="o">~</span> Year<span class="p">,</span> data <span class="o">=</span> pred<span class="p">,</span> lwd <span class="o">=</span> <span class="m">2</span><span class="p">,</span> col <span class="o">=</span> <span class="s">"red"</span><span class="p">)</span>
<span class="c1">## plot the 1995-2010 close-up</span>
plot<span class="p">(</span>Annual <span class="o">~</span> Year<span class="p">,</span> data <span class="o">=</span> grecent<span class="p">,</span> type <span class="o">=</span> <span class="s">"o"</span><span class="p">,</span> ylab <span class="o">=</span> ylab<span class="p">,</span>
main <span class="o">=</span> <span class="s">"Global mean temperature anomaly 1995-2010"</span><span class="p">)</span>
lines<span class="p">(</span>yhat <span class="o">~</span> Year<span class="p">,</span> data <span class="o">=</span> pred<span class="p">,</span> lwd <span class="o">=</span> <span class="m">2</span><span class="p">,</span> col <span class="o">=</span> <span class="s">"red"</span><span class="p">)</span>
layout<span class="p">(</span><span class="m">1</span><span class="p">)</span>
```

to produce this figure

What does this tell us about underlying mean temperature of the recent few years? The linear trend doesn’t fit the data all that well; in fact, one could argue for a flat trend the latter half of the 1995–2010 period, and when viewed in the context of the data leading up to 1995, the fitted trend appears to substantially *underestimate* the increase in global mean temperature. A non-linear trend migth do a better job.

By focussing on only the very recent observational period, we have neglected to consider the constantly evolving mean temperature over the past 160 years. It has always seemed to me somewhat perverse that climatologists and climate sceptics alike ignored the decades of data before the recent period and then argued the toss about the sign and strength of a linear trend in the recent period, often claiming that there is insufficient data to ascribe a statistically significant pattern in defence of their particular point of view.

To me, it seems far more sensible to fit a model to the entire series of data, but to use a regression technique that can model local features of the series rather than the global trend over the entire data. That model could then be used to ask questions about whether temperatures have increased in a statistically significant manner for whatever period one liked. In a follow-up posting, I’ll demonstrate one such modelling technique that builds upon the concepts introduced here of using regression models with auto-correlation structures for the residuals.

##
Update 17 July 2011

In trying to keep the original posting simple, I overlooked one important aspect of the model fitting; to get the best estimates of the correlation parameters we should be using REML, not ML. This was raised in the comments, and I promised to take a look and update the post. This is important as it does change the output of the models, and if you are the sort of person who gets hung up on p-values that change will have you pulling your hair out or grinning like the proverbial Chesire cat, depending on your point of view on climate change…

Load the data as per the original posting above and fit the trend and no trend models using GLS. For this we **do** need to fit with ML:

```
require<span class="p">(</span>nlme<span class="p">)</span>
gg0 <span class="o"><-</span> gls<span class="p">(</span>Annual <span class="o">~</span> <span class="m">1</span><span class="p">,</span> data <span class="o">=</span> grecent<span class="p">,</span> method <span class="o">=</span> <span class="s">"ML"</span><span class="p">)</span>
gg1 <span class="o"><-</span> gls<span class="p">(</span>Annual <span class="o">~</span> Year<span class="p">,</span> data <span class="o">=</span> grecent<span class="p">,</span> method <span class="o">=</span> <span class="s">"ML"</span><span class="p">)</span>
anova<span class="p">(</span>gg0<span class="p">,</span> gg1<span class="p">)</span>
```

Now, update `gg1`

so the fitting is done via REML and fit models with AR(1) and AR(2) residuals (the `method = “REML”`

bits are redundant as this is the default, but are included to make it clear that we want REML fitting):

```
gg1 <span class="o"><-</span> update<span class="p">(</span>gg1<span class="p">,</span> method <span class="o">=</span> <span class="s">"REML"</span><span class="p">)</span>
gg2 <span class="o"><-</span> gls<span class="p">(</span>Annual <span class="o">~</span> Year<span class="p">,</span> data <span class="o">=</span> grecent<span class="p">,</span>
correlation <span class="o">=</span> corARMA<span class="p">(</span>form <span class="o">=</span> <span class="o">~</span> Year<span class="p">,</span> p <span class="o">=</span> <span class="m">1</span><span class="p">),</span> method <span class="o">=</span> <span class="s">"REML"</span><span class="p">)</span>
gg3 <span class="o"><-</span> gls<span class="p">(</span>Annual <span class="o">~</span> Year<span class="p">,</span> data <span class="o">=</span> grecent<span class="p">,</span>
correlation <span class="o">=</span> corARMA<span class="p">(</span>form <span class="o">=</span> <span class="o">~</span> Year<span class="p">,</span> p <span class="o">=</span> <span class="m">2</span><span class="p">),</span> method <span class="o">=</span> <span class="s">"REML"</span><span class="p">)</span>
```

A quick call or two to the `anova()`

method suggests that the models with AR(1) or AR(2) residuals probably are not justified/required

```
<span class="gp">> </span>anova<span class="p">(</span>gg1<span class="p">,</span> gg2<span class="p">,</span> gg3<span class="p">)</span>
<span class="go"> Model df AIC BIC logLik Test L.Ratio p-value</span>
<span class="go">gg1 1 3 -13.29388 -11.37671 9.646942 </span>
<span class="go">gg2 2 4 -12.66126 -10.10503 10.330630 1 vs 2 1.367377 0.2423</span>
<span class="go">gg3 3 5 -14.31948 -11.12419 12.159740 2 vs 3 3.658220 0.0558</span>
<span class="gp">> </span>anova<span class="p">(</span>gg1<span class="p">,</span> gg3<span class="p">)</span>
<span class="go"> Model df AIC BIC logLik Test L.Ratio p-value</span>
<span class="go">gg1 1 3 -13.29388 -11.37671 9.646942 </span>
<span class="go">gg3 2 5 -14.31948 -11.12419 12.159740 1 vs 2 5.025597 0.081</span>
```

If we were to stick with model `gg1`

then the trend is, marginally significant, with zero not included in a 95% confidence interval on the estimate slope of the fitted trend

```
<span class="gp">> </span>confint<span class="p">(</span>gg1<span class="p">)</span>
<span class="go"> 2.5 % 97.5 %</span>
<span class="go">(Intercept) -40.521627038 -2.48493179</span>
<span class="go">Year 0.001433605 0.02042816</span>
```

Part of the reason that the AR(1) model is no longer significant is that the estimate of the autocorrelation parameter is very uncertain; there is so little data with which to estimate this parameter any more precisely (output truncated):

```
<span class="gp">> </span>intervals<span class="p">(</span>gg2<span class="p">,</span> which <span class="o">=</span> <span class="s">"var-cov"</span><span class="p">)</span>
<span class="go">Approximate 95% confidence intervals</span>
<span class="go"> Correlation structure:</span>
<span class="go"> lower est. upper</span>
<span class="go">Phi -0.2923927 0.3280283 0.7541095</span>
<span class="go">attr(,"label")</span>
<span class="go">[1] "Correlation structure:"</span>
<span class="go">....</span>
```

and given this uncertainty, we might wish to include the AR(1) to allow for autocorrelation in the residuals even if the estimate we get is not significantly different from zero. At which point you’d be looking at an *in*significant trend estimate

```
<span class="gp">> </span>confint<span class="p">(</span>gg2<span class="p">)</span>
<span class="go"> 2.5 % 97.5 %</span>
<span class="go">(Intercept) -48.116090814 3.46065310</span>
<span class="go">Year -0.001535917 0.02422018</span>
```

This exercise just reinforces the futility of trying to identify significant trends, or lack thereof, with such a small sample of data. Using the full data set, an additive modelling approach and the use of derivatives suggests that temperatures were significantly increasing well after 2000. See this post for details.

**leave a comment**for the author, please follow the link and comment on their blog:

**From the Bottom of the Heap - R**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.