What is Modern Statistical Workflow?

Jim Savage

Or: It’s not the methods or languages you use; it’s how you do them

Building good statistical models is hard. Unfortunately, most statistics/data science/econometrics training focuses on statistical models and the algorithms that we use to estimate them, at the cost of not covering other important topics. One such important topic is workflow: how to structure the process of your analysis to maximise the odds that you build useful models.

The following is the workflow that I try to force on myself. Often I take shortcuts, only to get stuck and come back to the workflow. When friends email me because their models aren’t doing what they’re meant to, I tell them to stick to this workflow. Put simply, this workflow is the best way to learn a wide range of modeling techniques and build better models.

A more elaborate explanation I wrote earlier in the year is here.

1. Plot your data

Plotting should always be your starting point, even if you are a whiz modeler. A few good reasons: first, you are either trying to discover relationships in the data or explain them away; your eyes will be able to detect these very easily. Your data might have problems, like missing values or incorrectly-coded observations. Sometimes your variables have different orders of magnitude. Plotting your data first is the easiest way of avoiding these pitfalls.

You should always produce density plots of your outcome variables. Doing this will provide guidance for the sort of model you should build.

2. Write down what you know both in and out of sample, and what you only know in sample.

As I wrote in my rules of thumb, you need to have a very clear idea of what the random variables are in your model. Random variables are simply those variables whose values are in part due to chance. For most modeling purposes, our random variables are the things we don’t know for sure out of sample. These might include model parameters, latent variables, predictions, etc. [Edit: I don't mean literally that model parameters are random; it's our understanding of them which has uncertainty].

A common problem occurs when the modeler uses a random variable as a predictive feature in a model but does not explicitly model it. For instance they might build a model that looks a bit like

Salest=f(Weathert,Day of the weekt)+errort$$$$

Next, when called on to make a prediction, the modeler uses forecasts for weather to generate predictions for sales, but without taking into account the uncertainty around weather (a random variable). As for the day of the week, this is not a random variable and so we don’t have to model it. Any forecasts conditioned on random variables without taking into account their uncertainty will be far too precise.

The reason we should write down what the random variables are in the model is because this is precisely what we are going to model.

3. Build a generative model of those things that you don’t know out of sample.

In this step, we ask ourselves: what is a plausible process that could generate the outcomes that we observe? For instance, if we think that a normal linear regression model with coefficients β and covariates X and residual standard deviation σ is suitable, our generative model would be

yi∼N(Xβ,σ)$$$$

Or we might consider a normality assumption to be too strong, and use a “fat-tailed distribution” instead

yi∼Student's t(ν,Xβ,σ)$$$$

Or perhaps our outcome yi comes from two distribtuions, each with a different probability (as in this post). Or it could be binary, or count data, or strictly positive data, or multimodal data, etc., in which we would choose different distributions still.

Note that the examples above are extremely simple models—you should almost always start with simple models and build up in complexity. As your model grows in complexity, the value to performing the fake data exercise in steps 4 and 5 grows.

After defining the generative model, you should assign some priors to all the unknowns—in this case, the parameters ν, β, and σ. These priors should give weight to plausible values of the parameters of the model, and no weight to impossible values. For instance, ν is restricted to be >1 and σ has to be positive. Priors for those parameters should not put weight on values outside that range.

4. Draw some data from the generative model with some known parameters drawn from the prior.

We have a generative model for our data—a way to simulate plausible values for y given X—and priors for the parameters.

The next step is to draw some values from a prior, which we treat as being “known” values of the parameters. After doing this we have values for X, and “known” values for ν, β and θ, so we can simulate some fake data by drawing observations from the generative model in step 3.

Why should we simulate some fake data? First, it gives us an idea of whether our model puts weight on impossible outcomes—we don’t want to use a model that does that! But more importantly, this (often skipped) step makes us be very explicit about all the assumptions in the model, and guides us to the estimation in the next step.

5. Estimate the model on the fake data. Can you recover the parameters? Is the model identified?

