Monday, September 22, 2014

Trend since 1998—significant??

I had a question sent to me about the trend since 1998.  As many of you know, my last post included an analysis which showed that the linear regression trend since 1998 was statistically significant.

Trends versus start year.  Error bars are the 95% confidence intervals.
My questioner asked if I had accounted for autocorrelation in my analysis.  The short answer is "No, I did not."  The reason?  According to my analysis, it wasn't necessary.

Here are my methods and R code.

#Get coverage-corrected HadCRUT4 data and rename the first two columns
CW<-read.table("http://www-users.york.ac.uk/~kdc3/papers/coverage2013/had4_krig_annual_v2_0_0.txt", header=F)
names(CW)[1]<-"Year"
names(CW)[2]<-"Temp"

#Analysis for autocorrelation—I check manually as well but so far the auto.arima function has performed admirably.
library(forecast)
auto.arima(resid(lm(Temp~Year, data=CW, subset=Year>=1998)), ic=c("bic"))

The surprising result?
Series: resid(lm(Temp ~ Year, data = CW, subset = Year >= 1998))
ARIMA(0,0,0) with zero mean    

sigma^2 estimated as 0.005996:  log likelihood=18.23
AIC=-34.46   AICc=-34.18   BIC=-33.69
 I was expecting something on the order of ARIMA(1,0,1), which is the autocorrelation model for the monthly averages.  Taking the yearly average rather than the monthly average effectively removed autocorrelation from the temperature data, allowing the use of a white-noise regression model.

trend.98<-lm(Temp~Year, data=CW, subset=Year>=1998)
summary(trend.98)
Call:
lm(formula = Temp ~ Year, data = CW, subset = Year >= 1998)

Residuals:
     Min       1Q           Median        3Q          Max
-0.14007  -0.05058   0.01590    0.05696    0.11085

Coefficients:
                    Estimate       Std. Error    t value    Pr(>|t|) 
(Intercept)   -19.405126   9.003395    -2.155     0.0490 *
Year             0.009922      0.004489     2.210     0.0443 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.08278 on 14 degrees of freedom
Multiple R-squared: 0.2587,    Adjusted R-squared: 0.2057
F-statistic: 4.885 on 1 and 14 DF,  p-value: 0.04425
The other surprise?  That the trend since 1998 was significant even with a white-noise model.  Sixteen data points is not normally enough to reach statistical significance unless a trend is very strong.

Temperature trend since 1998

Sunday, September 21, 2014

The "no warming" claim rises from the dead yet again.

Like a movie vampire, this one keeps coming back no matter how many stakes are driven through its heart.  I've covered this one (here, here, and here).  Bluntly: There is absolutely no evidence that global warming has stopped.  For global warming to stop, the Earth's energy balance must be either zero or negative.  The most recent estimates for the energy imbalance are generally between +0.5 W/m2 and +1.0 W/m2.  The only way the Earth is not going to warm while it is gaining energy is if the laws of thermodynamics magically do not apply.  If the Earth is gaining energy, some part of it, somewhere, must be getting warmer.  The heat must go into either melting ice, warming the oceans, warming the land, or warming the atmosphere (or some combination thereof).

Tuesday, September 16, 2014

WUWT and how NOT to test the relationship between CO2 and temperature

WUWT published a piece by Danle Wolfe which purports to measure the correlation between CO2 and global temperature.  As you can probably predict, Wolfe's conclusion is that there is no relationship.
"Focusing on the most recent hiatus below, both visually and in a 1st order linear regression analysis there clearly is effectively zero correlation between CO2 levels and global mean temperature."
 Unfortunately for Wolfe, all he's produced is a fine example of mathturbation as well as an example of forming a conclusion first then warping the evidence to fit.

Thursday, September 11, 2014

James Taylor versus relative humidity and specific humidity

It appears that the relative humidity and specific humidity continues to trip some people up.  Yes, I'm thinking of the screed James Taylor wrote on Forbes.com on Aug. 20.  In his article, Taylor trumpets two "facts".  First, that relative humidity has declined and second, that specific humidity isn't rising as fast as global climate models predict.  Since climate models assume that relative humidity has stayed constant, Taylor then claims that models are overestimating global warming.  Unfortunately for Taylor, his "facts" don't check out.

