# Statistics

## Causal inference in ecology – Counterfactuals

Causal inference in ecology – links to the series

• A causal factor (“cause” for short) is something that is predictably related to a particular outcome. For example, fertilizing crops generally increases their yield, so fertilizer is a causal factor related to yield. The way I think about it, a causal factor need not always lead to the outcome. It’s enough if it merely increases the probability of the outcome. For example, smoking doesn’t always lead to lung cancer among those who smoke, but it does increase the probability that you will suffer from lung cancer if you smoke.
• Causes precede effects.2 That’s one reason why teleology is problematic. A teleological explanation explains the current state of things as a result of, i.e., as caused by, something in the future, namely a purpose.3
• Effects may have multiple causes. The world, or at least the world of biology, is a complicated place. Regardless of what phenomenon you’re studying, there are likely to be several (or many) causal factors that influence.

The last point is one of the most important ones for purposes of this series. When we are investigating a phenomenon,4 we’re trying to discern which of several plausible causal factors plays a role and, possibly, the relative “importance” of those causal factors.5

To make this concrete, let’s suppose that we’re trying to determine whether application of nitrogen fertilizer increases the yield of corn. That means we have to determine whether adding nitrogen and adding nitrogen alone increases corn yield. Why the emphasis on “adding nitrogen alone”? Suppose that we added nitrogen to a corn field by adding manure. Then increases in the amount of applied nitrogen are associated with increases in the amount of a host of other substances. If yields increased, we’d know that adding manure increases yield, but not whether it’s because of the nitrogen in manure or something else. Why does this matter?

From very early on in our education we’re taught that “correlation is not the same as causation.” We want to distinguish cases where A causes B from cases where A is merely correlated with B. Yet, as David Hume pointed out long ago, experience6 alone can only show us that A and B actually occur together, not that they must occur together (link). One way of distinguishing cause from correlation is that causes support counterfactual statements. They provide us with a reason to believe statements like “If we had applied nitrogen to the field, the corn yield would have increased” even if we never applied nitrogen to the field at all. The only reason I can see that we could believe such a statement is if we had already determined that adding nitrogen and adding nitrogen alone increases corn yield.7

How do we determine that? Randomized controlled experiments are the most widely known approach, and they are typically regarded as the gold standard against which all other means of inference are compared. That’s where we’ll pick up in the next installment.

1. As I warned in the introduction to the series, I am not an expert in causal inference. The terminology I use is likely both to be imprecise and to be somewhat different from the terminology experts use.
2. Philosophers have argued about whether backward causation is possible, but I’m going to ignore that possibility.
3. Biologists sometimes use teleological language to explain adaptation, e.g., land animals evolved legs to provide mobility. It is, however, relatively easy (if a bit long-winded) to eliminate the teleological language, because natural selection shows how adaptations arise from differential reproduction and survival (link).
4. Or at least this is how it is when I’m investigating a phenomenon.
5. I’ll come back to the idea of identifying the relative importance of causal factors in a future post.
6. Or experiment.
7. If there are any philosophers reading this, you’ll recognize that this account is horribly sketchy and amounts to little more than proof by vigorous assertion. If you’re so inclined, I invite you to flesh out more complete explanations for readers who are interested.

## Causal inference in ecology – Introduction to the series

If you’ve been following posts here since the first of the year, you know that I’ve been writing about how I keep myself organized. Today I’m starting a completely different series in which I begin to collect my thoughts on how we can make judgments about the cause (or causes) of ecological phenomena1 and the circumstances under which judgments are possible. Before I start, I need to offer a few disclaimers.

• Any evolutionary biologist or ecologist who knows me and my work knows that it’s not uncommon for my ideas to represent a minority opinion. (Think pollen discounting for those of you who know my work on the evolution of plant mating systems.) I make no claim that anything I write here is broadly representative of what my fellow evolutionary biologists and ecologists think, only that it’s what I think. Please challenge me on anything you think I’ve got wrong, because I’m sure there will be things I get wrong, and the easiest way for me to discover those errors is for someone else to point them out.
• I had a minor in Philosophy as an undergraduate and there is an enormous literature on causality in the philosophy of science. I’ll be using a very crude understanding of “cause.” I don’t think it is wildly misleading, but I’m certain it wouldn’t stand up to serious scrutiny.2
• I’ll be thinking about causal inference in the specific context of trying to infer causes from observational data using statistics rather than from inferring causes controlled experiments.3 I’ll be using an approach developed in the 1970s by Donald Rubin, the Rubin Causal Model.4
• There is a very large literature on causal inference in the social sciences. I’ll be drawing heavily on Imbens and Rubin, Causal Inference for Statistics, Social and Biomedical Sciences: An Introduction,5 but there’s an enormous amount of material there that I won’t attempt to cover. I am also pretty new to the concepts associated with the Rubin causal model, so it’s entirely possible that I’ll misrepresent or misinterpret a point that the real experts got right. In other words, if something I say doesn’t make any sense, it’s more likely I got it wrong than that Imbens and Rubin got it wrong.

