copyright © 2007 Paolo B. DePetrillo, MD
Let's see some Surface Temperature Analysis Station Data from Haskell (33.2N, 99.8W) The raw data can be downloaded here in a text file, and imported in Excel.
Haskell
At the far right we have two columns. "ANN" and "MISSANN". I didn't make up the first label. The "ANN" column contains mean temperature data in Celsius for the year. You will also notice a "MISSANN" column, which is the same data except for the fact that I removed data points from "ANN" which relied on extrapolated data. "MISSANN" is clean pure, well almost pure, data. What does this mean? It means that on occasion, a data point is missing so investigators estimate the value by extrapolation and the data is treated the same way as the measured values. I don't know the procedure for how these values were derived, since it is not referenced. I would rather place my faith in the law of large numbers rather than in my ability to estimate some value for a missing data point. This kind of study has literally millions of data points given all the locations. If there is a signal, we should find it regardless of a few missing points.
Here is the CO2 data. I have made the file local so as not to burden their server. co2tableWRJ
Examples of Mathematical Modeling of Local Temperature Times Series
First thing I do before anything else is look at graphics of data. It's pretty. It helps to correct typographical errors. And If you get really lucky, you don't have to run statistics. You'll get your answer by observation and common sense. Alas and alack, the universe throws in enough noise to make that impossible except for very well designed experiments, crafted by geniuses. Not here.
Here are the graphs of the data. The A column shows a simple mean fit and the B column shows a linear fit. The first row uses extrapolated data, the second row uses the raw data with extrapolated average temperature points removed.
A B
Temperature = 17.53 +/- 0.07 Temperature = 0.0047 +/- 0.0021 (Year) + 8.27 +/- 4.22
SSE = 63.08 SSE = 60.38
Temperature = 17.53 +/- 0.08 Temperature = 0.0036 +/- 0.0026 (Year) + 10.45 +/- 5.12
SSE = 42.89 SSE = 41.87
Does the better fit justify adding another parameter to the model?
Let's compare models with the corrected Akaike's Information Criteria
First, for the two using extrapolated data:
Compare models with F test
Model SS DF
Mean Fit Using Extrapol (null) 63.08 109
Linear Using Ext (alternative) 60.38 108
Difference 2.7 1
Percentage Difference 4.47% 0.93%
Ratio (F) 4.83
P value 0.0301
If Mean Fit Using Extrapol (the null hypothesis) were true, there would be a 3.01% chance of obtaining results that fit Linear Using Ext (the alternative hypothesis) so well.
Since the P value is less than the traditional significance level of 5%, you can conclude that the data fit significantly better to Linear Using Ext than to Mean Fit Using Extrapol.
And now let's see the same models using raw data
Compare models with F test
Model SS DF
Mean Fit Using Raw (null) 42.89 80
Linear Using Raw (alternative) 41.87 79
Difference 1.02 1
Percentage Difference 2.44% 1.27%
Ratio (F) 1.92
P value 0.1693
If Mean Fit Using Raw (the null hypothesis) were true, there would be a 16.93% chance of obtaining results that fit Linear Using Raw (the alternative hypothesis) so well.
Since the P value is not less than the traditional significance level of 5%, you can conclude that the data do not fit significantly better to Linear Using Raw than to Mean Fit Using Raw. Since there is no evidence to reject the null hypothesis, the preferred model is Mean Fit Using Raw.
Thanks to the folks at GraphPad for their nice Java AIC calculator.
Digression
I've looked at thousands of time series, so sometimes I see stuff that may not be obvious. For example, do you see in both graphs that there is both a fast temperature cycle over a couple of years, and a slower cycle lasting somewhere between 60 and 70 years? I'll draw it out for you. The red line goes through the busiest part of the graphs. To illustrate my point, I did a 3 point running average [cringe!] to smooth out some of the peaks and valleys in the temperature graphs above.
One of the simplest models for a time series is the average of all the points in the time series. In model language, we can express that as:
Temperature = P1 + error
where P1 is the average and error is an error term.
I noted previously that it looked like there may be a cycle expressed in the point cloud. To test this hypothesis, we have to figure out some kind of relationship between time, expressed in years, and temperature. An approximation is found in the trigonometric functions, and we will use the Cosine function that gives us a nice cyclical pattern. With any luck, it might fit the data better than a straight line. At the end of this section, I try to explain what makes least-mean square regression dangerous in this kind of analysis, and why I chose this method.
The Cosine function takes a few parameters. The technique varies the parameters associated with a Jordan matrix using the Gauss-Newton method, looking for positive and negative peaks in the response surface that minimize the error of the fit according to an objective function. In plain English, what we are doing is trying to find the value of the parameters in the function that minimize the error. These techniques offer an efficient way to do it, rather than looking at all possible combinations. That might take a while. Here is what the function might look like:
Temperature = P1 x (cosine(Year + P2) x 3.14159 / P3) + P4
The cosine function will vary from (-1) to +(1) but we have to give it the variable in radians. That is why we multiply by pi. We have the date given in years "Year" but we can't assume that the cosine function starts at whatever place has an output of 1! We have to offset it with a parameter, "P2" so that we vary the phase given what is an arbitrary set of increasing integers. Not quite arbitrary but we could have started at 2022 or whatever. What this number does is shift the wavy cosine function to the left or the right to help with the overall fit.
We have to divide by 1/2 the cycle length so that the function reaches a maximum and minimum within that span of time. That is where P3 comes into the picture. P4 is a mean, a point to which we add or subtract the amount given by the cosine expression to get our answer.. And let's not forget about P1, which tells us the amplitude of the cosine function.
Here is what it looks like when we unleash Gauss on the function. It seems to follow that point cloud. It's not a great fit, but it's better than a line. Lets see what kind of parameters we discovered. The error of the estimate, plus or minus, is in the parentheses. The lower, the better.
Temperature = 0.35[0.11] x cosine((Year -330.74[224.50] x 3.14 ) / 32.46[3.20] + 17.56[0.08]
About every 32 years there may be a temperature oscillation. Chances are pretty good that this has been noted in the literature, and a search indicates that in the Northern hemisphere, we are subject to the "Pacific Decadal Oscillation" that has a 1/2 cycle, warm or cold, around every 30 years.
http://jisao.washington.edu/pdo/
This is good, because like I said, the model should reflect reality at some point.
Now I know that there is also an ENSO cycle, La Nino or La Nina events, that occurs every couple of years. You can read all about it here:
http://www.cpc.ncep.noaa.gov/products/analysis_monitoring/ensocycle/enso_cycle.shtml
If we add another term to the function, one that has a 1/2 cycle of every few years, we may get a better fit. Actually, we can't "add" another term, but we can suggest starting at some point related to what we know about the world - which in this case is that ENSO periods are about 4 years long. In this case, choosing an initial value of "1" also works, it seems there is a strong reduction of error in the function at a 1/2 cycle of around 2 years.
Temperature = 0.35[0.11] x cosine(Year - 369.24[222.22]) / 33.02[3.17] ) - 0.28[0.10] x cosine(Year-493.31[11.43] / 2.03[0.02] + 17.55[0.08]
The low error with the 1/2 cycle of 2.03 years, and with the phase shift error of 11.43 gives me some confidence that this is modeling something real, like the ENSO, which is more regular than the Pacific Decadal Oscillation, hence the higher error on the phase shift as well as the cycling time of 33 years. That last term, 17.55 degrees Celsius, is what the average temperature would be in absence of any cyclical events.
Sunspots!
Here is a graph of the sunspot cycle.
And here is a graph of the temperature cycle and the sunspot cycle superimposed. I have taken the logarithm of the sunpot number and substracted 15 from the average temperature so they would both fit on the same scale for purposes of comparison.
To my eye, there seems to be a pattern where the sunspot number begins to fall from a maximum, followed by the temperature, the reverse occurring as well. The rate of change in sunspot number also seems to determine the magnitude of the temperature change. Here is one of the better models relating temperature and sunspot number. Here, we model a lagged sunspot number coupled to a sunspot cycle, which hopefully will reduce the error enough to justify entering it into the model. The ugly Sunspot term is simply the area under the Sunspot - Time curve using the trapezoidal rule. So here it is:
Full Model
T = + P1 x Cosine {(Year+ P2) *3.1415) / P3 "PDO" term
+ P4 x Cosine {(Year +P5) *3.1415) / P6 "ENSO" term
+ P7 x {Cosine {(Year +P8) *3.1415) / P9} Sunspot Cycle
* (LogSun,lag 0 years + 2*LogSun, lag 1 years + 2*LogSun, lag 2 years +2*LogSun, lag 3 years+2*LogSun, lag 4 years + LogSun, lag 5 years )"Sunspot number" term
+ P10 baseline temperature
value [standard error]
P1 = 0.4 [ 0.1]
P2 = -1209 [ 298 ]
P3 = 31 [ 3 ] "PDO" term cycle length
P4 = 0.3 [ 0.1 ]
P5 = -285[ 17]
P6 = 2.16 [ 0.02 ] "ENSO" term cycle length
P7 = -0.009 [ 0.003 ] * {Sunspot Number Term)
P8= -625 [41]
P9 = 4.47 [ 0.07 ] Sunspot cycle length
P10 = 17.58 [ 0.08]
T = predicted mean annual temperature
SSE 30.70 DFE 71 MSE 0.43 RMSE 0.6630.70
Compare models with the corrected Akaike's Information Criteria
Linear Fit Full Model
Sum-of-squares 41.87 30.70
Number of data points 81 81
Number of parameters 2 10
Akaike's Information Criteria (corrected, AICc) -47.14 -52.76
Probability model is correct 5.68% 94.32%
Difference in AICc 5.62
Information ratio 16.61
Full Model has a lower AICc than Linear Fit so is more likely to be the correct model.
It is 16.6 times more likely to be correct than Linear Fit.
Compare models with F test
Model SS DF
Linear Fit (null) 41.87 79
Full Model (alternative) 30.70 71
Difference 11.17 8
Percentage Difference 36.38% 11.27%
Ratio (F) 3.23
P value 0.0035
If Linear Fit (the null hypothesis) were true, there would be a 0.35% chance of obtaining results that fit Full Model (the alternative hypothesis) so well.
Since the P value is less than the traditional significance level of 5%, you can conclude that the data fit significantly better to Full Model than to Linear Fit.
Thanks to the nice folks at GraphPad for the model comparison calculations.
Finally, we ask if the Full Model is better or worse than the Full Model which includes a linear trend. In other words, if we remove the effects of the cyclical climate variations and sunspots on temperature by taking them into account with the Full Model,is there any unexplained variation that could be due to a linear trend?
In other words, we want to compare our full cyclical model to our full cyclical model with an added linear trend.
Full Model + Linear Trend
T = + P1 x Cosine {(Year+ P2) *3.1415) / P3 "PDO" term
+ P4 x Cosine {(Year +P5) *3.1415) / P6 "ENSO" term
+ P7 x {Cosine {(Year +P8) *3.1415) / P9} Sunspot Cycle
* (LogSun,lag 0 years + 2*LogSun, lag 1 years + 2*LogSun, lag 2 years +2*LogSun, lag 3 years+2*LogSun, lag 4 years + LogSun, lag 5 years )"Sunspot number" term
+ P10 baseline temperature
+P11*Year Linear Trend
I am not going to give the resulting parameters but for this model, but here is the "fit" data:
SS = 30.62 DFE = 70 MSE = 0.43 RMSE = 0.66
And here is the comparison:
Compare models with the corrected Akaike's Information Criteria
Full Model Full Model + linear Trend
Sum-of-squares 30.69 30.62
Number of data points 81 81
Number of parameters 10 11
Akaike's Information Criteria (corrected, AICc) -52.79 -50.21
Probability model is correct 78.39% 21.61%
Difference in AICc 2.58
Information ratio 3.63
Full Model has a lower AICc than Full Model + linear Trend so is more likely to be the correct model.
It is 3.6 times more likely to be correct than Full Model + linear Trend.
Compare models with F test
Model SS DF
Full Model (null) 30.69 71
Full Model + linear Trend (alternative) 30.62 70
Difference 0.07 1
Percentage Difference 0.23% 1.43%
Ratio (F) 0.16
P value 0.6904
If Full Model (the null hypothesis) were true, there would be a 69.04% chance of obtaining results that fit Full Model + linear Trend (the alternative hypothesis) so well.
Since the P value is not less than the traditional significance level of 5%, you can conclude that the data do not fit significantly better to Full Model + linear Trend than to Full Model. Since there is no evidence to reject the null hypothesis, the preferred model is Full Model.
Thanks to the nice folks at GraphPad for the model comparison calculations.
What about is carbon dioxide? Well, it does not enter into the model at this point because the parameter is small { 0.004 } and the error is high { 0.005 }, and the overall fit is worse. The carbon dioxide signal in this model is small, and the model does not have the resolution to answer the question at this time.
In this case, the simplest model was the linear fit model with two parameters. We need to make sure we have not simply added parameters and generated a fit, while losing all the model information. In other words, we are penalized for adding additional parameters. The error has to decrease by enough to justify the added parameters.
SSE = 44.30, which is WORSE than a simple mean fit.
The CO2 data was obtained from Hawaii. Previous to 1957, it was from ice cores. It is expressed in PPM. In this extremely simplistic model, there is no relationship between the increase in CO2 and temperature. Additional models would involve extra parameters and may not be explanatory. I always like to compare them to the simplest case.
There appear to be warm and cold cycles of the same degree (pun intended) both at early and later points in time. Since the carbon dioxide curve does not look cyclical, it is probable that those variations may be due to other factors.
Conclusions
For Haskell, TX, neither the raw or extrapolated data set robustly support an increase in temperature over time.
In Haskell, TX, since 1890, at least three factors may influence the mean annual temperature.
A cycle with a period of about 64 years, ?PDO
A cycle with a period of about 4 years ?ENSO
A cycle with a period of about 11 years (sunspot cycle) and the effect of average sunspot number .
The best fit is obtained when the number of sunpsots value is taken into account as the area under the sunspot - time curve lagged up to 5 years. This makes sense as we would expect a large system like climate to require some time before experiencing the effects of the solar cycle.
Sunspot cycles probably alter temperature, and should probably be tested in GCM models. OK, this is a stretch, but if other locations show the same pattern, it strengthens the argument.
The model predicts a gradual decrease in average temperature for the next 20-25 years in Haskell, TX, based on entering the "cold" half of the PDO cycle, as well as the cold part of the sunspot cycle. We'll see.
Discussion
You need to always compare your models to the simplest case, otherwise you run the risk of adding parameters that are not justified by the better fit. I want to emphasize this point because you will come across "least-mean square" linear fits that seem to "fit" data well in the form Y = P1X + P2 where P1 is the slope and P2 is the intercept. Unless these two parameter models are compared to the simplest case, namely the one parameter "mean fit" model of the form Y = P2, any conclusion about "goodness of fit" is suspect.
The linear model using extrapolated data suggests there appears to be a little bit of a warming trend in Haskell, TX since 1889. About 0.005 +/- 0.002 degrees per year. However, if the same models are run using raw data, there is no statistical support for choosing the more complex model, and the warming trend disappears. The average yearly temperature in Haskell, TX may be 17.53 C +/- 0.08 and a warming trend may not be evident. Take your pick. I do not like using extrapolated data, nor doing linear fit models without comparing the linear fit to the simplest case. So your mileage may vary.
I say "may not be" because just like we accept a confidence level of 5% to reject the null hypothesis that there is no warming, meaning that 5% of the time we reject the null hypothesis even though it is true, we also need have a confidence level for accepting the null hypothesis, though it is rarely mentioned in reports of many studies. In this instance, stating "there is no warming" without also giving a confidence level for this statement is just plain stupid!.]
Much of the biological scientific literature has been polluted by this unclear thought, with statements of the form " There is no difference in the means between group A and Group B" without giving a confidence level for the statement. It is just as illogical as stating "There is a difference in the means between group A and Group B" without giving a confidence level. The confidence level for accepting the null hypothesis is "beta." I'm not going into the specifics yet, but my confidence level for accepting the null hypothesis is also 5%, meaning that 5% of the time I will accept it even though it is not true. In the case of effects of carbon dioxide on warming, I cannot accept the null hypothesis that there is no effect, since the "beta" is > 5%. It is still an open question which, for data I have analyzed, can only be answered by more data. More data, better conclusions. Isn't that always the case?
In the past century, we have experienced 2 x 1/2 cycles of PDO " warm," and a little more than one 1/2 cycle of PDO "cold." Look at the graph above. It trends cold from 1890 to 1900 but part of the cycle is not present. It again trends cold from 1948 to around 1979. But there are two complete warm trend cycles - one from 1900 to 1938 and one from 1978 to present. This skews the data towards a warm result over this period of time if this is not accounted by the model. I wish we had a couple of centuries worth of CO2 and local temperature data, since this would help refine the model by including a number of these cycles. As time permits, I will perform the same analysis for several locations in the United States and present the results in these pages.
Data is from 1890 to present. This limits the conclusions to this time period, namely the magnitude of many of the parameters.
Data is incomplete, since no extrapolated points were used when data points were missing. This may have biased the results.
I have not ruled out an effect of carbon dioxide on temperature. In stat speak, I do not have the power to rule it out at some reasonable confidence level. And I don't reach a reasonable level of significance to rule it in. Hence, the question is not answered. Like everything else, it would have been nice to have a lot more data.
There are obviously other factors influencing the temperature, and of course it is possible that carbon dioxide may itself influence the length or phase of some of the cycles. That is an empiric question which needs expert study.
Some time I'm going to get to Haskell, TX. And for all the work done by dedicated folks at GISS and NASA, a profound sense of awe and thanks.
I will also pretty up this page at some point.
Next page