Causal Model Structures and Random Effects

Causal Inference from Different Model Structures

The world is a complex place. Often, when we have hierarchical data - particularly data where we have repeated observations of the same locations over time - deducing causal inference with confidence is a tricky business. While a default for estimation might be to throw in all of the predictors you’re interested in, and then slosh everything else into a site-level random effect, that doesn’t guaruntee proper causal inference. I’m going to discuss some reasons why drawing heavily on some principles as outlines by the recent Bell et al. 2017 paper which discussed the core nugget of the literature as well as some great review by notes from Marc Bellmare’s course notes on causal inference from observational data. Along the way I’ll talk a bit about shutting the back door (sensu Pearl) and how random effects relate to causal model structures and the back door problem.

Let’s start with a biological scenario. Let’s say that we have a system we’re looking at where we thing that snail abundance might be influenced by temperature. We sample the same sites for 15 years, and want to build a model where we look at whether temperature influences snail abundance as sampled by a site.

There are just two catches.

First, as I simulate this data, yes, snail abundance is influenced by temperature. However, the sites each have a different average temperature. This is something many of us who sample along biogeographic gradients are familiar with. Second, recruitment of new snails correlates with site average temperature. So, perhaps, sites that experience warmer temperatures on average also have higher recruitment due to site characteristics. There’s also some residual uncorrelated site effects unrelated to temperature or recruitment.

For fellow nerds keeping track, I made that and other figures with http://rich-iannone.github.io/DiagrammeR/

Or, here’s the simulation of the data itself. This is one way to do it, there are others.

set.seed(nchar("A causal model  odyssey"))

#the basic physical data
my_df <- crossing(tibble(site = letters[1:10], 
                  site_avg_temp = rnorm(10),
                  recruitment = rnorm(10, site_avg_temp)),
                  deviation_from_site_avg = rnorm(15)) %>%
  mutate(temp = 3*site_avg_temp + deviation_from_site_avg) %>%
  arrange(ord = rnorm(150)) %>% #randomness!
  group_by(site) %>%
  mutate(sample_time = 1:15) %>%
  ungroup()

Let’s see how temperature varies over time.

qplot(sample_time, temp, data = my_df, color = site, geom="line")

OK, now we’ll add recruitment which correlates with site-level temperature.

#A site-level covariate that correlated with site average temperature
#my_df <- my_df %>%
#  mutate(recruitment = rnorm(n(), site_avg_temp))

What does that correlation look like?

qplot(recruitment, site_avg_temp, data = my_df, color = site, geom="point")

Finally, snails are a sum of temperature and recruitment.

#OK, snail abund is a function of temperature and recruitment
my_df <- my_df %>%
  mutate(snails = rnorm(n(), temp + recruitment))

The Problem

This is all well and good, and if we had measurements of everything above, we could do wonders. But we don’t. All we have is temperature at each site i at each point in time j and snail abundance. And, in fact, while I’ve outlined recruitment as a site-level variabilty that would affect temperature, who knows what elses is there that correlated with temperature! The best we can do is something that has this sort of structure

This despite the underlying relationships in our nature (where we know all of the data) looking like

ggpairs(my_df %>% select(snails, temp, recruitment, site_avg_temp))

OK, so, how do we craft our model so that we can get a causally valid estimate of temperature on snails? One that accomodates the fact that temperature at the site level could correlate with other site-level covariates and makes sure that site-level covariates are not contaminating our estimate of the temperature effect due to non-independence between error estimates in the model and our estimates of the temperature effect?

There are a few options. Let’s start with what is a common first approach.