Although I will be thinking about causal inference in the context of observational data and statistics, I don’t plan to write much (if at all) about the problems with P-values, Bayes factors, credible/confidence intervals overlapping 0 (or not), and the like. If you’d like to know the concerns I have about them, here are links to old posts on those issues.

1. I’m calling the post “Causal inference in ecology” only because “Causal inference in ecology, evolutionary biology, and population genetics” would be too long.
2. There’s a good chance that a moderately competent undergraduate Philosophy major would find it woefully inadequate.
3. To be more precise, we don’t infer causes from controlled experiments. Rather, we have pre-existing hypotheses about possible causes, and we use controlled experiments to test those hypotheses.
4. In my relatively limited reading on the subject, I’ve most often seen it referred to as the Rubin causal model, but it is sometimes referred to as the Neyman causal model.
5. Reminder: If you click on that link, it will take you to Amazon.com. I use that link simply because it’s convenient. You can buy the book, if you’re so inclined, from many other outlets. I am not an Amazon affiliate, and I will not receive any compensation if you decide to buy the book regardless of whether you buy it at Amazon or elsewhere. By the way, Chapter 23 in Gelman and Hill’s book, Data Analysis Using Regression and Multilevel/Hierarchical Models has an excellent overview of the Rubin causal model.

## Exploring mixed models in Stan

I am about to begin developing a moderately complex mixed model in Stan to analyze realtinoships among anatomical/morphological traits (e.g., leaf thickness, LMA, wood density), physiological performance (e.g., Amax, stem hydraulic conductance), and indices of fitness (e.g., height, growth rate, number of seedheads). One complication is that the observations are from several different species of Protea at several different sites.1 We’re going to treat sites as nested within species.

Before I start building the whole model, I wanted to make sure that I can do a simple mixed linear regression with a random site effect nested within a random species effect. In stan_lmer() notation that becomes:

stan_lmer(Amax ~ LMA + (1|Species/Site))


I ran a version of my code with several covariates in addition to LMA using hand-coded stan and compared the results to those from stan_lmer(). Estimates for the overall intercept and the regression coefficients associated with each covariate were very similar. The estimates of both standard deviations and individual random effects at the species and site within species level were rather different – especially at the species level. This was troubling, so I set up a simple simulation to see if I could figure out what was going on. The R code, Stan code, and simulation results are available in Github: https://kholsinger.github.io/mixed-models/.

The model used for simulation is very simple:

$$y_k \sim \mbox{N}(\mu_k, \sigma) \\ \mu_k = \beta_0(species|site) + \beta_1x \\ \beta_0(species|site) \sim \mbox{N}(\beta_0(species), \sigma_{species|site}) \\ \beta_0(species) \sim \mbox{N}(\beta_0, \sigma_{species})$$

Happily, the Stan code I wrote does well in recovering the simulation parameters.2 Surprisingly, it does better on recovering the random effect parameters than stan_lmer(). I haven’t completely sorted things out yet, but the difference is likely to be a result of different prior specifications for the random effects. My simulation code3 uses independent Cauchy(0,5) priors for the standard deviation of all variance parameters. stan_lmer() uses a covariance structure for all parameters that vary by group.4 If the difference in prior specifications is really responsible, it means that the differences between my approach and the approach used in stan_lmer() will vanish as the number of groups grows.

Since we’re only interested in the analog of $$\beta_1$$ for the analyses we’ll be doing, the difference in random effect estimates doesn’t bother me, especially since my approach seems to recover them better given the random effect structure we’re working with. This is, however, a good reminder that if you’re working with mixed models and you’re interested in the group-level parameters, you’re going to need a large number of groups, not just a large number of individuals, to get reliable estimates of the group parameters.

