Anomaly Regression – Do It Right!

I have been meaning to do this post for several years, but until now I have not found a particularly relevant time to do it.  In his recent post at the Air Vent, Jeff Id makes the following statement:

Think about that.  We’re re-aligning the anomaly series with each other to remove the steps.  If we use raw data  (assuming up-sloping data), the steps in this case were positive with respect to trend, sometimes the steps can be negative.  If we use anomaly alone (assuming up-sloping data), the steps from added and removed series are always toward a reduction in actual trend.  It’s an odd concept, but the key is that they are NOT TRUE trend as the true trend, in this simple case, is of course 0.12C/Decade.

The actual situation is deeper than Jeff thinks.  The usual method used by climate scientists for doing monthly anomaly regression is wrong! Before you say, “Whoa! How can a consensus be wrong?”, let me first give an example which I will follow up with the math to show you what the problem is.

We first produce a series of ten years of ten years worth of noiseless temperature data.  To make it look realistic (not really important) we will take a sinusoidal annual curve and superimpose an exact linear trend (of .2 degrees per year) on the data:

We now calculate the anomalies and fit a linear regression line to the anomalies.  The following is a plot of anomalies with a regression line:

The resulting trend of .198 does not match the actual trend in the data despite the fact that there is no noise present in the series.  Although the difference in this case is not large (a poor oft-repeated justification), this is an uncorrectable bias error.  What is not as obvious is that this will have an impact on the autocorrelation in the situation as evidenced by the ACF of the residuals – an effect due to the methodology and not inherent to the actual “error” sequence:

But it gets worse. If we change our starting month (but still use ten full years of data), we get different slopes for each month:

Month	Intercept	Trend
1	-1.27983	0.198014
2	-1.28521	0.198931
3	-1.28971	0.199681
4	-1.29331	0.200264
5	-1.29595	0.200681
6	-1.2976	0.200931
7	-1.29822	0.201014
8	-1.29775	0.200931
9	-1.29618	0.200681
10	-1.29344	0.200264
11	-1.2895	0.199681
12	-1.28431	0.198931

The average is equal to .2 so all is not lost.   However, all of the errors are a result of how we chose to do our analysis and were not really present in the original data.

So what’s causing the problem?

Here comes the math

What most climate scientists still do not seem to understand is the need to spell out the parameters and how they relate to each other.  Unless such a “statistical model” exists, we do not have a good grasp of what is happening when we carry out an analysis and a reduced ability to recognize when something is not quite right.  The implicit model in this situation is

where

X(t) = temperature at time t

t = y + m, where y = year and m = month (with values 0/12 , 1/12, …,11/12)

µ_m = mean of month m

β = annual trend

ε(t) = “error” (which in this case are all zeroes).

Now, what happens when we calculate anomalies?  If A(t) is the anomaly at time t,

where the “bar” represents averaging over that variable.

Here,

and if we substitute this in the original equation we get

The important thing to realize here is the the month no longer appears with the trend – only the year (and NOT time) should be used to get the correct trend.  By using time, we actually are fitting the line

for which the intercept should be different for each month.  Since the usual anomaly regression fits a single intercept, the resulting trend is incorrectly estimated.

How can this be fixed?

Several fixes are possible.  The simplest is to use year rather than time as the independent variable in the regression.  This may seem counterintuitive, but it takes the anomalies and lines then up vertically above each other solving the problem that Jeff noticed.  However, this solution can still be problematical if the number of available anomalies differs from year to year.

A better solution is to do the anomalizing and the trend fit at the same time.  This corresponds to a single factor Analysis of Covariance where month is treated as a categorical factor with twelve levels and time (or year!) is treated as a numeric covariate.  This can be implemented in R using the function lm.

The script used in the post follows:

#we construct 11 years of data for later use
seasonal = 5 + 10*sin(pi*(0:131)/6)
#add trend of .2 degrees per year
trend = (0:131)/60

temps = ts((seasonal+trend)[1:120],start=c(1,1),freq=12)
time.temp= time(temps)

plot(seasonal,type="l")

plot(temps, main="Plot of Monthly Temperature Time Series", ylab="Degrees C", xlab="Time") #anomregfig1.jpg

#function to calculate anomaly
anomaly.calc=function(tsdat,styr=1951,enyr=1980){
 tsmeans = rowMeans(matrix(window(tsdat,start=c(styr,1),end=c(enyr,12)),nrow=12))
 tsdat-rep(tsmeans,len=length(tsdat))}

anom = anomaly.calc(temps,1,10)

#anomaly regression
reg1 = lm(anom~time.temp) #  intercept =  -1.180, slope = 0.198

plot(anom, main = "Anomalies",ylab = "Degree C",xlab = "Time") #anomregfig2.jpg
abline(reg = reg1,col="red")

#calculate autocorrelation of residuals
acf(residuals(reg1))  #anomregfig3.jpg

#set up revolving regressions
tempsx = ts(seasonal+trend,start=c(1,1),freq=12)
timex.temp= time(tempsx)
anomx.all = anomaly.calc(tempsx,1,11)

#function to cycle regression starting points and give coefficients
cycle.reg = function(dats) {
 coes = matrix(NA,12,2)
 coes[1,] = coef(lm(window(dats,start=c(1,1),end=c(10,12))~time(window(dats,start=c(1,1),end=c(10,12)))))
 for (i in 1:12) coes[i,] = coef(lm(window(dats,start=c(1,i),end=c(11,i-1))~time(window(dats,start=c(1,i),end=c(11,i-1)))))
coes}

(all.coef = cycle.reg(anomx.all))

#           [,1]      [,2]
# [1,] -1.279832 0.1980138
# [2,] -1.285205 0.1989305
# [3,] -1.289710 0.1996805
# [4,] -1.293305 0.2002639
# [5,] -1.295949 0.2006806
# [6,] -1.297599 0.2009306
# [7,] -1.298215 0.2010140
# [8,] -1.297754 0.2009306
# [9,] -1.296176 0.2006806
#[10,] -1.293437 0.2002639
#[11,] -1.289497 0.1996805
#[12,] -1.284314 0.1989305

#fit regression using year
year = floor(time.temp)
reg2 = lm(anom~year)
reg2 #[1]  -1.1          0.2

#setup anova
month = factor(rep(1:12,11))
time2 = time(anomx.all)
year=floor(time2)
dataf = data.frame(anomx.all,time2,year,month)

#function to test different month starting points
#using time as a covariate
#column 13 is trend

cyclex.reg = function(datfs) {
 coes = matrix(NA,12,13)
 for (i in 1:12) {dats = datfs[i:(119+i),]
 coes[i,] = coef(lm(anomx.all~0+month+time2,data=dats))} #line corrected
coes}

(anov1 = cyclex.reg(dataf))

#function to test different month starting points
#using year as covariate
cyclexx.reg = function(datfs) {
 coes = matrix(NA,12,13)
 for (i in 1:12) {dats = datfs[i:(119+i),]
 coes[i,] = coef(lm(anomx.all~0+month+year,data=dats))} #line corrected
coes}

(anov2 = cyclexx.reg(dataf))