(see also https://t.co/0uINKPb7EA)

The Models

my_df <- my_df %>%
  group_by(site) %>%
  arrange(sample_time) %>%
  mutate(delta_temp = temp - lag(temp),
         delta_snails = snails - lag(snails)) %>%
  ungroup()

Comparisons Between Approaches

library(nlme)

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:dplyr':
## 
##     collapse

#Now, on to models
ran_mod <- lme(snails ~ temp,
               random = ~1|site,
               data= my_df)

fix_mod <- lm(snails ~ temp + site, data = my_df)
 
mundlak_mod <- lme(snails ~ temp + site_avg_temp,
                    random = ~1|site,
                    data= my_df)
 
gmc_mod <- lme(snails ~ deviation_from_site_avg + site_avg_temp,
                random = ~1|site,
                data= my_df) 

diff_mod <- lme(delta_snails ~ delta_temp,
                random = ~1|site,
                data = my_df,
                na.action = na.omit)

ran_mod %>% broom::tidy(effects="fixed") %>% kable("html",3) %>% kable_styling

(fix_mod %>% broom::tidy())[1:2,]%>% kable("html",3) %>% kable_styling

mundlak_mod %>% broom::tidy(effects="fixed")%>% kable("html",3) %>% kable_styling

gmc_mod %>% broom::tidy(effects="fixed")%>% kable("html",3) %>% kable_styling

diff_mod %>% broom::tidy(effects="fixed")%>% kable("html",3) %>% kable_styling

my_df <- my_df %>%
  mutate(snails_2 = rnorm(n(), temp - recruitment))

ran_mod_2 %>% broom::tidy(effects="fixed") %>% kable("html",3) %>% kable_styling

(fix_mod_2 %>% broom::tidy())[1:2,]%>% kable("html",3) %>% kable_styling

mundlak_mod_2 %>% broom::tidy(effects="fixed")%>% kable("html",3) %>% kable_styling

gmc_mod_2 %>% broom::tidy(effects="fixed")%>% kable("html",3) %>% kable_styling

diff_mod_2 %>% broom::tidy(effects="fixed")%>% kable("html",3) %>% kable_styling

Restricted by Our Assumptions

note, what follows is a lightly edited commentary from Laura Dee, my guru/fellow traveller in these discussions.

Let’s compare and contrast the pros and cons of the different approaches. For instance, of all of the approaches here, the FE and first-differences make the weakest assumptions here to get rid of confounders at the site-level and to feel confident in removing bias for a causal temperature effect. But! they are less efficient at estimating the error (which, in small ecological datasets, could be an issue - although it’s not ginormous here) AND THE BIGGIE - we cannot study the between-site effects. They are not causally identified (and you can see how the recruitment and site-level temperature effect get swept together). The only resolution is more measurements of more covariates to clean the signal and shut the causal back door.

In contrast, RE are very efficient and allow us to study between effects, but they make two very strong assumptions for causal inference, that are likely not met in obs data most of the time. First, that there is no correlation between the RE and the causally estimated quantity (so, between site and temperature, in this case). If we are only modeling the error for prediction/forecasting, maybe this isn’t a big problem. But it’s potentially a large problem if we care about bias in the coefficient estimates - and being able to say that there is a causal connection here. The other approaches assume that the correlation between site and the causal variable of interest is soaked up in the hierarchical predictor. From the examples above, you can see how they improve performance, and have real causal meaning relative to the RE approach. However, if there is an additional correlated hierarchical predictor not accounted for that correlates with the causal driver of interest, they again lose causal identification. Hence a straight FE or Differencing model might be better. In contrast, including a hierarchical predictor can yield additional information and pull out a particular signal of interest.
The 2nd assumption of RE is that all units (here sites) follow the same distribution (same mean, variance) AND we have correctly specified what that distribution is. That is why RE’s are more efficient (we are adding an assumption to get this though – specifying the distribution and assuming one/that it is correct). Any site-level other predictors are randomly distributed between sites, thus ensuring exchangability. This is often fairly reasonable, but, when we get into things like spatial gradients, can be difficult if not impossible to justify, and thus require a different model (e.g., incorporating that site-level predictor, or a FE approach). Given the proliferation of space-for-time designs cropping up in Ecology, it’s worth a long think about.

In contrast, fixed effects and first differencing allows each unit to have their own mean – but does not require that they necessarily follow the same distribution. Thus, exchangability is no longer a problem, even if we have not included hierarchical predictors.

But what about the SE? Well, there are ways to deal with that not commonly used in Ecology currently - cluster robust standard errors. Let’s compare with clustered robust standard errors (different than the FE’s and you have to do this separately). Clustered robust standard errors are very common in Econ – they allow for arbitrary correlation within the unit and through time (i.e., within the clusters, allows anything in the block for the off diagonal in the covariance matrix). They do not require specifying a distribution – and therefore make weaker assumptions but are also less efficient. But, with most data, it’s not a huge efficiency loss.

As an example, we (Laura et al. - Jarrett had no part in this!) did this in this 2016 Proc B paper, and provide code for clustering the standard errors on large marine ecosystem (the unit of the cluster): http://rspb.royalsocietypublishing.org/content/royprsb/suppl/2016/08/17/rspb.2016.1435.DC1/rspb20161435supp3.pdf (written by Steve Miller - all hail!). There is also the sandwich package in R (Jarrett - mmmmmmm….sandwich….) – sandwich errors are various forms of robust standard errors, the title is “Robust Covariance Matrix Estimators” and the package has “Model-robust standard error estimators for cross-sectional, time series, clustered, panel, and longitudinal data.”

The use of robust standard errors something that is not common in ecology – one may well ask, why not?? BTW: We could probably use a good Sandwiches for Ecology paper in Methods in Ecology and Evolution.

So, there are two things going on here – what each approach does to eliminate bias in coefficient estimates (i.e., assumptions about how the different approaches deal with confounders) and the need to account for correlations in the error (very easy with few assumptions with clustered robust standard errors, and generally not a concern or focus of most analyses in econ because there are simple and reliable solutions).

Are we 100% Certain About that Causal Link with Temperature

So, is there a causal link between temperature and snails, given these better models? No! Why not? Well, while we’ve accounted for site-level covariates, we have not accounted for inter-annual covariates that covary with deviation from site-level temperature. Like, say, nutrients! Or algal cover! Or other things. Yes, we can jigger our model structure a bit more to account for this, but, given that we have one sample point per site per year, we quickly run out of DF and inflate our SE of any estimate. More samples per site per year can help, but, it rapidly becomes a balancing act of model complexity and reasonable inference. Moreover, we really have to focus on what is our question, and what is the right model to nail down that causal inference. Someties, we just do not possess the covariates necessary to answer a question! And that is a design problem, not something that can be solved by models alone.

Conclusion

Rather than leave things on a dark-ish note, I hope that this musing has showed you that to derive causally valid answers from observational data sampled over time, we do need to be a hair more cautious that perhaps we’ve been in the past. There are, however, some elegant and beautiful solutions with respect to model structure. Me, I’ve become a big fan of group mean centered models and difference models, as their meaning is always clear, and I can have a high degree of confidence in their estimates.