Statistical Rethinking Lecture B08 - Missing and censored data

All right, welcome back brave people. Uh we are in the eighth lecture of the experienced section of statistical rethinking 2026. We're going to talk about missing data and a particular kind of missing data, censored data, today. This is a continuation of the themes from the previous lecture where we talked about measurement error. Um I want to remind you that this these lectures, too, uh is case where we see the the value really of the Bayesian approach. Uh these are measurement error and missing data are routine problems in almost all realistic data sets. Um unless you are so lucky to work with the perfect experiment. Um you have missing data, you have measurement error. And these are things that are actually quite difficult to deal with in other frameworks, uh but there are very direct solutions in the Bayesian approach and that's nice. Um I'm This also just want to remind you quickly that in this section, the last lecture will be on the 20th, not on the uh uh what was it? 13th? I think it was going to be cuz I will be coming back from Cambridge. I apologize. Yeah, but I I will give you your final lecture. I promise. Okay. Um let's start with a meme. Uh a colleague of mine has a student who made this meme. This colleague teaches with my book and lectures. Uh this really she is one of the things I always say in this course. Maybe you've heard me say it in the book in various places is uh if you're feeling confused, it's just because you're paying it you're trying to understand or it is because you're paying attention or sometimes I say it like uh confusion is what thinking feels like. Right? And I believe those things, but I understand that they're not always super helpful. So, I also try to like walk the walk for you and get there, but uh this is just to remind you if you are feeling that way, you're not drowning. >> [laughter] >> Okay? Um you're not actually drowning. It's just that's what swimming feels like. Okay? So, it's uh we'll do okay. All right. Uh Again, continuing the themes from the previous week, we talked about measurement models as um an addendum or companion, that's a better word, to the generative model. Uh and uh we're going to focus on missingness and censoring today. Censoring is a particular type of missingness where we have partial information about the missing values. Uh and it this is a very routine um So, missing data is super common. Uh I want you to kind of flip the script here and think about most data's missing. Most of the time, right? Um predictions are a kind of missing data. Yeah? That's one way to think about it. Uh and retrodictions are another kind of missing data, but if you have a generative model, you can um develop um uh probability distributions for these things. Um and uh uh Often, we have partial information about missing data, like we know their constraints because we know the type of measurement even if we haven't observed it yes yet. Um uh You know this is true in a sense because you can you can choose likelihood distributions for data you haven't seen. And how do you do that? You use the constraints on the measurement, right? So, I tell you something's a count variable, you're going to say, "Oh, I could use a Poisson." Right? Now, where did that knowledge come from? It's not the values, right? It's missing data, but you know the constraints on the variable and that lets you put a probability distribution on them. And then from that probability distribution, you have a distribution of expectations for the missing values and that's the key insight. Is that um the probability model or likelihood is a prior for the observations. So, when the observations are observed, we call it the likelihood and when they're not, it's a prior. Yeah? I'll say that again. The likelihood is a prior for the observation or prior for the residuals, right? The conditional observations uh in the Bayesian perspective. Um Uh and then the second thing that we know about missing data is from the generative model, the missing values will have relationships with the other variables. So, if there are associations between things and you have a missing value, if you know the other values, then you can update the expectations for the missing values, make predictions about them from the structural model. Uh Okay. Uh as you know, I'm really into workflow, so we want to talk about the workflow for making decisions about missing data. In the first part of this lecture, I want to walk you through an abstract example of the different kinds of missingness and relate them to DAGs. Uh uh as an example to let you think through. And the reason is often missing data is ignorable. Yeah? Have you ever done this? Uh R will do it for you automatically. It will just It will just automatically ignore missing data. It'll drop all the cases. It does by default a thing called complete case analysis. That's often very justifiable, actually, but sometimes it's not. And so, I want to give you a principled workflow for deciding that. And as you guessed it, it'll be based upon a causal model of your variables. All right, there's no inference without assumptions about causation. Um but the good news is often, dropping cases with missing values is entirely justifiable. Uh it's perfectly fine, but you should justify it. Um sometimes it's not and we have to do something else like imputation. Um Uh and imputation uh sometimes it's not necessary, but it's beneficial. It'll give you more statistical efficiency cuz you get to use um you're not throwing away any data. You're using all the In a complete case analysis, you're throwing away things you have measured, right? Because not all the columns are missing. Uh if you do some imputation of the missing values, you can use the things you have observed. And sometimes imputation is necessary to close back doors, which create systematic bias in the estimates. Um as as in the censoring example I'm going to talk about in the second half of the lecture, right? Where it's not ignorable. The missingness is never ignorable in a survival analysis. Yeah? And I'll explain why using the abstract example in the next slides. Okay. Let's let's start with um a simple example. We imagine we have some students and they have produced some homework. Right? Um Let's say some of the homework goes missing. It's not turned in. And why? Because the students say that their dog ate the homework. Um Here's our DAG. Uh Let me reference uh the S is the students uh some feature of the student like their ability. We're trying to assess it through their homework. Uh we don't get to observe all the homework because it hasn't all been turned in. Instead, you observe a vector, that's what the the circle around it means, right? It's missing, it's unobserved. We don't have the full vector of homeworks. What we have instead is H star, which is a partial vector. It's a mix of observed values and missing values. Yeah? That's what we get quite often in real data. Um and the reason is that the students claim that dogs ate their homework. All right? So, now just think about your own research area and translate dog into some other missingness mechanism. But uh the I said I told you it would be an abstract example. You can map this onto your own ideas. Okay. Um Now, how do you decide if we want to estimate the um abilities of the students for whom we do have homework, is it fine to do that with the evidence we have at hand here? Uh now, your intuition says, "Yes, of course it is." But we we we we want to go beyond intuition, right? We want to justify it. And in this case, you know the rules of d-separation, you are pros. At d-separation, right? Yes? No? The faces are saying no. >> [laughter] >> Um we'll work on that. Yes. You are. You are pros at d-separation. You can do this. You can do the backdoor criterion with your eyeballs, right? And so, that's what we're going to do here. Are there any backdoor paths? Um We we want to estimate uh um the the uh effect of the student abilities on the quality of their homework. So, we look for biasing paths connecting H to S. And there are no backdoor paths in this DAG. So, you're good. The H star is a proxy of H. Um But there are no there is no biasing paths in this. So, you're in good shape. This is a case where dropping the missing values and just doing complete case analysis is perfectly fine. You understand? Yeah? In a sense, it's because the dogs are eating homework at random. We haven't measured anything about the dogs in this DAG that guides them to eat some homework instead of other. This is where I'm going next. Does that make sense? It's just it's some dogs has nothing to do with the features of the students. It's just totally random. There are some dogs that eat homework. Yeah, but it has nothing to do with the features of the students, which is what we're trying to estimate. Um but that's not that isn't necessarily the case. You can guess where I'm going now. Let's go back to our original DAG and let's make it so that uh it is the features of the students like how much they study, which influences the dogs eating homework. Like maybe this could go either direction, whichever sort of fantasy story you prefer. It could be that the really good students ignore their dogs and don't take them for walks and the dogs get angry and they eat homework. Right? Um that'll mislead you in one direction. Uh or could be that um could be the other direction, that the students that don't study uh I don't know. I can't come up with a version here. There's something about that which which uh influences uh the dogs' behavior and makes them hungry for paper. Um Play with the dog. Like play with the dog. Yeah. Exactly. Yeah. Thank you. Play with the dog with the homework. Exactly. Uh what's important in this DAG is it's it's not the homework itself that's influencing the dog. We're going to do that next. But it's not the homework itself or the quality of the homework. That's where we're going next. Okay. This situation is has a backdoor path. Do you see it? Yeah? Again, you're pros at d- d-separation. You can see it. There's a backdoor path now. We're going to need to close that, but we can. Yeah? Um Uh if we can if if we uh uh can close that backdoor path by conditioning on on um uh uh things that influence the dog. Yeah, the the missingness example. So, this is not hopeless. Um what is hopeless? Let's go back to the original DAG. And imagine now that that what influences the dog, what influences missingness, is the is the homework itself, the features of the homework. How would this happen? Students who do bad homework feed it to their dogs. Yeah, this is what the teachers think every time they hear dog ate my homework, right? Um if that's true, this is really bad situation for inference. Right? When the outcome of interest is influencing the missingness, um inference is extremely challenging. It's not hopeless, though. If you have a mechanistic model of how that works, you can do inference sometimes. And in fact, the example I'm going to do today is an example of that case. I chose a worst-case example to show you. That worst-case example is survival analysis. Uh but in survival analysis, we can actually save it uh because we have a model of how the missingness works, at least a sufficient model of how the missingness works. Um um Yes, so wait wait though, I should have done the the summary here, biasing path. Yeah, dog eats conditional on homework itself. This is this is the bad case. Um so, here's a summary slide. Um number one, you could say dogs eat at random, to completely at random. Uh in in a real data analysis case, this would be that the missingness is just randomly sprinkled through your table. Yeah, maybe it's uh your research assistant um dropped some things one day, and so there's a gap somewhere, but it has nothing to do with the values. That's what's important. It has nothing to do with the measurements themselves. The process is not guided by the things you've measured, or by um unmeasured things in the DAG, right? Like unmeasured confounds. Okay? That's what's important. Um sometimes that's credible. It's not impossible. I don't want to demean that option. It it it's very credible sometimes. Like I I have a collaboration right now where that's clearly the case. Cuz there was a particular It's a time series study, and there's a particular time where some of the people just took the wrong measurements. Yeah, it's uh soil depth. And there was about a year period where they measured it all at the wrong depth. That sucks. Right? [laughter] But But it has nothing to do with the measurements themselves, or the outcomes, or any of that. It's just somebody not following instructions. Yeah? Um the second case being dog eats conditional on the cause. Uh uh you can if you can condition on the cause, then you um uh you can rescue yourself in this case. Uh and then the final case, the worst-case scenario, uh dog eats the homework itself, depending upon the features of the homework itself. And I say this is usually hopeless, unless you can model the dog. And this is what we're going to do in survival analysis. Um Is this okay? Does this make some sense? All right, this These examples are in the book, too. And in the book, I simulate from these DAGs, and run regressions, and show you the impact of the missingness scenarios. So, I encourage you to take a look at that, okay? I don't have time in this lecture to go through all that, but uh that's what I do in the book. Um Okay. So, now we're going to switch to the idea of when you need to deal with missingness, you're not just going to ignore it, what do you do? And as Bayesians, we do Bayesian imputation. Um Bayesian imputation is not multiple imputation. You may have heard of that. Multiple imputation is a frequentist version of Bayesian imputation. It was developed by Bayesian, Don Rubin, uh as a quick way to do Bayesian imputation uh before Markov chains were common. Uh but it's the theory that justifies multiple imputation is Bayesian. Um Bayesian imputation is the real thing. Now we have Markov chains, we have the full power. Right? We just go to work on it. Uh so, it often in the uh these cases, we will we may need to impute just for efficiency. Um when we are going to deal with uh missingness, we've got two computational options, two strategic modeling options. And I'm going to show you to both of them uh in the in the context of a survival analysis, uh because they're both useful. Um the first is is Bayesian imputation, although I'm going to show you this one second. Uh this is when we we compute posterior distributions for the missing values. We replace the missing value with a parameter. Why? I'm going to say this multiple times in this lecture, because in the Bayesian framework, um uh uh the only difference between data and parameter is that one's been observed. Yeah? Let's say that's the only difference. There are calculation differences, sure. But that's the philosophical difference. And so, the same probability uh definition uh for a variable, say for your observed data, also applies when they're unobserved. And in the first case, when you've observed the data point, we call it a likelihood, or the data model. Uh when it hasn't been observed, it's a prior, but it's the same probability expression. Why? Cuz it has to be. It's the same generative model. Just because you haven't observed something, that doesn't change the generative model. Yeah? Again, follow the axioms. The axioms will not betray you. Yeah, they have no strategic agency. >> [laughter] >> They're not out to get you. The universe might be, but the axioms are your friends. They're reliable. Uh they may not give you what you want, but they'll give you what you deserve. Right? So, uh the second case, and we we very often do this, I do this quite a lot, is marginalization of the unknowns, of the missing values. And this is where we we like average over the distribution of the missing values, so we don't have to compute them. And this is computationally much more expedient when you don't care about the missing values themselves. Uh um this is not always better, though. So, I'm going to show you both, and talk about the tradeoffs. And And there are many cases where doing the imputation is better. Uh and actually uh the marginalization is sometimes prohibitively difficult for numerical reasons. Uh and the imputation is nearly always possible, but slow. Right? So, these are your tradeoffs. Um but I'm I'm going to have a summary of that uh towards the end. Okay. Um Right, I think I already said this. Yeah, cause once you have a causal model of all the variables, um uh including the missing ones, uh this implies directly a strategy for imputation. There are there are computational obstacles, things can be slow, it can be hard to to make the coding efficient. Um you have to manage a bunch of index variables. Uh all of that is kind of a nightmare sometimes, but that's in your software engineering. And we can do that, we can conquer that. We do it incrementally, and we test, right? We know how to do that. Um and sometimes imputation is is unnecessary and difficult for some algorithms, in particular discrete parameters when they're missing. Um Hamiltonian Monte Carlo will not sample discrete parameters. Uh uh but there is no Markov chain algorithm which efficiently samples discrete parameters, because the jumps can be very big. And it's very hard for Markov chains to do efficient sampling when jumps are big. Um but we can marginalize discrete parameters, and then impute them post-sampling. Uh and that's that's what we often do. There's an example of that in the book. Um Uh and sometimes imputation is actually easier than marginalization uh to set up when you're developing the model. As in particular, the example I'm going to show you with censored observations, the imputation approach is direct. Like it's just like the same model. Like almost nothing changes. Whereas, you have to work muck around with a different uh data model for the censored observations when you do marginalization. So, I'm going to show you that as we go. Um The example is, of course, cats. Uh oh, by the way, yeah, I love this cat picture. Don't you Isn't this just like the best cat? >> [laughter] >> I love Anyway. So, uh everybody who has followed me for long knows that I adore cats. So, here's the cat example. Uh these This is a data set I downloaded from the Austin, Texas Animal Control Center uh some years ago now. I should update this data set. It's much bigger now, I guess. Uh but it used to be the case that you go to their website, and they had every animal that was admitted and released. Um and it's not just cats. There's a whole Noah's ark of stuff, uh lizards, everything uh that they're that they're dealing with there. Uh skunks, just like all kinds of stuff. So, I have downloaded 22,356 cats that were admitted. And what we're going to be analyzing And I've pared down the data set. The data set has a ton of columns, right? This is a very big uh thing. Um we're going to be analyzing the number of days until adoption. And the question is whether black cats are discriminated against. So, as you'll know, in many societies, black cats are considered unlucky. Yeah? But from natural selection's perspective, they're the optimal coat color. Right? So, I'm going to say I'm not sure that's true, but that's in the urban setting, probably it is, right? Any other color is is less camouflage. And uh um so, let's ask empirically, uh do black cats in this partic- in Austin, Texas Animal Control wait longer to be adopted? Um and there are many censored records here. What that means is there are a lot of cats, both black cats and other coat colors of cats, that at the time when I downloaded this database, had not yet been adopted. But that doesn't mean they won't be adopted. And it definitely doesn't mean that they were adopted. It's a missing value, cuz the value of interest is how long they wait until they're adopted. Does that make sense? Um this is non-ignorable missing data, because the missingness is caused by the outcome. So, you have to imagine counterfactually, we could time travel into the future, and know when this cat would be adopted. The reason we didn't observe it adopted in the past is because of the value of the time when it was adopted. Does that make sense? Yeah? So, good times, right? So, this is the bad case. This is the bad place. But, mechanistically we can model this and save the day. In fact, you have to. You cannot ignore these values and I'm going to show you that you get a biased estimate. Sometimes a really like twofold bias in the rate of estimates of of adoptions. Um Okay. Oh, I want to say yes, this this is the placeholder cover for a book I have just finished with uh with Andrew Gelman and Aki Vehtari and a whole slew of of co-authors who have contributed a bunch of stuff on Bayesian workflow and this black cat adoption is a case study in the second. The second half of this book is case studies of Bayesian workflow and data analysis and this is in there. So, um Uh and this will be this book will be freely available as a PDF online uh in I don't know, hopefully 3 months, 4 months, something like that. And hopefully we'll have a cover that looks a little nicer at that point. But, it's done. I swear this is done. Feels good to be free. >> [laughter] >> No, it's This is a fun project anyway. I'd say that it's The second half of the book is all case studies. That's There's a huge amount of value in just going through these case studies and Aki is going to have a website up with all the code uh as well. Okay. Um here's part of the data. Remember, there's 22,000 more than 22,000 cat records uh in this and they come in lots of colors. Yeah, this is the best censored cat image I could find on on image search. This is very good. So, so um Uh so, what I've done here is I think I just randomly sampled 100 cats. Uh each each row here is a cat and the length of the line is the days until adoption. The observed adoptions uh occur where there's a point on the end, a little circle on the end. And when there's not a circle on the end, that's a censored value. So, I show you two there. Black here means a black cat and orange means every other color. It doesn't mean orange. Uh there are orange cats in this data set. The The column for for coat color is got all kinds of crazy stuff in it. Like I had to actually just figure out all the stuff I would collapse to black because you know, it's somebody gets the cat and they have to fill in the data form and they look at the thing and they just type something in, right? It's This is the fun of data science. Um so, there's like smoky black, right? So, I had to collapse I made a bunch of decisions that you aren't seeing here, uh but I can share the the original data with you sometime if you're a masochist. Okay. So, this is what the data look like and we want to use these data to estimate rates of adoption and we want to compare the black cats to all the other coat colors. That's the goal. Um Uh of course, there's a workflow. Um what we've talked about so far uh uh is the general phenomenon. We've got a basic estimate here. We're going to need a generative model uh of this. We're going to build this up or I'm going to build this up. Um And then the censoring model, which is the measurement model. Uh and I'm going to of course do the proper software engineering thing and incrementally develop the features for you. We're going to go through three different statistical models as we incrementally develop this. Okay? And then we'll we'll um push predictions out to actually in each step. I'm also going to do some prior predictives in here. I should have put that on this workflow. Okay. Uh let's do a generative model. We're just going to simulate some cat adoptions. Uh Now, I love this example because now we get to a generative model which is not a linear model, right? Cuz we've got a mechanistic process here. We're going to do it day by day. It's discrete rather than continuous time because that's the way the data set is set up. Right? Obviously, there's a time of day the cat was adopted, but we don't have to. We just have the day. So, I'm going to do a generative model that's by day. Um so, it's discrete time survival analysis. The two functions I've written here, the first one called cat adopt, is a recursive function for just simulating cat adoptions. And the idea is there's a constant probability on any given day the cat is adopted depending upon its coat color. I have um the inputs here is is uh to this function is just what day we're on so far and uh that the probability of adoption. Um I've got a flag in here if the day is greater than 1,000 to stop, right? Just in case it runs forever. Um but that doesn't occur as long as the probability is is reasonable sized. Uh and then I generate a random number and if it's greater than probably of adoption, um then uh uh Wait, how's this work? Yeah, then we we keep the process going um and otherwise uh they're adopted and we return the day. So, that's a recursive function. Does it Some people Some of you do a lot of programming recognize this idea. There's a probability distribution for that called the geometric distribution and there's a function in R for that which would which would replace this function completely. But, the I like the recursive function cuz you can feel the generativeness of it, right? You see the process. But, uh there's an analytical expression for the distribution this function produces. But, I I wrote the recursive one cuz I thought it would be cooler cuz you could see the causality better, right? Instead of me just asserting this has a geometric distribution and I leave it as an exercise to the student. But, uh this distribution this produces is called the geometric distribution. Um and then the second function, simcats1, there will be simcats2 and simcats3, right? Uh yesterday in my department we talked about uh Lego Batman 3, right? I discovered there are at least three Lego Batman. I don't know how many there are, but so there will be three simcats. Um And what this does is it uses the cat adopt function to simulate a bunch of cats, however many you ask for, by default 10, and you pass in this vector of two probabilities, one for black cats and one for all the other coat colors, a probability of adoption. Um and it just simulates um a big vector of days until adoption with no censoring yet. There's no censoring, right? This is the pure generative model assuming all cats are eventually adopted. Oh, I should have said in the real data, uh some of the censoring is because the cat is still in animal control at the time I downloaded the data set. Other censoring is because the cat dies, right? Old cats can pass away uh in animal control. And there's also a code called escaped and missing. Cats, they can teleport, right? Anybody who's had a cat, you know this is true. If I bring raw chicken into my house, my cat immediately teleports to my feet. Right? I've seen it happen. Just like sparkles, he's there. So, some of some of these cats have teleported out of the Austin Animal Control and it's a code in the data set. From our perspective, I'm going to treat all those forms of censoring the same, but if you were you could model them as separate processes. It's like a competing process censoring model and that's a kind of survival analysis people do. Okay. Uh enough of that. Okay, let's build it up. Um we want a statistical model to do the estimation. Here's where the geometric distribution appears. Um we're going to ignore censoring for now. We haven't simulated it, so that's going to be fine. I'll say that the uh I assert the probability that a cat wait uh waits a cat I waits D days um and then is adopted. So, this is the probability of D sub I conditional on adoption equaling one. So, A sub I is going to be my censoring my censoring indicator. When it equals one, that means they were adopted. Yeah? Uh is equal to 1 minus P sub C where P sub C is the probability a cat with coat color C is adopted on any given day. Does that work with me? Yeah? Uh raised to the power of D sub I minus one. Why? Because that's the number of days they weren't adopted, right? So, each day they're not adopted is probability 1 minus P sub C. And it's so you've just got a long product of those. How many D sub I minus one of those of 1 minus P's cuz that's the number of days they weren't adopted. And then on the last day they're adopted and the probability of that is is P sub C. So, that's what's on the end. Does that make sense? Again, the axioms, they're your friends. Yeah? It's a joint event. We just multiply them. We just need what the event on each day, what's the probability of of what happened on each day. There were D sub I minus days where the event was was not adopted and then finally they were. So, the probability is 1 minus P to the D uh minus 1 times P. Good? Yeah? Um And uh let's do some prior predictive. Uh we've basically got the model, but we need a prior, right? So, let's do some prior predictives and I want a prior that does not have a bias in favor of any particular coat color. So, I'm going to assign the same prior distribution to all the coat colors. And I played around with this a little bit. I I'm not wedded to this prior, but I want something that is compatible with a huge range of probabilities of adoption. Some of these are really unrealistically fast uh as as we'll see. Um but, I chose a beta distribution um with alpha equal one and beta equals 10. That means it's a small value, right? So, this this is a beta distribution where the probability mass is towards the zero side. Yeah? Um it's like one success in 10 failures is what that means. If you're not experienced with beta distributions, one day you will have to be and you'll get used to it. Yeah, but that's what it means. The beta distribution is a probability distribution for probabilities, which is a fun thing to say, right? A probability distribution probabilities. But, sure, we can draw probabilities from probability distribution. Yeah? You can even draw probability distributions from probability distributions. Um There's an example of that in the book. Okay. So, good? Yeah? So, they're all over the place. This is a bunch of separate simulations. Um And these are just non-parametric survival curves. There's just the empirical decays of how many cats are still waiting to be adopted in that simulation at each point with separate simulations for black and non-black cats. Yeah, question? I have a question. Like doesn't this assume that um days of adoption are independent of each other if you write it Yeah. Yes, it does. The question was um doesn't this assume the probabilities of adoption are independent? It does transparently, I hope, assume that and that might not be the case. Yeah, absolutely. Um And there's also, I mean, there's a bunch of stuff in this data set that I have just not bothered to model. So, this is not If we were going I was going to really model this, I would gradually layer in other things. Like, for example, some of these cats are reclaimed uh went missing and are reclaimed by their original owners. And I'm coding that as adoption, but it's a different process, obviously, than like a stray, right, that is captured in an alley, and then people come in shopping for a new kitten. That's It's a completely different process, but that's a special code in this data set, reclaimed by owner. Uh these are Lots of these cats are chipped, right? So, uh I think it's the law in Germany you have to chip your your pets. Yeah, but it There are no laws in the states with respect to animals. No, there are I think like California you have to chip your cat as well. And um so, in those cases, they actually find the owners pretty rapidly. Uh then they come and it it's like a week or a couple days, reclaimed by owner, and I put those in here, and that's bad, right? Because this is a non-homogeneous process. Uh there's a bunch of stuff like that. And then non-independence, I'm not sure what you're thinking of, but I would gladly entertain hypotheses of such things. Um But just bear with me on the lack of realism so far. I want to teach you about the missingness and how to model it, and then you I think you folks know how to add in other complexities already, cuz you have the full power of generalized linear models at your disposal. If not, we'll give a talk at a conference someplace, and I'm sure our colleagues will Yes, absolutely. Okay. Good with the prior predictive? Um I encourage you to play around with this with the code, actually, and and see what you can get um uh different kinds of things by changing the beta. Okay. Um Yeah, also in the simulation, there are no cats that die. Like, they do have infinite lifespan. So, obviously, that's not realistic, either. Okay, the first survival model. Um uh We have synthetic data, so it let's test it on the test our survival model on the simulated cats. This survival model is just that, uh probability before. I'll show you some some uh code model code uh for it in a moment, but you folks are really familiar with Ulam code at this point. You just type that probability expression. That's the model. Right, what's the probability of D? I've already written it on the screen. You just type it as code. Uh So, this is just to show you that the model um is good at recovering uh the data generating values, which we know uh because we we chose them. And here is just on the left, I've got a case where black cats are adopted at a lower probability. I set it to 0.1 um than the other cats, where I set it to 0.15, and then we have posterior distributions. Uh and then on the right, I flipped it, right? Since I'm testing that that um it's not it's not prior sensitive. Uh which of course it won't be in this case, cuz the prior is symmetric, but still it's nice to You can always make a mistake. Trust me. I've made every mistake. Um Good? Yeah? This is just the, you know, responsible quality assurance workflow kind of thing I want to show you. Now we uh we've got a posterior distributions for the probabilities, but that's very hard to interpret. We want to push these out to predictions. So, again, I'm going to simulate using this using the same simcats1 function. I can just pass samples from the posterior back into that and simulate cats and plot their survival curves. This is very sample size dependent what the predictions look like, because if uh if you if you simulate a small population of cats and the probability difference in in the probability of adoption in in black and non-black cats isn't that big, there's going to be lots of stochasticity, which obscures it. So, that's what I'm trying to show you on the left. So, for simulated for cohorts of only 100 cats, there's lots of stochasticity in the decay, but you can kind of see that the black cats are pushing up up above the the orange trends, yeah? But there's a bunch of stochasticity. Uh Each of those trends is 100 cats waiting for them all to get adopted. Right? Um if on the right we do it for 1,000 cats, and now it's very easy to see that there's uh uh disadvantage. Now, this is simulated uh Oh, this is the real sample. Yeah, this is the real data. So, the black cats are adopted more slowly. It's not a huge effect, but after a long time, it's a big difference. Yeah? A small probability uh a small disadvantage on any given day after many weeks and months is a big is a big difference. Um This is This is quite sad. You can imagine those poor innocent black cats sitting there and watching the other cats leave. Poor things. Adopt a black cat. They're perfectly good cats. Uh I used to have a black cat. I loved that cat very much. Very good cat. Yeah. Very good cats. And yes, they're witches, but that just makes them better. All right. So, um Here's the censoring model. Let's put censoring into this. Uh All right. So, what we're going to assume, and this is very important, is that censoring is random with respect to the color of the cat. That's not true where we've got really bad problems. Um but uh uh you could imagine, for example, I don't know how this would even work mechanistically in this case, but that um black cats uh yeah, their data is mishandled or something like that. There would have to I can't come up with a mechanistic explanation. This is a credible assumption here. Um that that censoring itself uh the the fact of censoring is is random with respect to cat, but the color of the cat may influence the waiting times, and therefore the and the waiting time is what's causing the censoring. Yeah? So, it's not that the color of the cat doesn't I mean, the whole question is, does the color of the cat influence waiting times, but the the color of the cat is only indirectly causing the censoring. It's the waiting time itself that causes the censoring. You feel the difference there? You could see it in the DAG when you draw the DAG. I should have put the DAG on this slide. Okay, we're going to write simcats2 now. Simcats2 is the same function as before, but at the very bottom of it, I have added censoring. You'll see that there's a new argument in the in the head of the function called since, and I've set it to 30. That's the number of days I'm censoring it at. Any cat that is unadopted by day 30, that's when we stop observation. Right? Now, the it then the >> [laughter] >> Well, I did No, I mean, it's just the end of observation. That's Sorry, someone in the audience is is moving her hand across her neck in a execution motion. So, no, that is not what I intended. Um uh Yeah. >> [laughter] >> This is a happy day 30 event. I don't know. But uh Obviously, censoring in the real sample is not synchronized like this, right? It's like in the simulation, time zero is just when they come in, right? But in the real data set, the the dates are all staggered in the real data. I mean, they're not all censored at the same point like they are here. But So, this is not the the censoring model that's in the real data set, but it'll it'll suffice for us to develop a statistical solution, and it's easier to visualize. Is this good? Yeah? So, um here's what it looks like when you when you simulate from this. Um On the right there, we've got a cat uh index of cat across the bottom. I simulated 1,000 cats, and then the days uh they were there until they got adopted, and you'll notice that for this simulation, I set the censoring at 20. You'll see that big line of cats at 20. Those are all censored values. Okay? Um And so, again, we've got two strategies for dealing with this. This is I hope I've convinced you you cannot ignore this, but I'm going to show you what happens if you do. Yeah, but you really can't ignore this. Um Uh we've got two strategies. We can marginalize over these missing values. I'm going to develop that solution first. This is the most opaque one, but it's also extremely common. So, you're going to see it quite often. Uh And then second, we're going to do the brute force Bayesian approach, which I'm a huge fan of. Um Requires less thinking uh and more computation time. Uh And sometimes it's numerically more stable than the marginalization approach. Not in this particular case, but I've I've worked on survival analyses where the marginalization approach, even though I could write the correct functions, it was numerically unstable because of the the probability distributions involved. A gamma distribution was the criminal in that case. If you've ever worked with gammas, they have nasty tails. Really hard to do derivatives for parts of them. So, but then I just flipped it to imputation, and it it was great. Worked fine. Um So, it's nice to know how to do both of these. Marginalization. When we marginalize um in a survival analysis, we end up using something called the CCDF. The CCDF is the complementary cumulative distribution function. Um And I I talk more about what this is in the in the example in the book, uh and especially in the case study. Uh I'll I'll walk you up to understanding what it is. So, when we've observed an adoption, uh I've already uh defined that the probability of observing an adoption on day D is uh 1 - P to the D - 1 * P. Okay, that's what I call the data model, or you might call the likelihood. Yeah? Um Uh this is a a uh the probability of that event. It's a discrete probability. It's not a cumulative distribution. Yeah? It's the probability of the data you've observed. Um When we haven't observed something, we want to know the probability of getting up to that point. So, it's like what what's the probability that of all the cats not having been adopted that many days. Now we have cumulative uh probability distributions. But in this case, we want the complement of it, um because Think about this. What's the probability So, we observe a cat has not been adopted in his day D. Okay? Just think about that. Picture the cat. >> [laughter] >> Yeah? It's a cat. It's there. It's day D, it has not been adopted. You want to ask yourself, okay, what's the probability of this event? Well, what had to have happened is that um uh it was not adopted on every day up to this day. And you can write that, right? It's 1 - P uh to the D - 1. Yeah? Um and this is the complement of of the cumulative uh distribution, which would be the the dis the the a proportion of cats that um were are adopted by that day. And that'll be 1 minus this. Right? So, if you wanted to know how many cats will be adopted by day D, it's 1 minus all of the ones that aren't adopted, which is this expression that's on the bottom. So, this is the complement of the cumulative distribution cuz the cumulative distribution is 1 minus 1 minus P to the D minus 1. And this is the complement of it because it's 1 minus that. So, the 1 minus goes away. I know it's weird, but this is this is always weird. So, it's just so when you see this abbreviation CCDF, if you do Stan programming, you'll have to type that because you'll want the complementary cumulative distribution function. Um In this case, it's discrete, so we're going to write it uh literally. Uh but that's why you'll see people talk about the CCDF uh all the time. Uh but you can reason your way into it by just saying picture the cat, it's day D, it hasn't been adopted. What have to has to have happened on every day prior it must not have been adopted. And then you can write that directly. Um you don't have to go through the mental gymnastics of thinking about the cumulative distribution and the complement of it. Does that make some sense? Yeah? That's just I Sometimes I have to teach you weird notation even though I hate it uh because you're going to see it. Yeah, that we often deal with these issues and that's what the CCDF is. Okay, what happens when we ignore censoring? Um so, I get simulated data where we know the truth. This is how you test the impact of an assumption. Um so again, I've assumed that cats are uh black cats are adopted slower, uh probability 0.1 per day. Um uh this is this is order of magnitude larger than it is in the real data, by the way. Right? It's going to be much smaller in the real data. Um and uh non-black cats probability 0.15. Uh and the the solid densities there um uh the truth is in the um high probability region of the posterior. Yay. Bayesian updating works. Yeah, good times. Um uh and then the dashed ones are where I run the previous model uh ignoring censoring. Yeah? And you see what happens is it overestimates in both cases the probability of adoption. Let's think through just for a second why it overestimates. Because when we ignore the censored values, we're ignoring longer durations. And that leads us to think that the rate is faster. Cuz effectively, we've thrown away the cats that waited longest. Does that make some sense? Yeah? Uh there's different kinds of censoring. This is this is I think called right censoring. Um There's also left censoring and and and things like that, and they they may produce different biases. Uh but ignoring right censored data like this exaggerates rates. Um and you'll see it's a bigger effect for the black cats. Why? Because they wait longer. So, we throw away more of them. Yeah? So, it it not only biases all your estimates, it also biases the difference. Okay, have I convinced you? Yeah? That this is a So, this is a general problem. Um this is the kind of general problem you could actually prove analytically, and you can prove for simple survival models you can like write down expressions about what this bias is and stuff, which is cool. But you don't need to do that. You have the power of simulation at your disposal. You can always simulate out the consequences for these things. Um Okay. In the real data now, I fit these this model to the real data. And uh on the left, I've got ignoring censoring, and on the right we model it. I know these look the same. Uh so, I'm going to try to do a trick to help you compare them now. You just have to watch very closely. Are you watching? All right, here we go. Eh. All right. Did that work? Can I do it again? All right, again. Here we go. Yes. All right. So, I'm sorry. This is like the best I can could do. It's It's hard to visualize some of this stuff. Um So, yeah, it's the same effect. Uh when you model censoring, the durations are are it's higher up on the graph cuz they're longer. The waiting times become longer. When you ignore the censoring, it makes you think that the adoptions are faster cuz you've thrown out the records where the cats are not yet adopted. Make sense? Yeah? Um and again, in the beginning, the difference is very small, but by the time you're out to, you know, day 50, um uh it's it's a big difference in the proportion of cats remaining in both colors. Uh but the the difference between black and and non-black cats is is basically qualitatively the same in both models. Yeah, black cats are still at a disadvantage. Good? All right, I feel like I'm on time. Um Now, let's talk about the imputation model. We're going to do this again with a different computational strategy. Uh and and one of the things about learning this version of it is it will hopefully help you understand the Bayesian framework a little better in the way it provides solutions to measurement and missing data problems. Here's the idea. We're going to replace missing values with parameters, which seems like madness, uh but again, an unobserved value in Bayes is just that's all the parameter is. And sometimes those those parameters can't be observed because they're things like rates, which are not observable. Yeah? Maybe I should explain that a little more. Like you can't observe a rate. You can only measure things which let you calculate the rate. Rates are not observable, I assert. Yeah? There's a bunch of variables like that uh in in stats. But there are other uh parameters which are in principle observable, but we haven't observed them. Yeah? Things that we didn't measure, but we might be able to estimate from other things we measured. That happens all the time. Um and uh uh so, missing data is of that case where um we don't know it, but we replace it with a parameter. This is the same thing we did with measurement error last week. Right? We tried to measure something, and we partly failed. We've got an estimate of it. Uh the true value is a parameter. That is unknown, and we can get a distribution for that unknown true value. We can't know it for sure, but we can we can do better than the raw measurement. Yeah? And that's the idea. Um So, we we make a a vector of parameters, one for each missing value, one parameter for each missing value. And then we need a prior, right? Cuz you know you need a prior for every parameter if you know nothing else about Bayes. And uh Where does Where does the prior come from? You already have it. You have declared the generative model. And uh so, we're we're talking about something that is in principle observable, but hasn't been observed. The prior for it is the likelihood. I'll just say that again and hope it makes sense. The prior that you could have If you had observed this thing, the probability expression is the likelihood, which you have already written. You haven't observed it, but now if you say you were going to predict make a prediction, you'd use the same formula. You'd use the likelihood. And that's all That's the prior. The prior for Bayesians, the the the likelihood is a prior for the residuals. Yeah? That's all it is. Um and so, we don't need to do any further derivations. We've already got the probability expression. And that's what makes this so comforting as an approach uh and easy. Um Step three, we're going to sample 22,356 parameters. Okay. Good times. Uh Now, in previous decades, this would have been the end of the project. Uh but these days, this is just, you know, 20 minutes of your afternoon on a bad laptop. Yeah? So, it's not that big a deal anymore. Um and then profit, of course. Yes. Um Uh so, let's let me draw these ideas out a little bit more. See if I can finish on time. Uh so, for observed duration D, the probability of the data is just this um geometric distribution uh 1 minus P to the D minus 1 * P. Uh that's the likelihood or the data model. Um When it's missing, uh we have partial information about it. It's censored at a value D star, so we know that D is greater than D star. And we're going to use that. So, I I lied a little bit. When we define the prior, the prior is the same expression, but constrained to be greater than D. I'm going to show you in code how I do this in a second, okay? There's There's different strategies. I'm going to show you I think I've got an Ulam model and a Stan model implementation of this to show you the differences uh and how it works. I haven't showed you any model code today because you're pros, right? You're used to these ideas. But I'll show it to you in a second how it works. The script, of course, is on the website already. Yeah? For running all this whole all of the reproducing all the plots and everything. Um We know What do we know about the missing value? We know it's greater than D. Because on day D, the cat was still there. Yeah? It's like the counterfactual in the future, when will this cat be adopted? And you're you're going to make a you're going to bet on each day. That's the prior. The prior is your bet on each day. I don't usually like the betting framework for Bayes. There's a whole like philosophy of Bayes which is about betting, and I don't usually like that. But here is a metaphor I think it maybe it motivates the solution. Um That's what we know. And we also know that on each of those days greater than D, the probability is 1 minus P to the D minus 1 times P. It's exactly the same thing. Yeah. Does this make some sense? Yeah? Is this cool? Like, this These are routine problems, measurement error, missing data. Um this kind of direct solution is great because if you try to do this outside of Bayes, it it gets very difficult because you can't use inverse probability. You can't assign probability distributions to parameters and for if you're a frequentist. And that causes lots of problems. Um Bayesians are just like, "Yeah, trust the axioms. Let's see what happens." Right? We're very permissive. Uh the the Church of Bayes is very like, "We believe in the axioms and then we let a thousand flowers bloom." Right? It's that sort of thing. Of course, you've got to test and validate and be transparent about your assumptions, but we don't put philosophical constraints on what you're allowed to do with the axioms. Um Uh we're a big church. You can do frequentist calculations. No problem. Frequencies, those are things we calculate. Right? It's just It's perfectly fine. Um Anyway, does this make sense as a as a as a perspective? Um Okay. Here's the Ulam implementation that I came up with. Uh Uh two lines, the the first two lines in the formula list, there's some conditional notation I need to explain. I don't think I've used this before uh this year in this course, but this is in the book as well. Um when I do a couple of models in the second half of the book, where we've got um I'm using custom distributions which just inject code into the Stan. Uh this just lets me write any arbitrary likelihood that I like. And that's what I've done here. So, I'm saying um uh the probability um uh cat is when the cat is adopted, adopted equals one, is uh the log probability is what goes in custom. It's just the log of the geometric distribution for the days observed at that point. This is the same formula that I had on the previous slides. So, those are the observed cases. Right? So, this is just a way in Ulam to write an arbitrary likelihood. It just takes that code and dumps it straight in Stan. Um And then the second line is the censored cases, so adopted equals zero. Conditional on that, we use the same formula, but I've replaced days I, which remember is is is not been observed. Right? Cuz days till adoption is the thing that we want in this formula. We replace it with a parameter called days impute. And then the prior on days impute is just uniform with a lower bound days, which remember is a vector. So, for each cat, uh its prior for its imputed days starts on how long we observed it. And then I just put it an arbitrary upper limit of 2,000, which is big enough to deal with all the cats in the in the example. Um Uh and then I have the same prior on um uh on the colors as before. So, we've put the prior information we have, which is the likelihood, cuz that's the prior, that's the second line of this model formula, and then the third line puts the constraint that it be greater than the observed days so far. So, that makes sense? And this will do it. Um it's a lot of It's a lot of cats though. That's That's 22,000 uh parameters, and we declare them even for the ones that were observed to be adopted, but you know, the computer can handle it. It's It's no problem. And um uh this runs just fine. I'll show you the results in a moment. I think it might help to show you what the Stan code looks like. Um This is not the Stan code produced by Ulam. This is my manicuring Stan code version of it that's a little bit better because I get to inject the observed values back into the parameter vector. Um So, what's going on here? If you haven't worked with raw Stan code before, uh it's very much like Ulam code. It's just all Ulam does is automate these top two blocks, or sometimes more top blocks if you have them. Uh uh The top block declares the observed variables. We usually call those data. Yeah? Uh where N is the number of of cats uh in the data set, and then whether or not they're adopted, and then the days, and then their colors, and then the parameters are the unobserved variables. We have two We have a vector of length two for the probability of adoption P. And then we have a vector of length N, which is the number of cats, called days imputed, and the lower bound is days. So, this is just in raw Stan notation, the way to declare a uniform uh prior probability on this. But this has no upper bound except infinity. Right? Which is a little bit more comforting. Yeah. Uh that doesn't occur in the data set. Right? And we're not going to sample that either. Um And then the model block uh uses the likelihood again. You'll see here, we've got this days impute vector. And then um we use the same likelihood down here uh in the case of non-adoption. You see there's this if-then branch. Let me remove the boxes for a second. In the model block, if adopted I equals one, then we just do a target update using the geometric distribution. And then there's this trick where I copy the observed days into days imputed. Right? This is optional. It's just a trick cuz I'm going to show you in a graph later um which ones are observed or not. You don't have to do that. Um and then uh else uh adopted I equals zero, we use the same expression, but we replace days with days imputed, just like in the Ulam model. And that's it. That's Bayesian imputation. Very straightforward. Yeah. Connor. >> My My naive assumption would be that the parameter for days imputed would just be like days plus one would be the most probable. It's the most probable, but we want the probability of every Okay. Right? And it's There's a geometric distribution over the days above day D. Oh. Right? So, it does decay into the future, which is why infinity doesn't isn't relevant. Uh no cat will wait infinity days. Well, if you have infinite cats, uh at least one of them will wait infinity days. Yeah. >> [laughter] >> With probability one. But um but we don't need to worry about that problem. Uh But you're right. Your intuition is right, but we want the distribution. We don't want the the the most likely one. Yeah? Okay? And we're going to get posterior distributions. I'm not going to show them to you in a second. Okay? Well, I'm going to I'm I'm going to show them to you as intervals because there's a bunch of them. All right. So, so uh yes, um uh data visualization is my passion. >> [laughter] >> This is really excellent. So, this is the when you plot tracey for this model, this is what you get. Isn't that great? Okay, so every row here is uh is a days imputed and for a different cat, and there's uh yeah, uh a lot of them. Uh 22,000 of them, to be exact. Let me zoom in a bit. All right. Well, that just makes it worse, doesn't it? It's just really terrible. All right, what's going on here? Um So, across the horizontal is days. This is in unit days. This is synthetic data where I've made censoring always happen at day 20. So, and the points uh left of day 20 are the observed values, which I just copied into the vector of parameters. There was that line where I just said it's normal observed days with a standard deviation of 0.01. That just fixes the parameter estimate at the observed value. Just a trick so I can make this graph. Right? You don't have to do that. Um So, on the left of 20, those are the observed values that were simulated. And on the right are the posterior distributions for the days uh how long the cats would wait that had not yet who were censored. Right? And there's really only two posterior distributions here because there's one for each color type. There's sampling variation, so they look at their intervals look a little different, but you'll see that the posterior means at the points, there's really two categories cuz there's only two There's only two likelihoods. Right? If you stratified the cats by more features of the cats, you'd get more variation in these posteriors, which you're welcome to do. In the raw data, there's lots of things they know about these cats. There's age estimates, which almost certainly has an effect on this. Right? Uh other things uh that go on. Does this make sense? Yeah? So, this I think addresses your question a lot, Connor. It's Yes, on the far left, it's maximized there, but then it decays geometrically, literally, geometrically, and the mean is pretty far from the maximum. So, if you're going to say, "What's the expectation?" you'd be a fool to bet on D plus one. Does that make sense? Again, the gambling metaphor is not usually my preference, but maybe in this case it makes sense. It seems a bit macabre to bet on how long a cat will wait to adoption. This is like That's like prediction markets for that seem terrible. Yeah. Um Good? All right. Um All right, I'm kind of on time. This is amazing. Okay, so let me let me uh try to summarize this. Why impute? Um Marginalizing out missing values is is faster if you can do it, but it's opaque in the coding, and you have to do some extra derivation, and sometimes that derivation is difficult. Right? It Or the algorithm that's needed to do it can be quite opaque and difficult to get to work. Uh there's a family of extremely useful models called hidden Markov models, which are a kind of missing data model. Uh as in a hidden Markov model, the states of the system that you'd like to know are are not measured. They may not be observable, but sometimes they are observable and you just don't have them. And the data you have are called emissions from those states. But different states can produce similar emissions, but with different probabilities. Um used a lot in animal movement models. Like, where's the shark? Right? And you're getting periodic satellite pings of the radio-tracked shark. So, you don't always know where it is, and there's measurement error. Uh so, it's a hidden Markov thing, and you're trying to estimate not just the position of the shark, but what it's doing. Is it feeding? Is it resting? Things like that. Those are hidden states. So, hidden Markov models used a lot um not just with sharks. Uh I have a a colleague and co-author on the Bayesian workflow book. This is what she does. She does hidden Markov models of sharks. Um uh Vianney on the cover. And um uh it's a really cool there's a case study in the in the Bayesian workflow book for shark movement. Uh when it comes out you'll like it. But uh in those models there's an algorithm for doing the marginalization and it's a mess. Like it works. Like there's nothing wrong with it, but you don't want to write it yourself. Luckily, there's a Stan function that does it for you. Uh but that's that's a kind of obstacle if you have to develop something like that for yourself, it's non-trivial uh thing to do. Um and sometimes it uh marginalization is numerically unstable. As I said, I had a project where there were gamma distributed survivals and I just couldn't get it to be numerically stable. I flipped it to doing imputation and it ran no problem. Right? It seems crazy I'm going to invent, you know, tens of thousands of parameters and things will be easier. But, you know, the universe is weird like that. Yeah, so it's nice to know. Um imputation is slow, but it's very transparent. It's easier to read the code and understand what's going on once you get past the philosophical weirdness of of declaring parameters for missing values. Um it also gives you estimates for the missing values, which in marginalization doesn't. And sometimes that's the goal. Like you really want to know those estimates. In this case, it's not, but sometimes that really is the goal. Um ecologists uh use something called occupancy models quite a lot. You want to know if an animal is there in a patch. Um but the problem is that there's a detection probability. Many animals are cryptic. They don't want you to see them. So zeros are a mixture of absence and measurement error. That's called the occupancy model. It's the most basic statistical problem in ecology and it's a mixture model. It's already hard. Um sometimes the occupancy is the thing you want to know. Right? But you can't observe it and you need to posterior distribution for it. So imputing it is essential in a model like that. Does that make sense? Yeah. So for fertility models, no, you want to censor a woman at 20 years old? Yes, absolutely. Absolutely. Absolutely. Great example. Yeah, yeah, there's an example. Sorry, the microphone won't pick you up I don't think so. I'll I'll repeat it. Yeah, an example of fertility modeling if you have a woman that's only experienced part of her fertility. Yeah. Right? And you want to know her total fertility in her lifetime. This is a missing value. And and you need to impute that to do the regression for lifetime effects and things like that. Yeah, you've worked on a bunch of these problems, I know. Yeah, yeah. It's a great example. Um Sorry, I've gone over 2 minutes somehow in my enthusiasm. I am sorry. Anyway, I'm going to stop here. Um uh there's more examples of missing value stuff in the chapter in the book and in the uh previous year's lectures as well. And if you ever have a problem like this, you're welcome to come to me and ask for help on it. I'm very passionate about missing data. Thanks for your attention today.