The burning question for me therefore has been: If the stick guys' methods don't make sense and perhaps introduce bias, then what do we get when proper methods are used? What, exactly, is the truth of the matter?
And so I've worked through Steve's postings and the writings of others trying to find out what a proper method might be.
In a post from 18 months ago, Steve pointed out that the job of processing this data and eliminating the effects of various factors to come up with an estimation of the temperature anomalies over time could be done using well understood mixed effects statistical modeling algorithms. I thought it rather astounding that a statistical method seemingly completely different from Hansen's methods would come up with much the same results. Steve used this in modeling the combination of location data into grid cells, but why not apply it to the whole data set?
And so I tried that.
giss.dset0, you may recall is temperature data taken from different weather stations. In some cases these stations are all together in one location, as in 4 or 5 stations in a given city area. Hansen first combined these stations together for a location, then he did a lot of complicated smoothing and adjustment of the data in a series of steps. dset0 is raw data, dset1 includes some adjustments and combination of data of stations at given locations, dset2 includes more adjustments within a given location, and gridding of the data does more of the same, smoothing data over geometrical areas, between locations. With all of the adjusting and smoothing going on the introduction of some sort of bias seems likely. For example, Hansen always used the longest records to start with and then increased or decreased the shorter records to match the longer. But aren't the longer records usually from urban stations? Would this not cement the heat island effect into the combined location data?
And so, working in the R Project statistical package and following McIntyre's lead, I first reorganized the data in giss.dset0 into a simple long data frame with location information folded in from giss.info.
extr = function(giss.dset0, giss.info) {
year = 0; temp = 0; alt=0; ur=0; loc=0; lxs=0; sta=0
Ni = (1:length(giss.dset0))[!is.na(giss.dset0)]
for(i in Ni) #if(sum(is.na(giss.dset0[[i]])) > 0) next else
for(j in 1:length(giss.dset0[[i]])) if(is.na(giss.dset0[[i]][[j]])) next else
{
N = length(giss.dset0[[i]][[j]][[1]])
sta = c(sta, rep(names(giss.dset0[[i]])[j],N))
year = c(year, giss.dset0[[i]][[j]][[1]])
temp = c(temp, giss.dset0[[i]][[j]][[18]])
ur = c(ur, rep(giss.info$urban[i],N))
alt = c(alt, rep(250*(giss.info$alt.interp[i]%/%250),N))
tmp=substr(names(giss.dset0[[i]])[j],1,8)
loc = c(loc, rep(tmp, N))
lxs = c(lxs, rep(paste(tmp, "i", names(giss.dset0[[i]])[j],sep=""),N))
}
sta=sta[-1]; year=year[-1]; temp=temp[-1]; alt=alt[-1]; loc=loc[-1]
ur = ur[-1]; lxs=lxs[-1]
return(data.frame(loc=loc,sta=sta,year=year,ur=ur,alt=alt,lxs=lxs,temp=temp))
}
dset0.t = extr(giss.dset0, giss.info)
(No, I could not use monthly data without running out of memory.)
Then calculate random effects:
(load the lme4 package)
dset0.fm = lmer(temp~(1|loc)+(1|sta)+(1|ur)+(1|alt)+(1|year), data=dset0.t)
dset0.fm1 = lmer(temp~(1|sta)+(1|year), data=dset0.t)
(loc = location id, sta=station id with location, ur=urban vs suburban vs rural, alt= location altitude in 250 meter bands)
ranef(dset0.fm)$year and ranef(dset0.fm1)$year contain the random effects year over year. The temperature anomaly for each year, in other words. Most of the variance is dumped into other groups. Using more groupings in the model formula makes no difference in these estimates (dset0.fm vs dset0.fm1 R-squared = 0.999).
Here is your heat island effect:
> ranef(dset0.fm)$ur
(Intercept)
R -0.23716400 (rural)
S -0.04312743 (suburban)
U 0.28028908 (urban)
And here is the effect of altitude:
> ranef(dset0.fm)$alt
(Intercept)
0 10.0821297
250 9.5952519
500 8.2160810
750 7.4049196
1000 6.4774994
1250 5.4415439
1500 4.2931860
1750 3.4352673
2000 1.7974976
2250 -0.5588156
2500 -1.6522236
2750 -5.5811832
3000 -5.2383055
3250 -13.6373919
3500 -8.7024097
3750 -13.8671371
4000 1.2833311
4250 -6.0728011
4500 -3.3266110
4750 0.8495323
5000 -1.0044707
5750 -0.4843320
What that yearly random effect looks like compared to the GISSTEMP book numbers:
And here is the difference between the two:
R-squared = 0.73. Close enough for government work. The slope coefficients are nearly identical (a 0.00035 per year difference). And I'm sure that a statistician would do much better. Much of what Hansen did in step 1 and 2 and in gridding the data was smoothing. Application of a smoothing algorithm would no doubt take out the variation in the random effects estimation I got, but that would not improve the fidelity of the results.
dset0 annual averages are as close to the raw data as I can get without running out of computer memory. I went from that data to a final result with one command in R.
If I can come up with the same results with no spacial factors at all other than altitude and station grouping then does gridding mean anything? Can the analysis represent anything more than the 25% of the earth's surface the stations occupy?
Hansen is vindicated in a big respect. This shows that his methods, regardless how byzantine, introduce no significant amount of bias where temperature trends are concerned with regard to the station data available. The proof is in the pudding, and all that.
That lack of coverage of the earth's surface in station data, I'm guessing, is the reason for the difference in satellite, radiosonde, and station estimates of global temperature.
The amount of memory and processing power I have in my cheap, clearance priced desktop PC would have been duplicated by a computer the size of a city block, if at all, back in the days Hansen started doing this climate business. His multi step process of analyzing station data was perhaps an effort to perform an analysis that would have been impossible using the statistical methods I used here.
Update
Including another grouping factor, latitude in 5 degree bands, has no effect on the temperature anomalies, but it's interesting to look at the random effect for that group:
dset0.fm2 = lmer(temp~(1|ur)+(1|lat)+(1|alt)+(1|sta)+(1|year),
data=dset0.t)
> ranef(dset0.fm2)$lat
(Intercept)
-90 -48.0765255
-85 -30.6208086
-80 -30.5843329
-75 -26.8087508
-70 -17.5837928
-65 -10.6445832
-55 -0.4677458
-50 1.9302840
-45 4.6332613
-40 7.5261542
-35 10.1896612
-30 7.2221295
-25 16.5283583
-20 17.4532507
-15 18.3289638
-10 18.4701555
-5 18.5784181
0 18.7600821
5 19.5108772
10 19.8686941
15 19.7051168
20 17.5278219
25 15.0410194
30 11.2047241
35 7.3644670
40 3.2930061
45 -0.4509212
50 -0.3098546
55 -5.3064547
60 -9.0101544
65 -13.3917199
70 -16.3916299
75 -20.7265894
80 -22.7631131
And including that group gives more consistent altitude and urban estimates:
> ranef(dset0.fm2)$alt
(Intercept)
-1000 7.6009074
0 8.7884377
250 8.2316191
500 5.1136268
750 9.5794809
1000 5.9054092
1250 5.3426933
1500 4.0768528
1750 3.4410112
2000 2.0043036
2250 1.0247369
2500 -0.5251244
2750 -1.8542720
3000 -0.6135464
3250 -8.0035847
3500 -3.3222669
3750 -5.1468110
4000 -5.1510483
4250 -8.3572983
4500 -8.9318326
4750 -5.5925205
5000 -5.9627467
5750 -7.6483137
> ranef(dset0.fm2)$ur
(Intercept)
R -0.89551186
S -0.03941906
U 0.93493366
Here is the summary from lmer():
Linear mixed model fit by REML
Formula: temp ~ (1 | loc) + (1 | sta) + (1 | ur) + (1 | alt) + (1 | year)
Data: dset0.t
AIC BIC logLik deviance REMLdev
1274202 1274280 -637094 1274190 1274188
Random effects:
Groups Name Variance Std.Dev.
sta (Intercept) 0.510624 0.71458
loc (Intercept) 113.129532 10.63624
year (Intercept) 0.141931 0.37674
alt (Intercept) 29.856394 5.46410
ur (Intercept) 0.068179 0.26111
Residual 0.521113 0.72188
Number of obs: 548500, groups: sta, 13483; loc, 4514; year, 127; alt, 23; ur, 3
Fixed effects:
Estimate Std. Error t value
(Intercept) 5.20 1.25 4.159
Linear mixed model fit by REML
Formula: temp ~ (1 | ur) + (1 | lat) + (1 | alt) + (1 | sta) + (1 | year)
Data: dset0.t
AIC BIC logLik deviance REMLdev
1288347 1288426 -644167 1288337 1288333
Random effects:
Groups Name Variance Std.Dev.
sta (Intercept) 11.67732 3.41721
year (Intercept) 0.14162 0.37633
lat (Intercept) 330.71739 18.18564
alt (Intercept) 38.58698 6.21184
ur (Intercept) 0.84184 0.91752
Residual 0.52133 0.72203
Number of obs: 548500, groups: sta, 13483; year, 127; lat, 34; alt, 23; ur, 3
Fixed effects:
Estimate Std. Error t value
(Intercept) -1.544 3.427 -0.451
Update II
I get almost exactly the same results with the model used in the first update when I do the analysis with giss.dset2 (R-squared of temp random effects = 0.97)