James A. Rising

The America’s Water project, coordinated at Columbia’s Water Center by Upmanu Lall, is trying to understand the US water system as an integrated whole, and understand how that system will evolve over the next decades. Doing so will require a comprehensive model, incorporating agriculture, energy, cities, policy, and more.

We are just beginning to lay the foundation for that model. A first step is to create a network of links between station gauges around the US, representing upstream and downstream flows and counties served. The ultimate form of that model will rely on physical flow data, but I created a first pass using simple rules:

1. Every gauge can only be connected to one downstream gauge (but not visa versa).
2. Upstream gauges must be at a higher elevation than downstream gauges.
3. Upstream gauges must be fed by a smaller drainage basin than downstream gauges.
4. Of the gauges that satisfy the first two constraints, the chosen downstream gauge is the one with the shortest distance and the most “plausible” streamflow.

The full description is available on Overleaf. I’ve applied the algorithm to the GAGES II database from USGSU, which includes all station gauges with at least 20 years of data.

Every red dot is a gauge, black lines are upstream-downstream connections between gauges, and the blue and green lines connect counties with each of the gauges by similar rules to the ones above (green edges if the link is forced to be longer than 100 km).

This kind of network opens the door for a lot of interesting analyses. For example, if agricultural withdrawals increase in the midwest, how much less water will be available downstream? We’re working now to construct a full optimization model that accounts for upstream dependencies.

Another simple question is, how much of the demand in each county is satisfied by flows available to it? Here are the results, and many cities show up in sharp red, showing that their demands exceed the surface water by 10 times or more.

I put this here as a warning. While growing degree-days (GDD) are well-known as an effective model to predict yields, they don’t perform so hot at the country-scale.

I used mean temperature GDDs, between 8 and 24 degrees C, estimated at many locations from station data, and then using the weighted average by production within each country. Here are the results:

As you can see, for most crops, these GDDs aren’t even significant, and as frequently negative as positive. This defies a century of agricultural research, but the same data at a fine spatial scale seems to work just fine.

There’s a lot we don’t know about future climates. Despite continuous improvements, our models still do notoriously poorly with some fairly fundamental dynamics: ocean energy transport, the Indian monsoon, and just reproducing past climate. Worse, we have the computing resources to run each model a hand-full of times. We haven’t estimated the distribution of possible temperatures that even one model predicts.

Instead, we treat the range of models as a kind of distribution over beliefs. This is useful, but a poor approximation to an actual distribution. For one, most of the models are incestuously related. Climate is chaotic, no matter how you model it, so changing the inputs within measurement error can produce a whole range of future temperatures.

This is all motivation for a big innovation from the Risky Business report, on the part of Bob Kopp and DJ Rasmussen.

First, they collected the full range of 44 CMIP5 models (downscaled and adjusted to predict every US county). They extracted each one’s global 2080-2099 temperature prediction.

Then they ran an *hemisphere* scale model (MAGICC), for which its possible to compute enough runs to produce a whole distribution of future temperatures.

Here’s the result: a matching between models and the temperatures distribution for end-of-century global temperatures.

[From American Climate Prospectus, Technical Appendix I: Physical Climate Projections]

The opaque maps are each the output of a different GCM. The ghosted maps aren’t– at least, not alone. In fact, every ghosted map represents a portion of the probability distribution which is currently missing from our models.

To fix this, they made “surrogate” models, by decomposing existing GCM outputs into “forced” seasonal patterns which scale with temperature, and unforced variability. Then they scaled up the forced pattern, added back in the variability, and bingo! A new surrogate model.

Surrogates can’t tell us the true horrors of a high-temperature world, but without them we’d be navigating without a tail.

R doesn’t have a built-in method for calculating heteroskedastically-robust and clustered standard errors. Since this is a common need for us, I’ve put together a little library to provide these functions.

Download the library.

The file contains three functions:

• get.coef.robust: Calculate heteroskedastically-robust standard errors
• get.coef.clust: Calculate clustered standard errors
• get.conf.clust: Calculate confidence intervals for clustered standard errors

The arguments for the functions are below, and then an example is at the end.

get.coef.robust(mod, var=c(), estimator=”HC”)
Calculate heteroskedastically-robust standard errors

mod: the result of a call to lm()
     e.g., lm(satell ~ width + factor(color))
var: a list of indexes to extract only certain coefficients (default: all)
     e.g., 1:2 [to drop FE from above]
