Demonstration of Improved Anomaly Calculation

This post is an application and further explanation of  Roman’s work on improved anomaly trend analysis.  At this post, Roman gave an example of a perfect sine wave with a trend of 0.2C added to it.  By taking the monthly anomaly of the data, steps are introduced on an annual basis.   These steps cause a distortion of trend which are actually unnecessary, however, there is a greater understanding and a little bit of a tricky regression involved(in my opinion-not Roman’s) .  The programming of the regression is absolutely simple in R where the difference between the standard call first and the improved call second are below.

coes = coef(lm(anomalydata~ time))

coes = coef(lm(rawdatat~0+month+time))    #Improved calculation

The first line returns two coefficients for the anomaly regressed against time (standard form) y = ax + b where a(slope) and b(intercept) are returned.   The second example line above, also returns coefficients but there are actually 13 returned.  The function is again y= ax+b, however 12 b’s are returned for a single a.  One  best fit line is created for each month with the only difference between the lines being the offset!

In the second line of code above, “month” is a special type of indicator variable which instructs R to find the offsets separately from the slope.  It’s a sort of mask which separates the data by month for the intercept but not by slope.  If you haven’t seen it before, it probably will take some study to understand.  What it means though is simple, the function returns 12 offsets based on each month, yet only one slope which is based on the best simultaneous fit to all the months together.

Of course the best slope is dependent on the best offsets by month so solving together is necessary.

There is a code advantage to the improved call though, in that you no longer need to calculate the anomaly.   Figure 1, shows the fit to GHCN station number 60360.

Figure 1 - 12 offset least squares slope calculation

See we got the slope without ever actually calculating the anomaly and using a single line in R.  What’s more is, we’ve removed the anomaly step problem created by monthly averaging.  There are some interesting implications of this method as well.  Since the simple monthly average anomaly decreases the quality of the fit (increases the residuals), the improved method trend uncertainty, is narrowed but more interestingly is that there is a “true” anomaly value which can be extracted from this method.  The true anomaly, takes into account the trend information by month such that monthly steps created by the simple removal of monthly mean and the simple linear trend is identical to the improved calculation results.

There are dozens of demonstrations which could be done on this topic, I’ve used a single temperature station, but it might be fun to try with simple linear data or some simple waveforms.  Anyway, that is left for the reader to figure out.  Below, I’ve written an appropriate slope calculation function which returns the 12+1 coefficients, residuals and true anomaly from a temperature series as shown above.

######################################################################
##### slope plotting function with modifications for 12 month anomaly
######################################################################
get.slope=function(dat=dat)
{
	len=length(dat)
	month = factor(rep(1:12,as.integer(len/12)+1))
	month=month[1:len]	#first month by this method is not necessarily jan
					#the important bit is that there are 12 steps
	time2=time(dat)
	coes = coef(lm(dat~0+month+time2))
	residual=ts(rep(NA,len),start=time2[1],deltat=1/12)
	fitline=array(NA,dim=c(len,12))
	for(i in 1:12)
	{
		fitline[,i]=coes[13]*time2+coes[i]
	}
	anomaly=rep(NA,len)
	for(i in 1:len)
	{
		residual[i]=dat[i]-fitline[i,month[i]]
	}
	anomaly = residual+time2*coes[13]-mean(time2*coes[13],na.rm=TRUE)
	list(coes,residual,anomaly)
}

So, employing the station data above, the anomaly calculation using the standard method results in a plot which looks like this:

Figure 2 - Standard anomaly calculation

And the anomaly returned by the above improvement is shown in Fig 3:

Figure 3, Improved "True anomaly" calculation

Totally and completely different  — as you can see ;).   The next plot is a little weird, I subtracted the anomaly from Fig 2 from the residual of Fig 3.  The reason I  didn’t just subtract anomaly is that you get a huge up and down cyclic scribble which gives no feel for the magnitude of the  difference, engineers need to understand the magnitude.  Simple version – Think of this plot as the difference between Fig 2 minus Fig 3 with the slope of Fig 3 artificially added in for scale.  The wiggle itself is the error in the standard anomaly calculation corrected for by using the simple  improved line of code presented at the top of this post.

Figure 4 - Annual signal error in standard anomaly calculation.

In this case the increased trend from 0.2253 to 0.2256 C/Decade – or 0.0003 C/Decade.  Not much!  The confidence interval tightened from 0.0818 to 0.0817 C/Decade. Again, not terribly much.

There are of course examples which can be created that will have bigger differences between the methods.  Now that I’ve completed this though, I’ve been able to create complete global gridded trend without direct calculation of anomaly at any point.  I’ll post that next.