Monday, October 12, 2009

Global Cooling?

While I don't appreciate most of the pseudo-religious global warming hysterics, I can't say the recent BBC piece on 'global cooling' is anything but sensationalism. It's curious that they don't just show their readers a simple graph of the data.

I don't know about you, but the recent years that they are claiming show a cooling trend look well within the expected variability. The solid line on the graph is a simple linear regression, the inner dashed lines are a 0.95-level confidence interval and the outer dashed lines are a 0.95-level prediction interval. The prediction interval gives a range of where you would expect temperature measurements to fall based on the fitted trend and residuals. You can download the NASA temperature data yourself and have a look.

The R code to make this graph is shown below

temps = read.csv("temps.csv", header=TRUE)
attach(temps)

# fit two models, one for January to December annual averages, and one
# for December to November annual averages
m.1 = lm(JD ~ Year)
m.2 = lm(DN ~ Year)

# Make confidence and prediction intervals
m.1.cinterval = predict(m.1, level=0.95, interval="confidence")
m.2.cinterval = predict(m.2, level=0.95, interval="confidence")
m.1.pinterval = predict(m.1, level=0.95, interval="prediction")
m.2.pinterval = predict(m.2, level=0.95, interval="prediction")

# plot the data and the fits/intervals to a png file
png("temps.png", width=640, height=480)

plot(Year, JD, ylab="Degrees C")
lines(Year, m.1.cinterval[,1], lty=1)
lines(Year, m.1.cinterval[,2], lty=2)
lines(Year, m.1.cinterval[,3], lty=2)
lines(Year, m.1.pinterval[,2], lty=2)
lines(Year, m.1.pinterval[,3], lty=2)
title("NASA GISS Global Average Temperature")

dev.off()

In case you're interested here's the summary of the fitted model shown in the graph. Just looking at the graph shows there's probably something more complicated than just a linear trend going on though.

Call:
lm(formula = JD ~ Year)

Residuals:
Min 1Q Median 3Q Max
-0.317448 -0.093632 0.001046 0.087362 0.298468

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.9770910 0.5788968 5.143 1.01e-06 ***
Year 0.0056581 0.0002977 19.009 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1244 on 126 degrees of freedom
Multiple R-squared: 0.7414, Adjusted R-squared: 0.7394
F-statistic: 361.3 on 1 and 126 DF, p-value: < 2.2e-16

To be fair, the BBC piece does go into greater detail that it isn't really global cooling, that the Pacific Decadal Oscillation might be starting to kick in, so we'll probably see a little 20 year dip like that between 1940 and 1960. I'd say you need a few more years of data to actually be able to call it though (kind of like calling an economic recession).

1 comment:

  1. Just read an interesting article that does a trend analysis (with R of course) using several data sets rather than just the GISS data.

    What I like about it is that it uses openly accessible data and processes it in a transparent way.

    ReplyDelete