Statistical Rethinking Lecture B09 - Generalized Linear Madness

Okay, welcome back everyone. This is lecture B9 of Cisco rethinking 2026. Uh today I want to take you beyond the generalized linear model as wonderful as it is. Uh before I get into that, again I'm going to remind you deviation from the original plan schedule just for this section. Um your last lecture will be on the 20th of March, not the 13th of March uh because I'm going to be coming back from Cambridge on that day. I apologize this was this was uh poor scheduling on my part. Um but you miss a lecture on Friday the 13th. Yeah, if you're at all superstitious. No. Um okay. So, what is this about? Well, um I know some of you have seen this video before. This is a delightful video. Uh this video is a metaphor for modeling with uh linear and generalized linear models. The the thing about doing scientific modeling with uh linear regression, multiple linear regression, and generalized linear regression, generalized linear mixed models, and and so forth is you can model lots of stuff with them. They're extremely flexible. Absolutely, they're the square hole. You can stick all of the shapes through it. Right, just absolutely every time. But it's disappointing if you're a scientist and you're attending to the scientific details of your work because then when you look at the structure of the statistical model, uh even if it's generative, even if you can simulate predictions from it, it doesn't reflect the structure of the scientific topic you're studying, right? There's this weird gap and that friction there. So, that in the lecture today I want to go beyond this um and stop sticking all the shapes through the square hole. Okay? I'm just going to give you two examples. Uh there's a chapter in the book called generalized linear madness, same as this lecture, which has um two additional examples, I think if I remember it right, uh each of a different flavor, uh but I've only got time for two today. Oh yeah, so it's the the her expressions are just so good. It's like, "No, you monster." >> [laughter] >> Don't do it. Isn't that great? Okay. Um so, that's how I feel when I read journal articles. All right, though. I generalized linear models are great. Okay, GLMs, you can do a lot with them. Uh they're very flexible, and we can scale them up to really big data sets with lots of variable structure. Uh right, then they're like this. All right, this is a generalized linear mixed model. It's still a cat. Still all the same basic machinery, right? It's just bigger. It's a wonderful thing about cats. I'm definitely a cat person, uh is that no matter the size of the cat, they all behave the same. Tigers are just giant house cats. They have all the same play behaviors and grooming behaviors and everything. Isn't that a wonderful thing? Um and and you can also just stick GLMs together in interesting and functional ways, things like the fixed effects strategy, right? Where we put a bunch of them together. This is totally legitimate. There's nothing wrong with it. I searched through a lot of cat photos to make this meme for you. I hope you all appreciate. Um but this is not an argument against GLMs and GLM GLMMs and fixed effects models and the like. It It's actually an argument for them. They're broadly super useful. Uh That's a really great thing about them, but we can do better. Uh uh they're flexible uh machines for measuring partial associations uh combined with structural causal models and generative models and the rules of d-separation, we can make powerful causal inferences with them. Uh they're wonderful tools, essential parts of what we do in the sciences. Uh but there are scientific phenomena which don't fit neatly into these uh single equation model structures. Um and we want to have some way to express them. And when we get more scientific modeling into uh the scaffold that builds the statistics, sometimes we can discover new things that we don't discover otherwise, right? By by starting with something other than the measurement problem, but starting with the causal mechanisms themselves. And I want to show you some examples. So, let's go back to the height data. There going to be two examples in this. If I do the timing right, they'll each will be about half of the lecture. Uh, the first one is the old height-weight data that I started with in the other section, actually. So, at the start of the book, we start by modeling this is the first linear regression example in this section. We didn't start here, but this is this data set that Nancy Howell collected from Kalahari foragers in the 1960s, mostly in 1960s. And this is just human height-weight data. And there's some other anthropometric measurements as well. Um, and I used this in the very beginning of the beginner section, at the very beginning of my book, to introduce the linear model and how we can justify things like normal error and and such. Now, I want to return to this and think about modeling the relationship between height and weight across the human lifespan. Uh, but outside of a regression framework. Really going to the very basics of the geometry of the human body. So, this is allometry. If you're a biologist, that's the magic word that you would use. And I want to convince you that there's a lot to this. In particular, if I can gamify this a bit, I would ask you how many parameters do you think you need to fit this curve? Yeah? I'm going to do it in zero. Well, I need one just for the residual error around the curve. But for fitting the mean, I need zero parameters. Uh, because the allometry that relates height to weight is governed by nothing but the relationship of length to volume. And I want to derive that for you and show that what what how that works. Uh, Okay. All right, this is a human being, right? >> [laughter] >> You all recognize this. This is an Italian human being. >> [snorts] >> But, nevertheless, it's this is a representative human being. Right? And uh I'm going to approximate this human being. It's got some, you know, there's there's a geometry to the human body. It's got volume. There are limbs. Uh the point of this illustration uh when Da Vinci drew it was to study the um the geometry of the human body. Uh I'm going to approximate the geometry of this body with a cylinder. Cuz I'm a scientist. Right? Uh we're going to start here. Obviously, this is not good for all purposes, but I want to show you it gets you quite far just by doing this. And what I mean by approximate it with a cylinder is we're going to focus on two length measurements, uh the height and the radius of the human body. Okay? Uh and there's an equation you learned in secondary school and uh had no reason to use for the rest of your life, probably. Uh that relates volume to um uh that tells you the volume of a cylinder given its height and radius. Uh and it is that the volume is pi * r ^ 2 * height. Um now, this is not a geometry course. >> [laughter] >> Uh but uh if you want to recover the intuition you once had for this, you're just extruding a circle along a height h. So, the area of a circle is pi r ^ 2 and then we just multiply h to extrude it along the height of the uh of the uh of the man. Good? And that's all there is to it. Okay. Let's build up um to the relationship between height and weight. Uh let's take the radius now and relate it to height. So, uh for most human beings, uh there's a the height is a kind of regular proportion. Um uh the their radius of the cylinder, uh the diameter of their body is a regular proportion of their height. This varies by individual and it varies by population, but for a particular individual, we can just uh say that your your diameter or your radius is a proportion of your height. Good? So, I've replaced R with um another symbol, P. Uh this is progress, right? We have the same number of parameters we started with, but it feels like progress. Yeah? Um and then the next thing we need to do is relate uh volume to weight, and we do that with another parameter I'm going to invent called K, which I'm going to call density. It takes a volume and it says how that translates to weight or mass. Uh this can also be a property of the person or the population because your lean muscle mass will affect your density and it affects how much your volume relates to uh to weight, right? This is a thing they do in professional sports when they float the athletes in those tanks and measure their their body fat content, right? It turns out to be like 0.1% or something sickening like that. But uh uh but this is what K is going to measure. Are you with me so far? So, now we we can relate a weight measurement um uh to a height measurement with all these little parameters in between. Um and this is the final equation. Weight is K * pi * P squared, uh where P is a uh the proportion of your height that is your radius uh times your height cubed. Um cubed is nice. You can I hope you feel the tingle of being onto something cuz that's a volume. Right? Yeah, if it was H squared, that would be an area. H cubed is a volume. Now, we are approximating your weight. This is I don't know if it feels like it, but fundamentally, this is just a cube that is your height. Your body is now a cube adjusted by some stuff in front of it. Okay? Um All right. So, let's I should have zoomed in on this so far. So, we've got weight on the left, height on the far right, uh density and proportionality. Pi is not a parameter to be estimated. We know pi, right? Thanks to the ancient Greeks. Yeah, we know pi. Um, but we need to estimate uh K and P. Uh so we we're going to need to make a statistical model now um where uh weight is our outcome variable and height is our predictor. And this is our regression equation. Does it feel good? Yeah? There's no magic has happened here. I've just done the uh violent thing of approximating the human body with a cylinder. Uh but I want to show you this works unreasonably well. Not perfectly, but unreasonably well. Okay, what do we need to do next? We're going to need an error distribution to relate weight uh to the expectation because this equation that we've written is just the expected weight for each height and there's going to be an error distribution around it. So that's one of the things we need to uh decide. I'm going to postpone that choice for a moment and work on the priors first, but then we're going to come back to the error distribution, okay? Um, we're going to need a uh prior distributions for the proportionality uh uh parameter and a prior for the density. Um Choosing priors in generalized linear models is often uncomfortable, right? Because what did those things mean? You've got alpha, you've got beta coefficients on standardized measurements, it's really strange. It's much easier in a scientifically inspired model. These things mean something. And so the scientific background gives you reasonable ranges to start with and we can start with that and then do again, you know, uh everybody in here is an expert in prior predictive simulation, so we can do prior predictive simulations as well and see what we get. So here's my my recipe for setting priors. We're going to choose some measurement scales, and then we're going to simulate uh and then we're going to think, okay, in that order. You might think thinking uh should be number one, but I have a particular thing that I mean by think and you'll see when I get to it. Um Okay. What do I mean by choose measurement scales? Uh I haven't emphasized it, but it's always true that your measurements in regression equations have units on them. Yeah? Unless you've divided them out, which I keep encouraging you to do because units offend God, yeah, so to speak. Maybe not offend him, but they're they're our invention and they get in the way sometimes. Uh but in this case height is measured in centimeters in the data that that I've been giving you. Um and weight is measured in kilograms, and that tells us the units that exist on these parameters. Uh pi and P are unitless. All right, they're ratios. You May you believe me there so far? Proportions have no units. Yeah? Just like probabilities. They have no units. Uh but K is a density, and a density will have units like kilograms divided by a volume. Yeah, so kilograms over cubic centimeters. That that has to be the of K. That has to be true because otherwise you couldn't multiply uh K by uh H cubed and get kilograms on the left. All right, do you remember back to your secondary school education where you had to carry the units through? Yeah? Wasn't that horrible? Yeah, my son is going through this now cuz he's taking physics in in secondary school, and it's it's great fun. Um it's a it's a way to check your math, right? You got to get the units to balance. Uh so it's actually a really nice thing. You can do this with regression equations. All the same things hold there. Uh so in order for the left to be for weight to be on units kilograms, K has to have units kilograms per cubic centimeter. Um Uh Why is this useful for us? Well, um I want to get rid of these measurement scales actually. And that's the first thing I'm going to do. Uh they're like I said, they're an artifice of of the human species, the measurement scale, even as useful as they are. Um So, I'd like to divide them out. And how can we do that? We can divide the measurements, weight and height, by reference values. And that divides out the units. So, if I divide height by the average height, either in the sample or in some general population, doesn't matter. It's just some reference value that you like. Um that'll make that reference value one. And all other heights will be a proportion like that. And the units will be gone. The centimeters will be gone. It'll just It'll have the meaning The measurement will have the meaning proportion of the reference value. But not No information has been destroyed. We can That's That's a non-destructive transformation. We can always reverse it. Yeah? There are destructive transformations, but this is not one of them. Does that make sense? Okay. And we can do the same for weight. So, in this case, to make it simple, I'm going to divide all the heights by the mean height and all the weights by the mean weight. Um And that gets rid of all the units. Uh so, these curves look the same. This is just to show you that there's been no destruction or actual transformation. We've just changed the numbers on the axis and gotten rid of the units. Uh Nevertheless, work has been done here. Because now we're going to be It's going to be easier to think about um the parameters P and K, uh free from these arbitrary things like kilograms. And what is a kilogram? Something some Swiss person made up sometime some hundreds of years ago, right? Um Actually, I'm not sure it was a Swiss person, but I'm going to bet it was a Swiss person. Yeah? I don't know. I've When in doubt, guess a Bernoulli. That's like That's like my rule. Who did it? Some Bernoulli brother, right? There were like nine of them. They did a bunch of things. I'm not sure, though. >> [snorts] >> Okay. You with me so far? Yeah? Does this make sense? Okay. Now, we've still got our priors, um but now we can think about them as a proportion. So, uh P is a proportion, so it must be between 0 and 1. And it must be then less less than a half because the the radius of your body as a cylinder is less than half your height unless you're Gimli. >> [snorts] >> Right? If people still know who Gimli was. Sorry, I should check with these things. It's now It's now a classic movie. Um But, right, for most of us, uh uh it's going to turn out, and this is just a foreshadow, for small children that is not the case, right? Small children are never wider than they are tall. That's not true, but they're they're closer to being as wide as they are tall. And you're going to see that that's going to manifest in the fit. Uh uh which is a general feature I'll return to again when you have a scientifically inspired model like this that you build up to, the ways that the model misfits the data are informative scientifically in a way that it's not in a generalized linear model. Because generalized linear models are epicycles, remember? You can fit anything with them. Uh but the misfit is not a not an indication of physics. Yeah? Just like, you know, the Ptolemaic model of the solar system, uh when the orbits are misfit, it doesn't tell you you need more epicycles, right? It's There's not I mean, more epicycles would help, uh but that's not going to help you figure out the scientific basis. Uh here the misfit is way more informative. Okay, so we can put a prior distribution on P by thinking about what it is. Um K now, uh it's a density, so it needs to be positive, and it needs to be greater than 1. Um because after we we've rescaled, I had it on the previous graph. After we've rescaled the heights and weights, um uh the weight is is larger than the height after we've done the means. Yeah? For individuals. You'll see it again when the data comes back up. So, let's put in some some crude priors here, um which don't satisfy all all constraints. I'm just giving you an example of some weekly informative priors. But in these cases, you can argue for these scientifically based upon the meaning of them, which is not generally easy to do for generalized for for slopes and intercepts. Um So for P, I propose beta 25/50, uh which I know sounds ridiculous and weird, but that's because it creates the distribution on the screen there, uh the the one on top, um where it's all the probability mass is below a half, um but not much below a half. Uh but it's weekly informative. Anything There's a lot of uh positive probability between 0.2 and 0.5. Uh we'll start there and see what the curves look like. Um Uh and then K is exponential 0.5. Uh remember uh or learn, if this is the first time you thought of this, the the relationship between the rate of an exponential and the mean of an exponential is that each is the inverse of the other. So the mean of of an exponential with a rate of 0.5 is two, because that's 1 over 0.5. I know this is somebody made a decision a really long time ago to parameterize it this way, and that's just how it is, but the the the mean of an exponential is the inverse of its rate. Okay? Yeah? So I'm making it two cuz that's greater than one, and the density needs to be greater than one. I'm I've chosen a distribution though which will allow values less than one. This is um a a form of prior choice that I am a fan of. This is what's called weekly informative. There's a constraint, but I'm not making it hard. And I'm not making it hard because that gives me an opportunity to discover model misspecification. If I fit this to the data and K ends up being less than one, that makes me back up and check everything. If I force it to be greater than one, I don't have the opportunity to opportunity to detect my mistake. So that makes sense? This is called weekly informative priors. In the Gelman-Vehtari and me Bayesian workflow book that will come out later than this year, we argue very strongly for this in context of some case studies and such. It's just that you want you want to work in a way that gives you the opportunity to discover previous modeling mistakes. And I think hard boundaries on parameters will stop you from doing that sometimes. So, this is weekly informative. We're leaving information out like this needs to be greater than one, but it's strategic to to help us debug our workflow. Okay? If it if the posterior distribution ends up being greater than one, which it will in this case, yay. But, something else is going to go wrong here. Spoiler alert. Um and then then we'll work forward from there. I should pause for a moment and ask if there are questions. Does this make sense? I haven't done prior predictive simulation yet, but I'm going to. But, to do the prior predictive simulation, I need to pick a distribution, an error distribution as well. So, that's the next thing we're going to do. Remember cuz the prior predictive simulation uses the whole generative model, and I haven't picked the likelihood yet. I haven't picked the error model. So, that's next. Um what do we know about weight? Um it's always a positive real, right? You can't have negative weight. Yeah, agreed? There's the skeptical faces. Just think it through for a moment. Take as much time as you need. Right, it's like a philosophical question. Um Uh and the other thing about measurements like weight is the variance scales with the mean. There's not a constant error around the prediction equation for all weights. You buy that? Yeah, and this is this is common of of almost any growth process in nature. It's not that the scaling will be the same across all, but there's a there's a coefficient of variation if anything is maintained across it that the ratio between the variance and the mean might be maintained. Um might be. Uh but, that's at the best thing that might be maintained. So, there's a very natural distribution that arises from lots of growth processes and that's the log-normal. Um which we haven't worked with a lot in this in this section yet, have we? Did we have another example with log-normal? I don't think so. Log-normal super useful for growth processes. It's super natural. It's really easy to work with, too, because it's just the log of a normal. That's it. There's nothing else to learn. And indeed, it's parameterized like a normal distribution. The The parameters mu and sigma in it are the parameters of the normal that you get when you log the measurement. This is You have to keep this in mind because the mean The mean of the outcome is not mu. The mean of the outcome is a function of mu. And I'll show you what that function is in a moment. It is not simple because the variance scales with the mean. But I just want to get this stuff I've got a little bit of text. Yeah. Mu in log-normal is the mean of log, not mean of observed. Does that make some sense? Yeah. So you If you took all the weight measurements and you logged them, we're saying they'll have a normal distribution. But the raw weight measurements don't have a normal distribution. They have a log-normal distribution. This is not helping. This is I know this It's like inadequate language for this. Um and the mean if we parameterize a log-normal with mu and sigma, which are the parameters of the normal that you get when you log the measurements, that means that the equation for the mean of the raw measurements is a function of both parameters. And I'm going to show it to you in a moment, okay? It's a It's just true. Which means the mean The mean scales with both of them, which is what we want. Right? So this is all good. I know it seems annoying, but it's because it obeys the physics of growth processes. And this is how it works. You with me? Yeah. Um Okay. Yes. All right, here we go. Look, now we're ready to do the prior predictive simulation. Uh How does this work? Let me walk you through the code a little bit. Um I'm going to We're going to draw do 30 draws from the prior. I do 30 draws from the beta with 25 successes and 50 failures. That's what That's what how the beta distribution is parameterized. Yeah. Um and then 30 draws from the exponential for K. Then we've got the sigma parameter, which is the error on the log scale. And I'm going to give that this conventional exponential one that I've been using a lot. Um now I go through each each of the Well, I set up a sequence across the horizontal axis. That's what X sequence is, all the heights I want to simulate the expected weight across. I set up a an empty plot. That's what plot null does. Yeah. I'm rolling with base R here. I know for some of you this is painful. Like, why is old man McElreath doing base R? I'm sorry. This is just how it's going to be. All right. I'm not against GG plot. This is just, you know, this is the way you do this stuff in base R. I learned R in 1999, okay? That's just That's why it's like this. Yeah. So, we're going to party like it's 1999. And then I then I loop through each of the prior draws and I compute the mean. So, the first thing I do is I say mu is the log of um That's the mean equation, right? It's pi times KI cuz each I is one of the prior draws times PI squared. And then X sequence cubed. And X sequence are the height values, right? So, for each of them. Um and then I draw the line that goes with that X sequence on the horizontal. Now, here comes the equation for the mean of a log normal. It is E to the mu plus sigma squared divided by two. This rolls off the tongue, doesn't it? So, the mean of a log normal scales both with sigma and mu as a consequence of that equation. Yeah? And that's that's the thing you can just derive from the relationship. Um and again, this this is a property that we want for growth processes. It makes the mean scale with the variance. Yeah, it's a good thing. Um anyway, then we get the graph on the right um and the the horizontal and vertical dash lines there show you the reference values of one and one. Um so, the point is at the population mean combination of those. This is a very loose set of priors that cross a bunch of growth processes. Yeah? Um And in particular, we know that the the curve we want needs to pass through that point. Right? Because a regression line needs to pass through the person with the average weight and the average height. Or at least you should come pretty close to it. Otherwise, it's a bad regression line. Does that make sense? I'll say it again. A good regression line for most problems is going to come pass right through or come really close to passing through a person who has the average weight and the average height. Yeah? Doesn't have to exactly pass through it for a curvilinear relationship, but it should get pretty close or it's going to do a bad job elsewhere. Does that make sense? So, this is this is a super loose set of priors. Uh but we have a lot of data, so I suggest we go forward. Again, the point of prior simulation is not to pre-fit your your curve. Although in this case, as I'll show you in a moment, we basically can't cuz I told you we're going to go to zero parameters pretty soon. Yes, or one parameter. Sigma will survive. That's the only one that will survive. Okay? Are you with me so far? Is this I know I'm going slow, but I think these workflow steps are are useful. Um Now, we fit the model. Uh you don't have to use Ulam for this. You could use quap. Uh you could use uh you could just use maximum likelihood. Uh this is this is not a difficult fitting problem in in any means. Um Make the data list uh uh in this Ulam model. I don't think there are any surprises. DL norm is just the R name for log normal. Yeah? Um and you'll notice the uh exponential link function on the second line, I I write exp mu, and that's because mu is the mean of the of the normal you get when you log. And so the exponent of mu is the mean of the of the raw measurements. And so that's the relevant link function here. You're used to log links, now we have an exponent link. Right? But it's logical cuz that's what we need for the log normal. Cuz cuz we derived an equation that related weight to height on the measure on the raw scale. But the log normal wants a predictor on the log scale. So the exponent undoes that. Did that help? Does it make sense? Yeah? When you review these later, it'll make sense, I hope. Yeah. And I think that's the only surprise. Yes? Is there anything else here I should talk about? Is that good? Okay, good. I appreciate your nodding. It's the good audience participation. Um And all right, let's do a quick posterior predictive check and see how good this fits. Uh that ain't bad. Yeah? That's pretty good. Uh and all we did was assume a person is a cylinder. You'll notice the children are not predicted very well. Uh uh and I've already hinted why I think that is. Um and it's because they're uh the rela- their radius as a proportion of their height is is systematically different. Yeah. And so if we wanted to improve this model, we would adjust P by H, for example. Yeah? That would be a way forward. That's one of the advantages of scientifically inspired models like this is that the way it misfits the data gives you ideas about how to expand the model to do better. It gives you constructive ideas. Whereas with GLMs, you can just add more epicycles and have an arbitrarily good fit anytime you want. But that isn't going to necessarily help you make a better scientific model. Good? All right. Let me Let me quickly finish this example now. We're about halfway through. Now step three is thinking. You might think we've been thinking so far. Yes, we have, but now we're really going to think. Okay? It's like thinking squared. And so let's come back to our equation for the mean. Um And I want to show you something in the posterior samples now. So if we we take the posterior samples for the model we just fit and we plot K against P squared, this is the posterior distribution. Have you seen anything like this before? This This is all a curve that relates the two together. A curve. These two things share a lot of joint information. If you know one, you know the other. Do you see that? There's a curve. These posterior samples are on a curve. This is not a distribution. Well, it is a distribution. Everything is, but this is a curve. We can derive this curve from the equation on the left of the slide. Now we're going to do that. Um Okay. Because I have uh divided the height and weight by a reference value at the values one, I know that there's a person there. We want the curve to go through it, so we know that one is equal to We know that this is true, that one is equal to K pi P squared one. Yeah, approximately. This is the assertion that for the regression, the trend should go through the average person in both height and weight. And I know the average values, they're one cuz I set them. Okay? You with me so far? Um I can solve that equation for K and it tells me that K is 1 over pi p squared. And now I'm going to plot that equation on the posterior distribution on the right. Do you like that? >> [laughter] >> Yeah? I like that a lot. This is This is my happy place, right? So, we can get rid of parameters here. Uh, what are K and p squared doing? They're just handling measurement scales. Uh, and that's why they're bound together and they have they have so much joint information. Um, they're not they're really nothing that needs to be fit at all. So, now I can take uh, uh, K there. I can replace K with 1 over pi p squared in the original equation and the original equation becomes 1 equals 1 cubed. Which I hope you will agree is true. Isn't it nice when algebra works out that way? Okay. Um, so it turns out we don't need any of those parameters and so I'm going to get rid of them. Here's our new uh, regression equation. Uh, they I just say that e to the mu is h cubed. That's our whole regression. The only parameter that remains is sigma and I'm going to fit it again. How well do you think it'll resemble the previous one? It's the same. >> [laughter] >> Isn't that great? So, this is the magic of allometry. So much of the relationship between a length measurement of an organism and its volume and things related to its volume like its mass are governed fundamentally just by the dimensionality of those measurements. The fact that some of them are cubes of the others and all the other parameters are just to handle measurement scales, which are artifice that people use because they have to invent rulers. Yeah? Um, uh, like in the UK I think people still talk about their weight in stones. Is that true? [laughter] Yeah. Which I adore. Don't get me wrong. I think that's amazing. Um, but it shows how arbitrary that stuff is, right? Uh, Um, anyway, so, uh, this general principle I learned in high school a very long time ago from a physics teacher where you do dimensional analysis or dimensionless analysis, you try to divide out the the scales of measurement because that reduces down to fundamental ratios between theoretical parameters and those are the things you want to understand. And you can get rid of all the measurement scale artifice and simplify the problem. And that's what we've done here in the allometry example. Okay. So, let me try to summarize and then we'll do the second example. Um, most of the relationship between height to weight is just a relationship between length and volume, which in hindsight is maybe obvious, but when you start it's not. Yeah, um, again, uh, changes in body shape are likely to explain the poor fit for children. That's a way to expand the model. Uh, think about relationships to age. Um, and this is one of the big advantages of trying to build up a statistical model from the basic biology or physics of the problem first, not going immediately to linear modeling or generalized linear modeling, which again, I'm not against because it's productive. Uh, but this this approach is ultimately much better. Okay, good? Yeah? No questions? Good. All right, let's do the second one. Second one's completely different. Um, for the second half of this lecture, I want to talk about population dynamics. Uh, there's been lots of requests to talk about time series. Uh, Peter, who's not here today, uh, sorry, I just shamed him on a recording. I'm not trying to shame you, Peter. I'm sure you have a good excuse for not being here. Yeah, I would love again, half of Leipzig is sick today. Uh, he sent me a a very thoughtful email asking me about different ways to model time series and I responded to that and I said, "Well, I'm going to do something for you on Friday." So, here's the recording for you. Um, uh, time series data is really important, too, and those of you who know my professional work, I'm very in much into population dynamics. And uh that's one of the things I work on. So, um one of the things about modeling population dynamics is that the measurements aren't uh I'm going to walk our ways into this generally. So, the measurements you can take of the population are done with error. What we want to model are the true states of the population. And treat the measurements as emissions from those true states. So, we don't It doesn't make sense, even though it can be productive, it doesn't make sense scientifically to treat your previous measurement as a cause of your next measurement. Which is the kind of default thing to do in a time series analysis. Yeah? Do you understand what I'm talking about? Yes, there is You're wincing because You're ruining my life with this. Oh, sorry. I mean, no, I'm I'm helping you solve your your problems. Uh but this is a fundamental problem. So, um like I collaborate a lot with ecologists, and this is a standard problem in ecological time series. Uh there's the true state of the population, but we've only got proxies of the animal densities and their locations. And we're trying to model uh the true dynamics of the population, but all we've got are these uh noisy measurements of it. And we want to build models of that, and that's what um ecological statisticians do. Um so, in this case, I'm going to go back to kind of a a core thing in um population dynamics, and think and model the ecological dynamics between two species. This is a very classic ecology model. No one would fit this model to data these days. We have much better population dynamic models, but it's uh it's a nice toy model, and actually captures a lot of really interesting things about population dynamics. Uh so, the the the um uh estimate is going to be how different species interact, and how do interactions influence population dynamics. Uh trying to It's like a a community ecology uh relationship, and it's a predator-prey example. Uh the the predator here is a lynx. Uh these are Canadian lynx uh uh and hares, uh Canadian hares. Um and yes, there are more mammals in this pop in this community. Uh but this is a real historical data set. Um uh the data are thousands of pelts that were delivered to trading stations, um probably near near present-day Toronto, cuz that's all there is in Canada. Sorry, Canadians. But >> [laughter] >> I'm just inviting hate mail from Canadians. But um uh between 1900 and 1920. So number of pelts is a proxy because it depends upon hunting intensity and lots of other things, right? It's it's measured with error, absolutely. But there are these very intriguing cycles which have inspired generations of ecologists and led to lots of productive work um in population ecology for modeling other data sets, um including exper experimental data sets. You can do experimental versions of this in laboratories with insects. And there's been a number of studies like that which are really slick. Um there's a homework problem in my book where I give you an insect data set and ask you to apply this model to it. Yeah, so um uh but you'll see that what you see in lots of predator-prey systems is the the prey density goes up, the black curve is the hare, and then the uh there's a lagged response, it seems to be. Our eyes really want to see some causal relationship here cuz uh lynx eat hare. Where do lynx come from? They come from other lynx. Yeah? There was a time when biologists didn't know that. Just frame But now we all agree that if you want a lynx, you need a lynx. Yeah? Uh and if you want a hare, you need a hare. That's that's agreed. I'm getting funny looks, but there there were classic experiments to see where maggots came from. Right? They sealed a steak in a vacuum-sealed tube, right? This is a classic that I forget the French biologist who did this, but there was this idea of abiogenesis that maggots would spontaneously emerge from rotting meat. I mean, we you laugh, but No, no, famously in the way from eels. Yes. >> confused about how eels came to be, but I just didn't think it applied to links. Yeah, exactly. Exactly. Well, eels are strange. There was a comment about eels. Eels have a very strange life history. It's fascinating. Uh they're heroes of the life life history Olympics in some ways, but uh anyway, so there were these classic experiments and again, I forget forget which French biologist it was. I think it was a French biologist who literally blew glass around a steak. Right? So, and then watched it rot. And then you get no maggots. But if it's exposed to the air, maggots appear. So, where do they come from? They come from eggs. Right? That's where they come from. Okay. Uh but that said, uh we want to see causation here. We want to see that uh the links are eating the hares and that increases their uh their population, and then they exhaust the the the local hares, and then their populations go back down, and the cycle can continue until one or the other crashes. Yeah, and this is a classic predator-prey system. And there are many models of this, and I want us to work with the most basic one. Um Uh little bit of general information about time series modeling. The standard GLM approach to time series modeling is to use lag variables. So, this is the kind of dag you would draw if you were going to talk about a time series. Yeah, Sarah is nodding. I know you've done this, and I'm I'm a big fan of this. There's nothing wrong with this. Like, diagramming your receptions is good. Um this is what you call a cross-lag model, where the the measurements uh Y and X in the next time period are influenced by the previous measurements and some constant effects like Z, which are influencing everything. Yeah? Um and the lag variables is uh we're going to predict, say, T using Y T minus 1 and X T minus 1. That's the GLM approach. Uh you can do a lot of useful work this way. Um sometimes people will do even deeper lags, which are often necessary to get a good fit, like uh y - 2 and t and x - 2. This is non-causal nonsense, but it will let you fit curves. Right? Why is it non-causal? There I assert that something that happened two time periods ago cannot exert a causal effect on what happens next. Right? The only thing that affects what happens next is the present. You with me? I mean, even lag minus one's a little weird cuz it's only the present, but that's just a measurement issue. We don't have anything closer. Right? So, any of the causal effect of what happened two time periods ago is only transmitted through what happened one time period ago. Right? And then what happened now? Does that make sense? This is the Markovian assumption about how these about physical simulations. Does that make sense? That that real causation is memoryless in that sense. That all that matters is the state of the world now. And the state of the world now is conditioned on all the stuff that happened before that matters for causation. This is a little bit philosophical, but have I have I persuaded you? So, if you want to do lag minus two, fine, but you're giving up a causal inference at that point. Yeah? Um Okay. What's wrong with the lag thing? Well, aside from the weird deeper lag issues, uh it makes very strong assumptions about the functional relationships among the variables. Maybe they're not additive on some latent scale. Yeah? Um ignores measurement error. You can tack measurement error on here. You can make you can expand this DAG in monstrous ways to add the measurement error. Absolutely. That's You could do that. Uh but usually people don't. Um and there's very little science here, mostly statistics. And the whole goal of this is to push more science into it cuz cuz there are real dividends to be had there. So, let's model this uh in the classic population ecology way. Um what we know about closed populations, and again, this is axiomatic, is since hares only come from hares and lynx only come from lynx, the change in the hair population in uh with respect to time, that's what dh over dt means here, um is the births minus the deaths. And the same is true of the links. That's true. This is not an assumption. Well, it's an assumption, but it's true. It's arguably true. Migration would add extra terms to this. So, you can expand this to an open population, but in a closed population without migration, this is all there is. Okay? You with me? Um And now we can start modeling this. So, we've got we're going to need parameters for the birth rate and the death rates, and we can relate those parameters to the other species, and that's the causal hypothesis that we want. Uh so, this immediately becomes what in in GLM terms an interaction model. So, we get a parameter for the birth rate of hares, which is independent of the links because hares don't eat links. Usually. I suppose if they came across a dead one, they would. But uh in the same sense that deer eat baby birds when they come across them, right? They're facultatively carnivorous. But um uh but you know, the hares are eating some some um unconstraining vegetable resource is the idea here. Uh which is a weakness of the model cuz it's not necessarily true. Uh and then their death rate, however, is related to the density of links. There's a predation parameter, mortality parameter m sub h, which is um uh multiplies time the links density, and then their own density, right? Because if there are no hares, then they can't die. All right? So, l * m is the proportion of them is the probability any particular one dies. It's Does that make sense? Yeah? If you don't have any hares, you can't get negative hares. So, that the the physics of this equation guarantee you don't get negative population densities for the hares. Which is already we're better off cuz a lag model will happily produce negative population densities. Yeah? Good? Um same kind of stuff goes on for lynx, but it's flipped uh because the lynx is the predator. We've got um that their birth rate depends upon the number of hares uh times their own population size because only lynx can produce lynx. If there are no lynx, then none will be born. It doesn't matter how many hares there are. Yes? So, there's an interaction term here, but there's no main effect. And that's one of the things I want to get across to you. If you made a generalized lag model, there would be main effects of of of the um lynx population on their birth rate, which makes no biological sense. Yeah? Okay. I need We got 15 minutes. We'll finish this. All right. Now, we're going to add the measurement layer. Um I apologize this goes a bit fast. This is in the chapter in the book and I go a bit slower with it there. Um we're going to add log normal error onto it again uh because what we have measured are numbers of pelts. And again, it's as population related so the mean and the variance are related. So, we're going to use the log normal as a default. Um and the mean of the mean value is the log of the um uh of the hare population times some parameter P, which is the trapping intensity. Yeah? That's how much they're interested in trapping us. And I'm going to make this a constant here. If you really wanted to put a you know, devote your life to modeling this time series, you would try to figure out the amount of trapping intensity at different years. Yeah, because maybe there'll be proxies for that. Um and then uh uh wait, I should put this. Yeah. And then the observed lynx pelts on the side. Same idea, but uh with little L's for everything. Um now, at the bottom of this, we need to generate predictions for each year after the initial one. And the way to do this is to take those equations we wrote, the differential equations is what they're called. These are ordinary differential equations. Um the dH over dT and the dL over dT and sum them up uh through the time we need to get to the later time point. And when you sum up a a differential equation, you do something called an integral. Yeah? And some of you had integral calculus so and then forgotten it. But uh you remember doing this. Yeah? Um uh In this case, Stan is going to do this integration for us. So, we don't have to do it inside the Markov chain. It's going to do the integration. But the the logical structure of it is what's important to understand here is that we can get a prediction for any year after the first year of observation um by taking that first observation, h sub one here, we can get a prediction for h sub t by taking h sub one and adding adding up all the incremental changes that happen after it according to the model. Yeah? For any particular parameter values. And then we're going to use the Markov chain to get a distribution of those parameter values that fit the data. Just as usual. And this is a statistical model. Now, does it feel like one? No. Um but it absolutely is. Are you good with me? Yeah? I'm trying to read the facial expressions and they are every range of human emotion at this moment. Okay. Um Yeah, I I have an odd education in statistics cuz I started here in evolution in the evolutionary ecology theory and learned statistics after the fact. So, for me, it's very natural to do the something like this on top of a a non-statistical model. Uh but I I definitely appreciate that for most people it doesn't work that way. But I I learned evolutionary ecology uh population genetics before I learned statistics. And so, for me, it's it's very comforting to start with the the the ODEs. Um let's do prior predictive simulation. I'm not going to go through talking about all the priors. I spend a lot more time in the text on this. I encourage you to look at it. Uh it's hard to visualize the prior predictive simulation here because these curves are related to one another in the priors in powerful ways. It generates stronger um uh echoes uh and predator-prey dynamics depending upon the uh the birth and death rate and the predation parameters and so on. So, I'm just here taking some samples for you from it to show you that. Uh but in the book, I build them up. Like, what do the priors need to be like in order to produce uh cycles and prevent things from crashing and and and those sorts of things. In these nonlinear models, prior predictive simulation is really critical because you just there's no way to intuit how these parameters interact when you choose the values, okay? Um okay. This is not something um you can do in Ulam. Or I shouldn't say that. I bet you could with enough of custom injection of code, but we're not going to. I'm just going to do this as a raw Stan model, and I'm not going to explain every detail of it. Um but let me if you'll indulge me, I've got 10 more minutes. I want to walk through this for you and then show you what the fit looks like, okay? So, a Stan model has blocks and the the top block, which is optional, is called functions, where you can write arbitrary functions that do sub calculations for your model. And when you model ODEs in Stan, you use the function block to write a function that just calculates those uses those differential equations. We just write them into this function. So, that's what we've got here. Um there's a there's a header there that just passes in all the parameters and the data. Um and then the differential equations appear there on the lines near the end. You can recognize them, right? I factored out the H and L so they look a little bit simpler cuz I want to do less multiplication, so I always factor to the simplest form when I do math inside the computer. Yeah, if you were operations is good. Yeah. Um but that's just the differential equations. Good? Uh and then Oh, wait. Yeah, I actually Here we go. So, I had another animation where I say those are the differential equations. All right. Good. Um there's a transform parameters block where uh we're going to pass that function to a specialized function inside Stan which uh solves the differential equation set. And this is integrate ODE RK45 here uh is what it is. It's It's a way uh it's just one particular algorithm for doing the that integral. Um and then you pass in There's a particular syntax here uh which is explained in the Stan manual and you pass it all in and it does the hard work for you or magically. It's nice and comforting. Um every probabilistic programming language has functions like this for integrating ODEs. This is just such a common thing to do. Yeah, it's a big area of statistical modeling to do this. Um and that's the good news. You don't have to invent any of it. Uh uh good applied mathematicians have worked hard on those algorithms, tested them, and so on. Uh then finally, the actual model block is almost trivial. You define the priors and then you just um cuz in transform parameters, we've already got predictions for every uh every set of measurements in the whole time series, the all the measurements of the pelts all the way through. That's what the transform parameters block did by passing the integrate ODE cuz it gives us predictions for the whole time series. And so, the only thing that's left to do in the model block is deal with the measurement error. Right? Take the the the integration of the ODEs gives you predicted true population states and then we just have to translate that to the pelts. And that's what the log normal does. And that's it. Easy, right? Um uh this generalizes to bigger systems of ODEs. This is a huge area of modeling. There's this whole area of statistical modeling called pharmacokinetics where people study how drugs and hormones diffuse through the body. Right? Depending upon to reach target tissues. And that's a really essential way of doing drug design and understanding pain relief and therapies and like in chemotherapy and other things like that. And it's it's done with ODEs. You treat the body like a set of compartments with tissues and there are flow rates. That's a huge active area of applied research. And lots of people use Stan to do it for these reasons. You with me so far? Yeah, get the idea? Um Okay. So, we can do posterior simulations now. This is the figure from the book. It's not a great figure. I'm going to explain why and then I'm going to do something better. Um Each of these curves is a sample from the posterior and we've drawn it up. But there's a black curve that goes with each red curve. Right? Cuz they're common draws. Yeah, but you can't see which ones go with which. Um that's just very difficult to do. Yeah, but that's how it is. But this gives you an idea of this this is true population dynamics plus measurement error, which is why there's lots of spiking and jumping because the measurement error here is quite substantial the model thinks. Yeah? Um but you get an idea how it fits. And the points are the raw data. Um Uh so, I'm going to do the thing I do always in these cases and do an animation. Uh which also may not be super satisfying, but at least it's nice to look at. So, these are now tracking one another. They're common draws from the posterior. And you can see like one goes up, the other goes up. Sort of thing that you couldn't see in the other one. Yeah? But there's measurement error here, too. So, you can get them out of sync sometimes. Is it satisfying? Yes, you're loving this? Yeah, okay. Uh Yeah, what once I wrote the the generalized functions to make these animations, I kind of went silly and just started doing it for lots of stuff because it became easy, right? Um Yeah, there's a bunch of undocumented functions in the rethinking package for making these animations now that I use to make them. I should document them. Okay. Um and we can compare this to the prior just to show you that uh the posterior is a lot more contracted than the prior. The the prior relationships are much more general. Um many more flat relationships in the prior, for example. Uh and even stronger ones compared to what we've learned from it. Okay. I'm I'm kind of on time. That's amazing. Okay, what do I want you to get from this lesson? Um well, I wanted you to have an example of a workflow that includes uh uh ordinary differential equations as a way to model um time series dynamics. Uh ODEs, of course, aren't only uh constrained to time. Uh you can have ODEs in space. Yeah, so it's like whatever you want to do. Um Uh this is the beginning of a scaffold up in a project like this, not the end. Real ecologies are much more complex. They're uh links eat other things than hair. So, there are confounds here. Confounding has not gone away, but you can put them into it. You can add unmeasured populations or unmeasured species into these things. You can have shadow ODEs that are unmeasured. And you can generate the consequences of that just like we do with other models. Yeah, so now all of that anxiety about unmeasured confounding is not to be ignored now, but it you have a productive scientific scaffold to build it into as well. Um Uh this kind of stuff is still much simpler than real ecologies, but it puts a lots of great work um in applied management. Uh fisheries management, for example, has uh, going on more than 100 years now of effective use of modeling frameworks like this to really manage fisheries. Which isn't to say that they always do it well, uh, or the way that people want it, but there's lots of great success stories in this. And fisheries are a case where you can never see the population cuz they're underwater. Right? So, you're always dealing with massive measurement error from well, landings on on boats. Right? And counting fish. Nevertheless, they the fisheries people got serious about this, as I said, 100 years ago. Uh, and fisheries ecology is very mathematical for this reason, and it's a whole professional system where every uh, uh, every political entity that has a coast is got to employ some fisheries manager who has a strong training in this tradition. And there are whole handbooks on how to do fisheries management and things like that. And it's an it's an example of this framework. Um, uh, and it depends critically, sorry to go off on fisheries for a bit. It depends critically on the biology of the species. So, for example, lobsters, um, uh, uh, lobsters, you can harvest every adult lobster every year and not deplete the lobster fishery. Why? Because most of the lobster population exists as larva at any one single point in time. Right? This is true for most crustaceans. Most of the population are larva who have not settled out and developed shells. Right? And so, you harvest all the adults after they after the breeding time. Uh, let them lay their eggs. Yeah, have them fun lay their eggs. Harvest them all, eat them eat them all. This happens in Maine. Straight all the lobsters are harvested every year. Next year there will be lobsters. Because most of the population, like there's an order of magnitude more larva than there are adults. And this is true of clams and lots of other stuff in the ocean. Their their life histories are very strange and they spend almost their entire lives as larva floating around on the currents. Yeah? Is this news to some of you? But I think it's a fascinating thing about their life history. This is not true of cod, tragically. Uh and some of you know, the Europeans here know the the the tragic history of the North Atlantic cod. Um they have never recovered, right? Now, the North Atlantic is just stocked full of squid. Which I'm sorry are do not taste as good. I like squid, but they're just not as good as cod. Right? It's just not the same. Um but that's part of what goes on. You need to do You need to get more biology in here to deal with that sort of stuff. Um the life history stages uh uh uh you can put life history stages of the fish into these models as well, and people do that uh lots and all that helps with adaptive management of the population. Okay, that's my sermon about how this stuff is actually useful. Um Right. So, uh I've wanted to give you two examples of of putting things in the in the uh non-square hole all the time, uh and this is the sequel, right? I don't know if you've seen it. Um she's fearful now, right? So, that there is a big applied stats literature where people make bespoke models, scientifically inspired bespoke models, uh and and fit them to data, and there are real dividends to doing that. She says like, "Oh, it's so good." >> [laughter] >> Right? And uh uh uh In the book, I give you two other examples, uh another ODE example, which is a growth process example, and um an example where there are hidden latent states. It's uh uh it's a hidden Markov model, a simple hidden model uh in the chapter as well. So, you take a look at those. Um you can use those kinds of case studies to uh scaffold up to your problems as well. Right, this is we can get the elation here. She's going to be really happy in a moment. Yes. Okay. >> [snorts] >> Very good. Um Right. I know this this material is is difficult. It's the kind of stuff you really need to work through the case studies a bit and vary the assumptions and do simulations to understand. That's really the only way to get there. The best I can do in these lectures is give you a workflow to follow where I think you can responsibly engineer what you're doing and try to inspire you that there are real gains that come from departing from the generalized linear modeling framework. So, thanks for your attention and I wish you well in trying to adapt these tools to to your work.