Before taking your estimation model to real data, you should always try estimating the model on the fake data you simulated in step 4. Why do this? First, you know the values of the parameters for that data. So you should check to see that when you estimate your model on fake data you can recapture the known values. If your model is unable to recapture known parameters with fake data, it will definitely estimate the wrong parameter values using real data.

6. Estimate the model on real data

If your model is able to recapture known parameter values, it’s time to estimate the model on real data.

Often people jump to this step without performing 4 and 5 first, and get funny results. (I know this because some of these people are my friends and I receive a few emails a week about precisely this problem). Gelman’s folk theorem is “it’s probably not the computer. It’s probably your model.” and this is almost always the case here.

7. Check for model convergence

Sometimes, especially for big, loosely-identified models, you might not be getting very good samples from your posterior. One great tool in R for exploring pathologies in sampling is shinystan (available here), which provides a web interface to your MCMC fits.

If you have poor convergence or pathologies in sampling, this can often be fixed by reparameterizing your model. Reparameterizing a model is simply a way of expressing the same model in a way that your posterior has a more regular shape (and so is easier to sample from).

8. Posterior predictive checking and inference

Now we know that the model has been built well and is estimating fine. But was it the right model in the first place? Posterior predictive checking is a very useful method for answering this question. The aim is to check to see if, when simulated, the model generates predictions that have a similar distribution to the observed data, after taking into account uncertainty in the model parameters.

An example of this is below, from some recent work of mine modeling micro-loan repayments in Sub-Saharan Africa.

9 Iterate!

We have just built a fairly simple model, but now we have it working, and we probably have a good idea what’s wrong with it. At this point we can afford to go back to step 3 and build up a more complex model, knowing that if anything breaks (or a deadline approaches), we have a well-built, well-checked model ready to go.

My zen of modeling

1. Your model should have some theoretical basis.
2. Your model, when simulated, should produce outcomes with a similar density to the observed values. Similarly, your model should not place weight on the impossible (like negative quantities, or binary outcomes that aren’t binary). It should place non-zero weight on possible but unlikely outcomes.
3. Think deeply about what is a random variable and what is not. A good rule of thumb: random variables are those things we do not know for certain out of sample. Your model is a joint density over the random variables.
4. You never have enough observations to distinguish one possible data generating process from another process that has different implications. You should model both, giving both models weight in decision-making.
5. The point of estimating a model on a big dataset is to estimate a rich model (one with many parameters). Using millions of observations to estimate a model with dozens of parameters is a waste of electricity.
6. Unless you have run a very large, very well-designed experiment, your problem has unobserved confounding information. If this problem does not occupy a lot of your time, you are doing something wrong.
7. Fixed effects normally aren’t. Mean reversion applies to most things, including unobserved information. Don’t be afraid to shrink.
8. Relationships observed in one group can almost always help us form better understanding of relationships in another group. Learn and use partial pooling techniques to benefit from this.
9. For decision-making, your estimated standard deviations are too small; your estimated degrees of freedom are too big, or your have confused one for the other. Remember, the uncertainty produced by your model is the amount of uncertainty you should have if your model is correct and the process you are modeling does not change.
10. You always have more information than exist in your data. Be a Bayesian, and use this outside information in your priors.

Why you should be posterior predictive checking

Jim Savage

26 August 2016

This week at work I had an excellent problem come across my desk. A client, for whom we had done some modeling work, sent across a plot that looked a little like the one below (this is all fake data).

The client wanted to know if they should be concerned that the distribution of actual outcomes (the top panel) look nothing like the predictions (bottom panel).

If you’ve never played with models before and you saw this chart, your natural response might be to think that the model was terrible. And a terrible model might produce predictions that look nothing like the actual outcomes. Yet in this case, you would be mistaken in writing off the model due to the difference in plots. The predictions in the chart above actually come from the exact same model that I used to generate the “Actual outcome” data—it is impossible to get a better set of predictions than those summarized in the bottom panel. So what’s going on? Let’s step through it.

Often when modeling proportions, we get weird-looking distributions like the top panel in the chart above. Often there are spikes at 0 and 1, and some distribution of continuous outcomes between. In this particular example, we have a spike at 1 and a set of proportons roughly between 0.3 and 0.9, skewed towards higher values in this range.

