Ron had asked me to consider adding measurement uncertainty to the observed values in the graph earlier this week. It occurred to me that I could estimate the uncertainty due to measurement as follows:
- Peruse papers to discover the portion of the uncertainty in the annual average temperature measurement due to lack of station coverage. I found values for GISTemp in Hansen et al. 2010 and for NOAA in Brohan et. al. These were 0.05C for the (2s) value and 0.03 for the 1s value respectively. For Hadley, I read off the ±uncertainty from the annual time series.
- I transformed these into standard errors (i.e. 1s values), and then, assuming this error is white, I obtained monthly standard errors by multiplying the error for the annual value by 121/2. I then generated 120 months of noise with the monthly standard error, added the appropriate uncertainty to its corresponding to the observations , and computed trends for each error-padded observation. I repeated 1000 times then took the standard deviation of the of the error-padded trends for each observation. These were $latex \displaystyle s_o $ =( 0.028C/dec, 0.034C/dec, 0.033 C/dec) for GISS, HadCrut and NOAA respectively. For all future calculations, I took the largest of these to use as the uncertainty in measurements: that is I used $latex \displaystyle s_o = 0.034C/dec $
- I then estimate of the standard error for the reported value of an observation about the model mean to include the following three uncertainties: a) the uncertainty we would expect for “model weather” (i.e. $latex \displaystyle \sigma_m $), b) the uncertainty in our estimate of the model mean based on the average of ‘n’ runs $latex \displaystyle \frac { \sigma} {n} $, and c) the measurement uncertainty $latex \displaystyle s_o $. This resulted in a total standard error of:
(1) $latex \displaystyle \sigma_T^2 =\sigma_m^2 \frac{1+n} {n} + s_o^2 $So, to be clear, this standard error is supposed to encompass the region in which a trend computed from a single realization of earth data measured with the amount of error claimed by HadCrut- would fall if a model produces the correct trend and variance for weather. (Some caveats apply. for example: we really don’t know (i.e. $latex \displaystyle \sigma_m $) with 100% certainty. Accounting for the uncertainty would expand the error bars. I just don’t know how much.)
As before, I treated the standard error for model weather as known, and multiplied by the critical “z” value for the ±95% confidence intervals. I superimposed these on the previous values computed without including the uncertainty due to measurement error, $latex \displaystyle s_o $. The graph with a start date of 2001 and the new wider uncertainty intervals added on top of the old narrower ones is shown below:
If you squint, you can see there really, really are two sets of horizontal hashes on the error bars.
Now for a bleg: I don’t like showing both the narrower and wider error bars in the same color. Tomorrow, I’m flipping through a book with examples of plots with error bars. Any tips on making plots with error bars would be welcome. I would particularly like verbose examples.
You can use the polygon function. Let x0 and x1 be the bounds on the x-coordinate. For your particular graph above, these should be less than zero and more than 10 to ensure the error band completely covers the horizontal plot region. Let y0 and y1 be your lower and upper limits of uncertainty. Then the command polygon(x=c(x0,x1,x1,x0),y=c(y0,y1,y1,y0),col=’grey’) should produce a grey band. I’ve used this to show the 1-sigma and 2-sigma uncertaintly. To do that, plot the 2-sigma band first, then the 1-sigma band. Here’s an example:
rm(list = ls())
# Open a new plot window, define its limits and create a bounding box.
plot.new()
plot.window(xlim = c(0,100), ylim = c(-4,4))
box()
# Plot the 2-sigma confidence interval.
polygon(x=c(-5,105,105,-5), y=c(-2,-2,2,2), col = ‘lightgrey’)
# Plot the 1-sigma confidence interval.
polygon(x=c(-5,105,105,-5), y=c(-1,-1,1,1), col = ‘darkgrey’)
# Plot 100 random points N(0,sigma)
points(rnorm(100))
Chad–Thanks. I’ll use that for some future graphs. It’s great.
But for this graph, it’s not a good method because people can’t usually look at things, square two separate uncertainty intervals and take the square root in their head. So, showing the measurement ones as the grey bands and the model ones as uncertianty intervals gives the impression that the total uncertainty would be found like this:
When the correct relationship is:
(1) $latex \displaystyle \sigma_T^2 =\sigma_m^2 \frac{1+n} {n} + s_o^2 $
Compare the results when the ratio of
(1) $latex \displaystyle \frac{s_o} {\sigma_m} << 1 $
and you’ll see the problem. What happens is we make a pretty picture that misleads the person looking at the graph.
if you squint real hard you’ll see that the blue uncertainty intervals are really a slightly shorter one just inside a slightly longer one.Both are blue. What I would like to do is either
a) Make the horizontal bars on the longer uncertainty intervals larger or
b)