If you’ve spent any time using R, you probably know the name Hadley Wickham. He’s chief scientist at RStudio, the author of 4 books on R, and the author of several indispensable R packages, including ggplot2, dplyr, and devtools. I was reminded recently that several years ago, he wrote a very useful paper for the Journal of Statistical Software, “Tidy data” (August 2014, Volume 59, Issue 10, https://www.jstatsoft.org/article/view/v059i10).

If you are familiar with Hadley’s contributions to R, you won’t be surprised that tidy data has a simple, clean – tidy – set of requirements:

1. Each variable forms a column.
2. Each observation forms a row.
3. Each type of observational unit forms a table.

That sounds simple, but it requires that many of us rethink the way we structure our data, no more column headers as values, no more storing of multiple variables in one column, no more storing some variables in rows and others in columns. Fortunately, Hadley is also the author of tidyr. I haven’t used it yet, but given how bad I am at starting with tidy data, I suspect I’ll be using it a lot in the future.

## I’m sorry. P < 0.005 won't save you.

Recently, a group of distinguished psychologists posted a preprint on PsyArxiv arguing for a re-definition of the traditional signifance threshold, lowering it from P < 0.05 to P < 0.005. They are concerned about reproducibility, and they argue that “changing the P value threshold is simple, aligns with the training undertaken by many researchers, and might quickly achieve broad acceptance.” That’s all true, but I’m afraid it won’t necessarily “improve the reproducibility of scientific research in many fields.”

Why? Let’s take a little trip down memory lane.

Almost a year ago I pointed out that we need to “Be wary of results from studies with small sample sizes, even if the effects are statistically significant.” I illustrated why with the following figure produced using R code available at Github: https://github.com/kholsinger/noisy-data

What that figure shows is the distribution of P-values that pass the P < 0.05 significance threshold when the true difference between two populations is mu standard deviations (with the same standard deviation in both populations) and with equal sample sizes of n. The results are from 1000 random replications. As you can see when the sample size is small, there’s a good chance that a significant result will have the wrong sign, i.e., the observed difference will be negative rather than positive, even if the between-population diffference is 0.2 standard deviations. When the between-population difference is 0.05, you’re almost as likely to say the difference is in the wrong direction as to get it right.

Does changing the threshold to P < 0.005 help. Well, I changed the number of replications to 10,000, reduced the threshold to P < 0.005, and here’s what I got.

Do you see a difference? I don’t. I haven’t run a formal statistical test to compare the distributions, but I’m pretty sure they’d be indistinguishable.

In short, reducing the significance threshold to P < 0.005 will result in fewer investigators reporting statistically significant results. But unless the studies they do also have reasonably large sample sizes relative to the expected magnitude of any effect and the amount of variability within classes, they won’t be any more likely to know the direction of the effect than with the traditional threshold of P < 0.05.

The solution to the reproducibility crisis is better data, not better statistics.

## Not every credible interval is credible

Lauren Kennedy and co-authors (citation below) worry about the effect of “contamination” on estimates of credible intervals.1 The effect arises because we often assume that values are drawn from a normal distribution, even though there are “outliers” in the data, i.e., observations drawn from a different distribution that “contaminate” our observations. Not surprisingly, they find that a model including contamination does a “better job” of estimating the mean and credible intervals than one that assumes a simple normal distribution.2

They consider the following data as an example:
-2, -1, 0, 1, 2, 15
They used the following model for the data (writing in JAGS notation):

x[i] ~ dnorm(mu, tau)
tau ~ dgamma(0.0001, 0.0001)
mu ~ dnorm(0, 100)


That prior on tau should be a red flag. Gelman (citation below) pointed out a long time ago that such a prior is a long way from being vague or non-informative. It puts a tremendous amount of weight on very small values of tau, meaning a very high weight on large values of the variance. Similarly, the N(0, 100); prior on mu; may seem like a “vague” choice, but it puts more than 80% of the prior probability on outcomes with x < -20 or x > 20, substantially more extreme than any that were observed.

Before we begin an analysis we typically have some idea what “reasonable” values are for the variable we’re measuring. For example, if we are measuring the height of adult men, we would be very surprised to find anyone in our sample with a height greater than 3m or less than 0.5m. It wouldn’t make sense to use a prior for the mean that put appreciable probability on outcomes more extreme.

In this case the data are made up, so there isn’t any prior knowledge to work from. but the authors say that “[i]t is immediately obvious that the sixth data point is an outlier” (emphasis in the original). Let’s take them at their word. A reasonable choice of prior might then be N(0,1), since all of the values (except for the “outlier”) lie within two standard deviations of the mean.3 Similarly, a reasonable choice for the prior on sigma (sqrt(1/tau)) might be a half-normal with mean 0 and standard deviation 2, which will allow for standard deviations both smaller and larger than observed in the data.

I put that all together in a little R/Stan program (test.R, test.stan). When I run it, these are the results I get:

         mean se_mean    sd    2.5%     25%     50%     75%   97.5% n_eff  Rhat
mu      0.555   0.016 0.899  -1.250  -0.037   0.558   1.156   2.297  3281 0.999
sigma   4.775   0.014 0.841   3.410   4.156   4.715   5.279   6.618  3466 1.000
lp__  -16.609   0.021 0.970 -19.229 -17.013 -16.314 -15.903 -15.663  2086 1.001


Let’s compare those results to what Kennedy and colleagues report:

AnalysisPosterior mean95% credible interval
Stan + "reasonable priors"0.56(-1.25, 2.30)
Kennedy et al. - Normal2.49(-4.25, 9.08)
Kennedy et al. - Contaminated normal0.47(-2.49, 4.88)

So if you use “reasonable” priors, you get a posterior mean from a model without contamination that isn’t very different from what you get from the more complicated contaminated normal model, and the credible intervals are actually narrower. If you really think a priori that 15 is an unreasonable observation, which estimate (point estimate and credible interval) would you prefer? I’d go for the model assuming a normal distribution with reasonable priors.

It all comes down to this. Your choice of priors matters. There is no such thing as an uninformative prior. If you think you are playing it safe by using very vague or flat priors, think carefully about what you’re doing. There’s a good chance that you’re actually putting a lot of prior weight on values that are unreasonable.4 You will almost always have some idea about what observations are reasonable or possible. Use that information to set weakly informative priors. See the discussion at https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations for more detailed advice.

## Against null hypothesis testing – the elephants and Andrew Gelman edition

Last week I pointed out a new paper by Denes Szucs and John Ioannidis, When null hypothesis significance testing is unsuitable for research: a reassessment.1 I mentioned that P-values from small, noisy studies are likely to be misleading. Last April, Raghu Parthasarathy at The Eighteenth Elephant had a long post on a more fundamental problem with P-values: they encourage binary thinking. Why is this a problem?

1. “Binary statements can’t be sensibly combined” when measurements have noise.
2. “It is almost never necessary to combine boolean statements.”
3. “Everything always has an effect.”

Those brief statements probably won’t make any sense,2 so head over to The Eighteenth Elephant to get the full explanation. The post is a bit long, but it’s easy to read, and well worth your time.

Andrew Gelman recently linked to Parthasarathy’s post and adds one more observation on how P-values are problematic: they are “interpretable only under the null hypothesis, yet the usual purpose of the p-value in practice is to reject the null.” In other words, P-values are derived assuming the null hypothesis is true. They tell us what the chances of getting the data we got are if the null hypothesis were true. Since we typically don’t believe the null hypothesis is true, the P-value doesn’t correspond to anything meaningful.

To take Gelman’s example, suppose we had an experiment with a control, treatment A, and treatment B. Our data suggest that treatment A is not different from control (P=0.13) but that treatment B is different from the control (P=0.003). That’s pretty clear evidence that treatment A and treatment B are different, right? Wrong.

P=0.13 corresponds to a treatment-control difference of 1.5 standard deviations; P=0.003, to a treatment-control difference of 3.0 standard deviations, a difference of 1.5 standard deviations, which corresponds to a P-value of 0.13. Why the apparent contradiction? Because if we want to say that treatment A and treatment B, we need to compare them directly to each other. When we do so, we realize that we don’t have any evidence that the treatments are different from one another.

As Parthasarthy points out in a similar example, a better interpretation is that we have evidence for the ordering (control < treatment A < treatment B). Null hypothesis significance testing could easily mislead us into thinking that what we have instead is (control = treatment A < treatment B). The problem arises, at least in part, because no matter how often we remind ourselves that it’s wrong to do so, we act as if a failure to reject the null hypothesis is evidence for the null hypothesis. Parthasarthy describes nicely how we should be approaching these problems:

It’s absurd to think that anything exists in isolation, or that any treatment really has “zero” effect, certainly not in the messy world of living things. Our task, always, is to quantify the size of an effect, or the value of a parameter, whether this is the resistivity of a metal or the toxicity of a drug.

We should be focusing on estimating the magnitude of effects and the uncertainty associated with those estimates, not testing null hypotheses.

## Against null hypothesis significance testing

Several months ago I pointed out that P-values from small, noisy experiments are likely to be misleading. Given our training, we think that if a result is significant with a small sample, it must be a really big effect. But unless we have good reason to believe that there is very little noise in the results (a reason other than the small amount of variation observed in our sample), we could easily be misled. Not only will we overestimate how big the effect is, but we are almost as likely to say that the effect is positive when it’s really negative as we are to get the sign right. Look back at this post from August to see for yourself (and download the R code if you want to explore further). As Gelman and Carlin point out,

There is a common misconception that if you happen to obtain statistical significance with low power, then you have achieved a particularly impressive feat, obtaining scientific success under difficult conditions.

I bring this all up again because I recently learned of a new paper by Denes Szucs and John Ioannidis, When null hypothesis significance testing is unsuitable for research: a reassessment. They summarize their advice on null hypothesis significance testing (NHST) in the abstract:

Whenever researchers use NHST they should justify its use, and publish pre-study power calculations and effect sizes, including negative findings. Studies should optimally be pre-registered and raw data published.

They go on to point out that part of the problem is the way that scientists are trained:

[M]ost scientists…are still near exclusively educated in NHST, they tend to misunderstand and abuse NHST and the method is near fully dominant in scientific papers.

## Designing exploratory studies: measurement (part 2)

Remember Gelman’s preliminary principles for designing exploratory studies:

1. Validity and reliability of measurements.
2. Measuring lots of different things.
3. Connections between quantitative and qualitative data.
4. Collect or construct continuous measurements where possible.

I already wrote about validity and reliability. I admitted to not knowing enough yet to provide advice on assessing the reliability of measurements ecologists and evolutionists make (except in the very limited sense of whether or not repeated measurements of the same characteristic give similar results). For the time being that means I’ll focus on

• Remembering that I’m measuring an indicator of something else that is the thing that really matters, not the thing that really matters itself.
• Being as sure as I can that what I’m measuring is a valid and reliable indicator of that thing, even though the best I can do with that right now is a sort of heuristic connection between a vague notion of what I really think matters, underlying theory, and expectations derived from earlier work.

It’s that second part where “measuring lots of different things” comes in. Let’s go back to LMA and MAP. I’m interested in LMA because it’s an important component of the leaf economics spectrum. There are reasons to expect that tough leaves (those in which LMA is high) will not only be more resistant to herbivory from generalist herbivores, but that they will have lower rates of photosynthesis. Plants are investing more in those leaves. So high LMA is, in some vague sense, an indicator of the extent to which resource conservation is more important to plants than rapid acquisition of resources. So in designing an exploratory study, I should think about other traits plants have that could be indicators of resource conservation vs. rapid resource acquisition and measure as many of them as I can. A few that occur to me are leaf area, photosynthetic rate, leaf nitrogen content, leaf C/N ratio, tissue density, leaf longevity, and leaf thickness.

If I measure all of these (or at least several of them) and think of them as indicators of variation on the underlying “thing I really care about”, I can then imagine treating that underlying “thing I really care about” as a latent variable. One way, but almost certainly not the only way, I could assess the relationship between that latent variable and MAP would be to perform a factor analysis on the trait dataset, identify a single latent factor, and use that factor as the dependent variable whose variation I study in relation to MAP. Of course, MAP is only one way in which we might assess water availability in the environment. Others that might be especially relevant for perennials with long-lived leaves (like Protea) in the Cape Floristic Region rainfall seasonality, maximum number of days between days with “significant” rainfall in the summer, total summer rainfall, estimated potential evapotranspiration for the year, and estimated PET for the summer. A standard way to relate the “resource conservation” factor to the “water availability” factor would be a canonical correspondence analysis.

I am not advocating that we all start doing canonical correspondence analyses as our method of choice in designing exploratory studies, this way of thinking about exploratory studies does help me clarify (a bit) what it is that I’m really looking for. I still have work to do on getting it right, but it feels as if I’m heading towards something analogous to exploratory factor analysis (to identify factors that are valid, in the sense that they are interpretable and related in a meaningful way to existing theoretical constructs and understanding) and confirmatory factor analysis (to confirm that the exploration has revealed factors that can be reliably measured).

Stay tuned. It is likely to be a while before I have more thoughts to share, but as they develop, they’ll appear here, and if you follow along, you’ll be the first to hear about them.

## Designing exploratory studies: measurement

On Wednesday I argued that we need carefully done exploratory studies to discover phenomena as much as we need carefully done experimental studies to test explanations for the phenomena that have been discovered.1 Andrew Gelman suggests four preliminary principles:

1. Validity and reliability of measurements.
2. Measuring lots of different things.
3. Connections between quantitative and qualitative data.2
4. Collect or construct continuous measurements where possible.

Today I’m going to focus on #1, the validity and reliability of measurements.

If there happen to be any social scientists reading this, it’s likely to come as a shock to you to learn that most ecologists and evolutionary biologists haven’t thought too carefully about the problem of measurement, or at least that’s been my experience. My ecologist and evolutionary biologist friends are probably scratching their heads. “What the heck does Holsinger mean by ‘the problem of measurement.'” I’m sure I’m going to butcher this, because what little I think I know I picked up informally second hand, but here’s how I understand it.