A popular way of modeling these distributions is to break up the problem using a generative model. A generative model asks “what is a plausible mathematical structure that could give rise to the data that we observe.” One way of breaking up the problem would be to consider there actually being two processes.

1. The first process determines whether an observation is 1 or less than 1.
2. If the proportion is not 1, the second process would assign relative probabilities to outcomes in the range (0,1).

Let’s get specific. If β1, β2${\beta }_2$, and β3${\beta }_3$ are coefficient vectors, μ1, μ2, and μ3 are intercept terms and X is a design matrix containing our predictive features, then we might say that the data generating process for an outcome Yi$Y_i$ could be

Step 1:

I(Yi=1)∼Bernoulli(Logit−1(μ1+Xiβ1))$$$$

where I() is the indicator function, and Step 2, in the case that Yi<1:

Yi∼Beta(μ2+Xβ2,μ3+Xβ3)$$$$

For those who’ve not used them before, the Beta(α,β) density is a highly flexible distribution that takes two positive parameters and assigns weight to observations between 0 and 1, with no probability of outcomes outside that interval. The wikipedia entry is helptul.

A nice thing about using a model like this to model proportions is that it assigns some weight to 1 for all values of μ1+Xiβ1, and weight between 0 and 1 otherwise. If you had a spike at 0 also, you could simply replace the first model with a multinomial logit model, and it’d have the same benefits.

Simulating a draw from the model

When we build generative models, we consider the world as observed to be just a single draw from a universe of different potential draws. Our generative model is just a probability density over those potential draws. So we can generate one draw from the model as so (in R).

# Load the relevant libraries
library(reshape2); library(arm); library(ggplot2); library(dplyr)

# Set some known parameter values

# We'll just assume a single covariate x

# parameters for the logit
beta11 <- -1.5
beta12 <- 0.5

# parameters for the beta
beta21 <- 5
beta22 <- 0.5
beta31 <- 2
beta32 <- -0.3

# Number of observations
N <- 1000

# Generate our predictor x
x <- rnorm(N)

# Generate some draws from the beta distribution for the continuous outcome
continuous <- rbeta(N, beta21 + beta22*x, beta31 + beta32*x)

# Generate the probabilities for the binary outcome
prob_binary <- invlogit(beta11 + beta12*x)

# Generate the draws for the binary outcome
binary <- rbinom(N, size = 1, prob = prob_binary)

# The outcome as defined
outcome <- ifelse(binary==1, binary, continuous)

Now we can plot the fake outcome data that we’ve created.

data_frame(outcome) %>% 
  ggplot(aes(x = outcome)) +
  geom_histogram(alpha = 0.3) +
  ggthemes::theme_economist() +
  ggtitle("One draw from our generative model")

So now let’s generate predictions for each value of x for the value of y. How might we do this? Two methods could work.

• The first method works when we have analytically tractable models. We can simply calculate the expected value of y for each value of x using basic algebra.
• The second method is what we do when we have larger, more complex models for which analytical expected values might be hard to derive. In this approach, we’d simulate the model hundreds or thousands of times for each value of x, and take the average of the simulated outcomes.

In this toy example, we will do the first. It happens that the epected value from a Beta distribution with parameters α,β>0 is αα+β. The expected value from a Bernoulli logit model is simply the probability, here Logit−1(μ1+Xβ1). So for a given value of x, our prediction is just the weighted average of the two models:

E[yi|xi]=(1−Logit−1(μ1+xiβ1))×μ2+Xβ2μ2+Xβ2+μ3+Xβ3+Logit−1(μ1+xiβ1)
$$$$

Now if we plot this density, we get something that looks nothing like our data density. Why? The expected value averages across two outcomes in proportion to their probability, so only those observations with an extremely high probability of having a proportion equal to 1 will have an expected value anywhere near 1.

expected_value <- (1 - prob_binary)*((beta21 + beta22*x)/(beta21 + beta22*x+ beta31 + beta32*x)) + prob_binary

data_frame(expected_value) %>% 
  ggplot(aes(x = expected_value)) +
  geom_histogram(alpha = 0.3) +
  ggthemes::theme_economist() +
  ggtitle("Predicted values given x") +
  xlim(0, 1)