estimator: an estimator type passed to vcovHC (default: White's estimator)
Returns an object of type coeftest, with estimate, std. err, t and p values

get.coef.clust(mod, cluster, var=c())
Calculate clustered standard errors

mod: the result of a call to lm()
     e.g., lm(satell ~ width + factor(color))
cluster: a cluster variable, with a length equal to the observation count
     e.g., color
var: a list of indexes to extract only certain coefficients (default: all)
     e.g., 1:2 [to drop FE from above]
Returns an object of type coeftest, with estimate, std. err, t and p values

get.conf.clust(mod, cluster, xx, alpha, var=c())
Calculate confidence intervals for clustered standard errors

mod: the result of a call to lm()
     e.g., lm(satell ~ width + factor(color))
cluster: a cluster variable, with a length equal to the observation count
     e.g., color
xx: new values for each of the variables (only needs to be as long as 'var')
     e.g., seq(0, 50, length.out=100)
alpha: the level of confidence, as an error rate
     e.g., .05 [for 95% confidence intervals]
var: a list of indexes to extract only certain coefficients (default: all)
     e.g., 1:2 [to drop FE from above]
Returns a data.frame of yhat, lo, and hi

Example in Stata and R

The following example compare the results from Stata and R for an example from section 3.3.2 of “An Introduction to Categorical Data Analysis” by Alan Agresti. The data describes 173 female horseshoe crabs, and the number of male “satellites” that live near them.

Download the data as a CSV file. The columns are the crabs color (2-5), spine condition (1-3), carapace width (cm), number of satellite crabs, and weight (kg).

The code below calculates a simple model with crab color fixed-effects, relating width (which is thought to be an explanatory variable) and number of satellites. The first model is a straight OLS regression; then we add robust errors, and finally crab color clustered errors.

The analysis in Stata.

First, let’s do the analysis in Stata, to make sure that we can reproduce it in R.

import delim data.csv

xi: reg satell width i.color
xi: reg satell width i.color, vce(hc2)
xi: reg satell width i.color, vce(cluster color)

Here are excerpts of the results.

The OLS model:

------------------------------------------------------------------------------
      satell |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       width |   .4633951   .1117381     4.15   0.000     .2428033    .6839868
       _cons |  -8.412887   3.133074    -2.69   0.008    -14.59816   -2.227618
------------------------------------------------------------------------------

The robust errors:

------------------------------------------------------------------------------
             |             Robust HC2
      satell |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       width |   .4633951   .1028225     4.51   0.000     .2604045    .6663856
       _cons |  -8.412887   2.939015    -2.86   0.005    -14.21505   -2.610726
------------------------------------------------------------------------------

Note that we’ve used HC2 robust errors, to provide a direct comparison with the R results.

The clustered errors.

------------------------------------------------------------------------------
             |               Robust
      satell |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       width |   .4633951   .0466564     9.93   0.002     .3149135    .6118767
       _cons |  -8.412887   1.258169    -6.69   0.007    -12.41694   -4.408833
------------------------------------------------------------------------------

The robust and clustered errors are a little tighter for this example.

The analysis in R.

Here’s the equivalent code in R.

source("tools_reg.R")
tbl <- read.csv("data.csv")

mod <- lm(satell ~ width + factor(color), data=tbl)

summary(mod)
get.coef.robust(mod, estimator="HC2")
get.coef.clust(mod, tbl$color)

And excerpts of the results. The OLS model:

               Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -8.4129     3.1331  -2.685  0.00798 ** 
width            0.4634     0.1117   4.147 5.33e-05 ***

The coefficient estimates will be the same for all three models. Stata gives us -8.412887 + .4633951 w and R gives us -8.4129 + 0.4634 w. These are identical, but with different reporting precision, as are the standard errors.

The robust errors:

(Intercept)    -8.41289    2.93902 -2.8625  0.004739 ** 
width           0.46340    0.10282  4.5067 1.229e-05 ***

Again, we use HC2 robust errors. Stata provides 2.939015 and .1028225 for the errors, matching these exactly.

The clustered errors:

                Estimate Std. Error  t value  Pr(>|t|)    
(Intercept)    -8.412887   1.258169  -6.6866 3.235e-10 ***
width           0.463395   0.046656   9.9321 < 2.2e-16 ***

Stata gives the exact same values as R.

Sol, Amir, and I have been slaving away over a report on the business-case for fighting climate change.  And it was released this morning!  The media outlets give a sense of the highlights:

Take a look!  Here’s the media report:
http://riskybusiness.org/report/overview/executive-summary

And the scientific report (what we actually helped write):
http://rhg.com/reports/climate-prospectus

I’m in a group studying Bayesian techniques for environmental research, under Upmanu Lall at the Earth Institute. We recently did a session on the packages available in R for modeling Bayesian Networks. Here’s the spreadsheet that came out of the meeting:

NINO 3 is a measure of El Niño/La Niña (ENSO) intensity. It’s often said that ENSO has a period of 3-7 years. Why is it so hard to identify a single period? This is a job for SSA!

The scree plot shows that there are a ton of significant eigenvectors– this is a very complicated signal– but I want to focus on the first four, which are a shade more significant than the rest. They come in two pairs, which for SSA means that they represent two sinusoids. Here are the eigenvectors.

How well do these two capture the ENSO signal?

Not great, but it gets most of the ups and downs. Just not the peaks.

So what are the periods of these two? The first two are at 5.8 years, and the second two at 3.5 years. To note: those are relatively prime, so the two frequencies are always going to be going in and out of sync. So it’s not so much that ENSO has a clear frequency (it doesn’t), nor that that frequency is 3-7 years (because what would that mean anyway). It has two main frequencies, like a dial tone.