Monday, September 8, 2014

R code for my Seasonal Trends graph

I had a request for the code I used to generate the graphs in my Seasonal Trends post.


While it looks complex, the R code for it is very simple IF you have the data ready.    I'm assuming that you already have the temperature dataset you want as an R object (I have several datasets in an object I simply call "monthly": GISS, Berkeley Earth, Cowtan-Way, HadCRUT4, UAH, and RSS, along with the year/decimal month Time, Year, and numeric Month).  The code I used to generate the graph is as follows:
#Create separate datasets for each season

monthly.79<-subset(monthly, Time>=1978.92 & Time<2013.91)
DJF<-subset(monthly.79, Month=="12" | Month =="1" | Month=="2")
DJF$Year_2 <- numeric (length (DJF$Year))
for (i in 1:length (DJF$Year) ) {
        if ( DJF$Month [i] == 12) {
                DJF$Year_2[i] <-   DJF$Year [i] + 1
        }
        else {
                DJF$Year_2[i] <-   DJF$Year [i]
        }
}
MAM<-subset(monthly.79, Month=="3" | Month =="4" | Month=="5")
JJA<-subset(monthly.79, Month=="6" | Month =="7" | Month=="8")
SON<-subset(monthly.79, Month=="9" | Month=="10" | Month=="11")

#Calculate the seasonal average for each year

DJF<-aggregate(DJF$BEST, by=list(DJF$Year_2), FUN=mean)
MAM<-aggregate(MAM$BEST, by=list(MAM$Year), FUN=mean)
JJA<-aggregate(JJA$BEST, by=list(JJA$Year), FUN=mean)
SON<-aggregate(SON$BEST, by=list(SON$Year), FUN=mean)

#Check for autoregression

library(forecast) #for the auto.arima function

auto.arima(resid(lm(x~Group.1, data=DJF)), ic=c("bic"))

auto.arima(resid(lm(x~Group.1, data=MAM)), ic=c("bic"))

auto.arima(resid(lm(x~Group.1, data=JJA)), ic=c("bic"))

auto.arima(resid(lm(x~Group.1, data=SON)), ic=c("bic"))

 #Construct the plot

plot(x~Group.1, data=DJF, type="l", col="blue", xlab="Year", ylab="Temperature anomaly (ÂșC)", main="Seasonal Climate Trends", ylim=c(-0.1, 0.8))
points(x~Group.1, data=MAM, type="l", col="green")
points(x~Group.1, data=JJA, type="l", col="red")
points(x~Group.1, data=SON, type="l", col="orange")

#Add the trend lines

abline(lm(x~Group.1, data=DJF), col="blue", lwd=2)
abline(lm(x~Group.1, data=MAM), col="green", lwd=2)
abline(lm(x~Group.1, data=JJA), col="red", lwd=2)
abline(lm(x~Group.1, data=SON), col="orange", lwd=2)

legend(1979, 0.8, c("DJF", "MAM", "JJA", "SON"), col=c("blue", "green", "red", "orange"), lwd=2)
#Get the slopes

summary(lm(x~Group.1, data=DJF)
summary(lm(x~Group.1, data=MAM)
summary(lm(x~Group.1, data=JJA)
summary(lm(x~Group.1, data=SON)
 That's all there was to it.  I just repeated this code, modifying only the period of the initial subset to create the second graph (Monthly.2002<-subset(monthly, Time>=2001.92 & Time<2013.91) and the related seasonal subsets.  To the person who requested my code: Hope this helps.

Monday, September 1, 2014

One hundred years ago today...

...the last passenger pigeon, a female named Martha, died in the Cincinnati Zoo.  A species that once had an estimated population size of 3 billion was destroyed in roughly 50 years by a combination of habitat loss and overhunting.  The story of that extinction is being told in numerous articles on this centenary (i.e. in Nature, Audubon Magazine, National Geographic) and at museums like the Smithsonian Institute which tell the story far better than I could here.  The Audubon Magazine article, in particular, is well worth reading as it details the history of the extinction.