So if I can’t validate my model by comparing actual and predicted distributions, what can I do?

In three words, posterior predictive checking. This method involves simulating outcomes from your model, incorporating any uncertainty you have about the parameters of the model after estimating it. The steps

1. Estimate your model. If you are a Bayesian (you should be), and performed MCMC, you will have hundreds or thousands of draws from the joint density of p(μ,β), in proportion to their posterior probability. If you are a frequentist, you should have performed some magic to get many draws from the parameter space (I’m not sure this is actually possible). Common techniques would be bootstrapping or making distributional assumptions about the parameters.
2. For each parameter draw from your model, for each value of X, draw an outcome from the model. These draws make up your posterior predictions.

A worked exercise

Let’s use the realized x and y to estimate this model in Stan and generate posterior predictive draws.

# A Stan program to estimate the model described above

beta_logit_model <- "
data {
  int N; // number of observations
  int P; // number of explanatory variables
  int N2; // number of observations to predict
  vector[N] Y; // the proportion paid of expected at some period ahead
  matrix[N, P] X; // The explanatory variables (make sure no variables known in the future!)
  matrix[N2, P] X_new; // explanatory variables for out-of-sample
}
transformed data {
// we an integer Y and a real Y (Y_adj is only here for flexibility of extension)
  vector[N] Y_adj;
  int Y_int[N];
  
  for(i in 1:N) {
    if(Y[i] == 1.0) {
      Y_adj[i] = 1.0;
      Y_int[i] = 1;
    } else {
      Y_adj[i] = Y[i];
      Y_int[i] = 0;
    }
  }
}
parameters {
  matrix[3, P] beta;
  vector[3] mu;
}
model {
  // priors 
  to_vector(beta) ~ normal(0, .2);
  mu ~ normal(0, 1);

  // likelihood
  Y_int ~ bernoulli_logit(mu[1] + X*beta[1]');
  
  for(i in 1:N) {
    if(Y_int[i]==0) {
      // We constrain alpha and beta to be positive by taking exp()

      Y_adj[i] ~ beta(exp(mu[2] + X[i]*beta[2]'),exp(mu[3] + X[i]*beta[3]'));
    }
  }
}
generated quantities {
  vector[N2] Y_predict;
  int tmp;
  for(i in 1:N2) {
    tmp = bernoulli_rng(inv_logit(mu[1] + X_new[i]*beta[1]'));
    if(tmp ==1) {
      Y_predict[i] = 1.0;
    } else {
      Y_predict[i] = beta_rng(exp(mu[2] + X_new[i]*beta[2]'), exp(mu[3] + X_new[i]*beta[3]'));
    }
    if(is_nan(Y_predict[i])) {
      Y_predict[i] = 1.1;
    }
  }
}

"

Now we estimate the model

library(rstan)
options(mc.cores = parallel::detectCores())

beta_logit_estimation <- stan(model_code = beta_logit_model, 
                              data = list(N = N,
                                          P = 1, 
                                          N2 = N,
                                          Y = outcome,
                                          X = matrix(x, N, 1),
                                          X_new = matrix(x, N, 1)),
                              iter = 600)

And finally plot the posterior replicatons

predictions <- as.data.frame(beta_logit_estimation) %>% 
  select(contains("Y_predict")) %>% 
  melt()

predictions %>% 
  ggplot() +
  geom_line(aes(x = value, group = variable), colour = "orange", stat = "density", alpha = 0.1, adjust = .8) +
  geom_density(data = data_frame(outcome), aes(x = outcome), colour = "black", adjust = 0.8) +
  ggthemes::theme_economist() +
  ggtitle("Actual outcomes and posterior predictive replications") +
  annotate("text", x = 0.2, y = 5, label = "Density of actual outcomes", hjust = 0) +
  annotate("text", x = 0.2, y = 3.5, label = "Posterior replications", colour = "orange", hjust = 0) 

And there we have it. For each value of x, we have generated many plausible outcomes. It is these we should be comparing against the outcome, not the average of the many draws.