In September 2023 I taught a week-long course on statistical workflow at the Nelson Mandela African Institution of Science and Technology (NM-AIST), a public postgraduate research university in Arusha, Tanzania established in 2009.
The course was hosted by Dean Professor Ernest Rashid Mbega and the Africa Centre for Research, Agricultural Advancement, Teaching Excellence and Sustainability (CREATES) through the Leader Professor Hulda Swai and Manager Rose Mosha.
Our case study was an experiment on the NM-AIST campus designed and implemented by Dr Arjun Potter and Charles Luchagula to study the effects of drought, fire, and herbivory on growth of various acacia tree species. The focus was pre-data workflow steps, i.e. experimental design. The goal for the week was to learn some shared statistical language so that scientists can work with statisticians on their research.
Together with Arjun and Charles, with input from Drs Emmanuel Mpolya, Anna Treydte, Andrew Gelman, Michael Betancourt, Avi Feller, Daphna Harel, and Joe Blitzstein, I created course materials full of activities. We asked participants to hand-draw the experimental design and their priors, working together with their teammates. We also did some pencil-and-paper math and some coding in R.
Course participants were students and staff from across NM-AIST. Over the five days, between 15 and 25 participants attended on a given day.
Using the participants’ ecological expertise, we built a model to tell a mathematical story of how acacia tree height could vary by drought, fire, herbivory, species, and plot location. We simulated parameters and data from this model, e.g. beta_fire = rnorm(n = 1, mean = -2, sd = 1) then simulated_data …= rnorm(n, beta_0 + beta_fire*Fire +… beta_block[Block], sd_tree). We then fit the model to the simulated data.
Due to difficulty in manipulating fire, fire was assigned at the block-level, whereas drought and herbivory were assigned at the sub-block level. We saw how this reduced precision in estimating the effect of fire:
We redid the simulation assuming a smaller block effect and saw improved precision. This confirmed the researcher’s intuitions that they need to work hard to reduce the block-to-block differences.
To keep the focus on concepts not code, we only simulated once from the model. A full design analysis would include many simulations from the model. In Section 16.6 of ROS they fix one value for the parameters and simulate multiple datasets. In Gelman and Carlin (2014) they consider a range of plausible parameters using prior information. Betancourt’s workflow simulates parameters from the prior.
Our course evaluation survey was completed by 14 participants. When asked “which parts of the class were most helpful to you to understand the concepts?”, respondents chose instructor explanations, drawings, and activities as more helpful than the R code. However, participants also expressed eagerness to learn R and to analyze the real data in our next course.
The hand-drawn course materials and activities were inspired by Brendan Leonard’s illustrations in Bears Don’t Care About Your Problems and I Hate Running and You Can Too. Brendan wrote me,
I kind of think hand-drawing stuff makes it more fun and also maybe less intimidating?
I agree.
More recently, I have been reading Introduction to Modern Causal Inference by Alejandro Schuler and Mark van der Laan, who say
It’s easy to feel like you don’t belong or aren’t good enough to participate…
yup.
To deal with that problem, the voice we use throughout this book is informal and decidedly nonacademic…Figures are hand-drawn and cartoonish.
I’m excited to return to NM-AIST to continue the workflow steps with the data that Dr Arjun Potter and Charles Luchagula have been collecting. With the real data, we can ask: is our model realistic enough to achieve our scientific goals ?
Shira:
I like hand-drawn sketches too. Not always, though!
Those hand-drawn sketches better not ruin it for the rest of us ! Always good to have counterexamples.
Unrelated question: When doing fake data simulation, do you prefer fixing parameters, looking at range, or simulating from a prior ? Does it depend on whether the purpose is teaching versus research ? In Section 16.6 of ROS you fix one value for the parameters and simulate multiple datasets. In Gelman and Carlin (2014) you consider a range of plausible parameters using prior information. Betancourt’s workflow simulates parameters from the prior.
Shira:
When simulating from the prior you can do simulation-based calibration checking as laid out in Cook et al. (2006) and Talts et al. (2020). You should also take a look at Modràk et al. (2024), which has our latest thoughts on the topic.
I often find it helpful to simulate with parameters (or, for a multilevel model, hyperparameters) set to fixed values, and then you can check that the assumed parameters are approximately recovered. We have some thoughts on this in section 4 of our Bayesian workflow article and in this new paper on posterior recalibration.
I like to start by simulating just once because then I can understand the model better, and when there’s a problem it often shows up right away. Simulation from the prior can make sense when you’re working with a model that is already understood and you want to build some trust in your Bayesian inference computation.
Thanks, Andrew ! I hadn’t seen the more recent papers, really cool.
In terms of SBC(-like) computation checks:
I like thinking about the spectrum from weaker prior you intend to use to fixed truth point, as you write in section 4 of your Bayesian workflow article:
“simulation-based calibration and truth-point benchmarking are two ends of a spectrum.”
Similarly, in your new paper on posterior recalibration you say:
“supposing a particular ‘true’ value of theta could be considered as a special or limiting case of the use of a sharper prior for evaluation than for modeling.”
QUESTION 1: When using a sharper prior for generating the SBC samples, when would you use this sharper prior (6) also in the computed posterior (and therefore get calibration) versus using the weaker prior in computing posterior (4) (and therefore not get calibration) ? Which is what you suggest in section 4 of your Bayesian workflow article when you write “SBC could draw parameters from widened posterior fit to real data” ?
QUESTION 2: In terms of Design Analysis (not computation) checks of what to expect in our experimental design in terms of uncertainty interval width, Type M and S errors, etc.: thoughts on drawing parameters from the prior versus looking at plausible ranges versus fixed values ?
Re: Shira’s QUESTION 1 on when to use sharp/weak prior,
As I think there are great rooms for improvement in dynamic model update, please let me focus on this. I’d use weaker prior if I sense increase of heterogeneity in data generating process (DGP). Two examples: first is when DGP is expected to be vertically differentiated. For instance, if you’re a marketing firm segmenting customers, knowing the fact that GPT diversifies customer desires (I don’t know this is true) and not updating your model would make you feel uncomfortable. However, you’re not ready to add one more layer. In this case, you can have the effect of injecting heterogeneity by widening your prior, until you finally commit to one more layer.
Second is when DGP is expected to be differentiated horizontally. Keeping track of transmission rate and recovery rate of pandemic may be hard due to virus mutation and different biological and behavioral response across countries.
In both cases, the more is at stake in failing to capture these dynamics, the wider prior should be. Marketing firm in luxury brand industry would use wider prior than that in commodity industry. Transmission rate of a Covid model that kills people is more likely to have wider prior than that of chronic wasting disease that harms deer.
If you flip this, if the object you’re trying to model evolves slowly or you have strong control over it, it’s safer to use sharper prior. Andrew may have better suggestion but you may want to skim BDA chp.17. Further thoughts on dynamic model in the context of SBC is in https://discourse.mc-stan.org/t/could-sbc-rank-uniformity-be-attacked-by-unstable-equilibrium-point-in-generator/29960. Happy to get your feedback!
Thanks so much, Angie !
Your thoughts on when to use sharp/weak priors for modeling a DGP are super interesting. I did not intend to ask about this (so it is bonus wisdom !). My focus was only the calibration checks (e.g. SBC) of the Bayesian computation. As Michael Betancourt notes in the Stan discourse discussion you link to, SBC doesn’t require that models are close to the DGP.
In QUESTION 1, I was referring to this new paper on posterior recalibration:
https://www.stat.columbia.edu/~gelman/research/unpublished/calibration_toy.pdf
They describe prior calibration checking as simulating from:
(2) prior: theta^l ~ p(theta)
(3) observation model: y^l | theta^l ~ p(y | theta^l)
(4) posterior, assuming the prior and observation model above: theta_post^s | y^l ~ p_software(theta | y^l)
They note that when (2) is replaced with a new sharper prior (6), but (4) remains the same (i.e. continues to assume the old weaker prior (2)), that we would not get calibration.
But why not also update the computed posterior (4) to assume the new sharper prior (6) ? Then we would get calibration, right ?
This was very interesting! Thanks for the description. Regarding the last bit about “the voice we use throughout this book is informal and decidedly non-academic” – I agree; that is nice approach to writing about statistical modeling. I try to do this when I write descriptions of my workflow for researchers that I do statistical work for.
Thanks, jd ! Can you share any of these descriptions of your workflow ? I’d love to learn from them.
Unfortunately, I cannot share the markdown files, but I will try to describe what I meant by the comment:
As an example, I recently analyzed a laboratory experiment that involved serial dilutions. In the experiment, a laboratory process produces a large quantity of pure spores; repeated samples are taken from that spore stock (or a single sample is taken, diluted and repeated samples taken from those dilutions) and placed onto certain materials; the materials are then submerged in solution and shaken to remove the spores; samples from this solution are taken and then serially diluted many times and then finally put onto growth plates; the number of spores on the growth plates are counted. The goal is to estimate the number of spores in the original stock and the difference between the number placed on and removed from the material.
I wrote the experiment into Stan, using an approach where the samples from the original stock are latent variables in the model. Everything maps from the observed plate counts back to the latent samples. (It’s actually a bit more complex than I describe here, but that seems enough info)
When describing the above model to the PI’s, I decided to write up the report in (what I hope) was a more fun and easy to read format. Some of the section titles include “Dilutions of grandeur”, “Serial Errors”, “A distribution solution”, “Poisson the perfect poison?”, “Thinking forwards”, “Thinking backwards”, “Roadmaps”, “Becoming a Stanimal”, etc. Throughout the report, I try to use a friendly language style that is incorporating – “we”, “us”, etc – and I use conversational language, “Let’s think about our simple serial dilution in the above figure from the perspective of identifying the observed and unobserved variables in our experiment. What do we observe? Only the plate counts of spores at the end of the serial dilutions. What do we not observe? The actual number of spores at any point in the experiment other than those on the growth plates. What do we want to know? The actual number of spores in the stock sample.” When describing technical aspects of the model, I aim for a conceptual approach. For example, when I describe different probability distributions, I write the descriptions and figures in terms of the kinds of numbers these distributions produce given the parameters. I use lots of histograms. Before ever describing a model, I walk through a very simple simulation of some growth plate counts in R, given that we know the latent spore stock (“thinking forwards”). Only after the simulation, do I describe how we might estimate the latent parameters given the observed plate counts (“Thinking backwards”). To me, this is more intuitive, and nothing about the statistical model seems like magic. I try to use very plain, conversational language. I include a lot of pictures. I have a full diagram of the laboratory experiment drawn out as a picture, and each snippet of Stan code is placed next to the part of the picture that it corresponds to.
Admittedly, it would be impossible to give a full non-technical ‘conversational’ explanation of something like a probability density function, but that was not the point of this report. The point was that the researchers clearly see the modeling process, a conceptual of how the model works and how it relates to the experiment, and all of the various assumptions that are made. And all of that in a friendly manner.
Thanks again, jd !
I love your conversational language: the humor and warmth, the clear and motivating questions-and-answers… wonderful.
I also love the forwards (probability) versus backwards (statistics) language.
I especially love putting pictures directly next to Stan code snippets. So much easier than papers where you have to scroll pages and pages from the text to the corresponding figure !
Thanks for sharing this! I’m eager to try something similar in our quantitative methods courses now.
Thanks so much, Shane ! I would love to see your course materials.
I really loved this. Ben Orlin has been doing awesome stuff with hand-drawn math art at https://mathwithbaddrawings.com/, recommended to anyone who dug this post.
Thanks, Skip Intro !
How is this the first time I’m seeing Ben Orlin’s work ?? Just bought his book and started to look thru it. Aside from titling the chapter on statistics “The Fine Art of Honest Lying”, I’m quite liking it. :)