speaker-0: My guest today is Matis Hollanders, a postdoc at the University of Newcastle in Australia. Yes, there is a Newcastle in Australia, where he spends his days analysing camera traps and song meter datasets and building Bayesian models that take that messy ecological data seriously. The reason I wanted Matis on the show is his R package, Ocaru. It does something I find. Conceptually clean and refreshingly honest. Instead of collapsing camera trap data down to did we see the species today, yes or no, it models the counts directly. That choice unlocks a whole family of richer questions about detection rates, multi-species co-occurrence, and how all that varies in space and time, and the package handles it all with, you guessed it, hierarchical Gaussian processes. Global, local, shrinkage priors, and some of the newest stand features around zerosum and simplex constraints. So, as you know now, all of these are a bunch of show favorites. So if you have ever wondered how Bayesian thinking actually helps ecologists get more out of their data, this is the episode for you. This is Learning Bayesian Statistics, episode 159, recorded May 13. Let me show you how to be a good basier change your predictions after taking information in. And if you think it now be less than amazing, let's adjust those expectations. What's a Bayesian is someone who cares about evidence. Welcome to Learning Bayesian Statistics, a podcast about Bayesian inference, the methods, the projects, and the people who make it possible. I'm your host, Alexandora. You can follow me on Twitter at Alex underscore andora, like the country, for any info about the show. Learnbase stats.com is Lab Plus2B. Show notes, becoming a corporate sponsor, unlocking Vision Merch, supporting the show on Patreon. Everything is in there. That's learnbasetats.com. If you're interested in one-on-one mentorship, online courses, or statistical consulting, feel free to reach out and book a call at topmate.io/slash alex underscore andora. See you around. Folks and best Asian wishes to you all. Mathies Hollanders, welcome to Learning Basion Statistics.

speaker-1: Thank you.

speaker-0: Yeah, take it. How is my Dutch pronunciation?

speaker-1: It was perfect. You nailed it.

speaker-0: Well, we'll see if at some point I learn I learn Dutch. I think I will not until like unless I have to go leave there at some point, to be honest.

speaker-1: It's hard to justify learning another language, especially if you have no reason to speak it on a regular basis.

speaker-0: Yeah, yeah, exactly. ⁓ especially since ⁓ Dutch people actually speak so like such a good English, you know. So d all the Dutch people I've met in my life speak perfect English. So I'm like, Well, you know no incentives.

speaker-1: It's pretty well the last country you'd actually need to learn the language to be able to communicate with the people there. So

speaker-0: Yeah, yeah. Yeah yeah. I have ⁓ I have a friend who lives there right now. He's German, so actually it's very close, like German and Dutch, or super close, but he still hasn't learned Dutch because everybody speaks English.

speaker-1: Well to to be fair, I I don't speak a word of German. I I don't reckon the two languages are all that close together, but ⁓ maybe it just shows some of my skills of learning languages. I mean

speaker-0: I you live in Australia. So you know. Yes. You're not l you're not living in Germany, so the like the the comparison here is not ⁓ you know it's not a a perfectly ⁓ controlled trial. but actually we were talking about that just before recording, so you're the first you're the first person I meet in I in in I'm the first person you meet who has whose last name is also the name of a country because your name means From from Holland, basically, from what you told me, right?

speaker-1: That's right, that's right. Yeah, I think it's it's fairly unusual.

speaker-0: Okay, well, I think it fits well because ⁓ we've got the World Cup going on. and so so perfect. Today today game is ⁓ is Andorra Andorra Netherlands. So you know. That's perfect. ⁓ that is perfect. And actually let's ⁓ let's dig into into your background, ⁓ because you do you do interesting stuff. So maybe tell us what you're doing nowadays and How you ended up working on wildlife data analysis?

speaker-1: Right. So ⁓ I'm currently doing a postdoc at the University of Newcastle in Australia, where I'm looking at ⁓ analyzing basically camera trap data sets. before that I was working on song meter data sets, which ⁓ it's microphones hanging up in the bush that record information continuously, basically. And there's a lot of overlap between those two data types. So that's what I've been focusing on the last couple of years, I'd say. ⁓ before that I was doing my PhD and ⁓ I was doing mark recapture analysis, which is basically just another type of hidden Markov model that I'm doing right now. So I think a lot of a lot of the modeling that I've done has kind of been ⁓ very similar but answering different questions.

speaker-0: Okay, yes, so I already see ⁓ some some similarities with other episodes we've recorded. The one with ⁓ Aaron McNeil obviously comes into mind, ⁓ where ⁓ actually also worked on on that model with him. He's leading a fantastic team and they are they are doing modeling to infer the trade of shock meat. around the world and and so they are using a bunch of very very cool Bayesian ⁓ methods to ⁓ to do that. It's an extremely complex model, one of one of the most complex I've I've ever worked on. ⁓ but that that has been super fun to work with them on that. And that was episode one twenty seven. I'll put that in the ⁓ in the related episodes for this one. ⁓ and also the one with ⁓ Vianellios Baraja. She is doing a lot of hidden Markov models. So ⁓ this one will ⁓ will hit close to home for you, Matt. ⁓ and and so people who wanna wanna hear more about these kind of models, I will put that in the in the show notes. and that was episode 14 with Yanni. So this is ⁓ this is one of the one of the oldest ones, damn. ⁓ so Lots of ⁓ hidden Markov models, things like that. So fun fact for ⁓ listeners, I actually so piz. I I saw on your email signature that you were at the University of Newcastle. and ⁓ as an ignorant ⁓ person who just ⁓ relies on their priors, I thought it was Newcastle, England. And so and so that that was That was interesting ⁓ when you told me you were in Australia. I was like, ⁓ wait.

speaker-1: I have to beg you to change the time, otherwise I would have been here out in the middle of the night.

speaker-0: Yeah, yeah, the original time was one AM for you. So, you know. but actually it's better. Now that I'm in California, it's actually much easier to record with ⁓ Australians than it is to record with ⁓ with people in Europe. So I I thank you for that because otherwise now I have to wake up very early to record the episodes. so so that's good. And actually how did you end up in in Australia? ⁓ is that related to ⁓ to the the kind of models and data and And animals you have access to or that was just kind of random.

speaker-1: It was definitely related to the animals mostly. ⁓ I've been interested in wildlife my entire life and ⁓ I'm actually from the Netherlands, but I spent some time in the US, which really opened my eyes to what it's like to live in a place with lots of wildlife and I always knew I didn't want to stay in the Netherlands and I'd shortlisted Australia as a place I'd like to live and during my masters I just looked for PhD opportunities ⁓ in Australia just to get me over here. because I just wanted to live here to to live here. The the specific choice of PhD was a bit of an afterthought. ⁓ but yeah, that's that's how I ended up here. And ⁓ yeah, so my work didn't really bring me here, but during my PhD I became very interested in the work I'm doing now.

speaker-0: Yeah, yeah, yeah. Okay. So it's no, it's really you ⁓ you really fell in love with Australia in a way.

speaker-1: That's right, yeah. And it was I think to be quite honest, it wasn't the best to do my PhD for those reasons. I don't think I did a PhD for the right reasons and I found it quite stressful at the time. ⁓ but here we are, I landed on my feet, so ⁓ can't complain.

speaker-0: Yeah, yeah, yeah. It sounds like it. and yeah, I still have to to visit Australia. I visited ⁓ New Zealand a few a couple years ago, actually with Aaron, ⁓ because ⁓ we were working on the model over there and that was absolutely incredible, and I've been meaning to come back since then and also come through ⁓ through Austral Australia at that time. So yeah, definitely wanna do that at some point.

speaker-1: Well, maybe before you come to Australia I can convince you to come to the Yucatan Peninsula in Mexico for the International Statistical Ecology Conference in January twenty twenty seven.

speaker-0: ⁓ yeah, you definitely can. That's that sounds very fun. ⁓ and actually much closer to home. So you know, ⁓ for sure. Let me ⁓ like ⁓ yeah, send that my way and ⁓ and I'll see ⁓ how I can make that happen. That'd be fun to to record some ⁓ some episodes over there actually. That would be ⁓

speaker-1: Fantastic conference. It's every two years. ⁓ four years ago I think it was in South Africa, so I could really justify going. ⁓ the time after that it wasn't Wales, so I couldn't really. But now that it's in Mexico again, ⁓ where we have neotropical rattlesnakes and lots of margaritas, I think it'll be it'll be easy to go.

speaker-0: Yeah. Yeah yeah that sounds ⁓ that sounds perfect, yeah. Yeah, send send the links my way for sure, please. I will. And ⁓ and what's about Bayes, actually, because you do a lot of Bayesian models, but what drew you to Bayesian stats specifically?

speaker-1: Well, that's that's also easy to answer. ⁓ it kind of starts with being at uni. ⁓ during my master's and bachelor's, I think I took three or four ⁓ classical statistics classes and I I never I never understood what was going on. I always found it incredibly not intuitive. And therefore when I started doing my PhD and I knew I was gonna have to analyze these ⁓ kind of complex mark recapture data sets. I needed to have a proper analysis. So I was like, okay, I need to actually sit down and learn this stuff. And it was the introduction to Bayesian statistics in win bugs, ⁓ by Mark Carey, one of the like introductory jags, bugs, books for doing analysis. And it ⁓ it's When I feel like I need to do something, I can either rush myself and get frustrated or I sit down and just go as slow as I need to go until I get it. And that's what I did with this book. I went chapter by chapter, like hand copied all the code. It never ran, always took it twenty times to rewrite it, even though it was only like ten lines of code, ⁓ because I literally had no idea what I was doing. But as I went, I became very interested in it and ⁓ Yeah, that that was definitely the the beginning where it all started probably six years ago now. And since then it's completely changed the way I think about everything because because I've always written stuff in jags or bugs, like you have to write the generative model out, which I'm sure I don't need to tell you how ⁓ how good that is for forming an intuition around what you're doing or a getting an actual mechanistic understanding of what a model is doing. ⁓ so that's definitely where it started. I have I have to credit Mark Carey really for ⁓ all of those books that that they wrote for ecological modelling.

speaker-0: Wow. Yeah. Yeah that ⁓ that makes sense and and yeah, I these models are very complicated. So starting with those, this is ⁓ yeah, that's hardcore. Well done because ⁓ this is clearly not the easiest.

speaker-1: It's they are both complicated and not complicated in the sense that ⁓ they're they're complic they're only complicated or not only, but one of the main reasons they're complicated is because the data is messy, because we're not perfectly observing animals, because we're not surveying at set times all the time. So we have all this like almost like infrastructure around the model that is complicated that you need to account for. But It's also not complicated in the way that a lot of these complicated models can be thought of as just some flavor of GLM still. ⁓ which I'm sure we'll get into it in a bit. ⁓ but yeah, very I I'd say most of the complexity comes just from how messy the data is that you need to accommodate. But under the hood, you're kind of still doing classical things like GLMs and linear models on different model components.

speaker-0: Yeah. Yeah yeah. ⁓ and they are they're super fun also to to work on, to be honest. and and that's also because they are not that easy. ⁓ that you always have something to learn. There's always a complexity or nuance you didn't know about, and so there is always something you can you can add on top of it and that that really makes a difference. ⁓ and actually if we so let's set the the stage for the for the listeners who don't work in ecology. ⁓ again I recommend listening to the two Irons and Vienna's episodes because we go into the into the the details of capture recapture models, especially with VNA. ⁓ but you focus a lot and a at at least today we'll focus a lot on occupancy models. So yeah, set the stage for the listeners about what that is and what problem that's actually solving.

speaker-1: So occupancy models are ⁓ are used to accommodate for the fact that we don't always detect a species when we go and survey for it. So a lot of animals are rare, so you might only see it once every few times that you go and visit a site. So instead of just modeling the naive detections, like let's say we have a hundred sites and we go and visit them once and we find the species at fifty percent of sites. We would never think that the species only occurs at 50% of sites. We only detected them at 50% of sites. So what occupancy models do is they add this repeated measures component where you survey every site multiple times, and that lets you estimate a detection probability. Where, for instance, a common, highly detectable species, you see it almost every time you go and survey at a given site. So if you're never getting it at a site, it's probably not there. ⁓ whereas if there's a rare species where the sites where you did get it, you only get it one out of every three visits, every five visits, every ten visits. So that's a rare or hard to detect species. And therefore some of the sites where you never detected the species, they might be there. You just didn't pick them up because basically what it's telling you is that if the species is there, You're still unlikely to see it.

speaker-0: Hm. Yeah. Yeah, that's a very good ⁓ I think that's a very good way of explaining it and that makes it very intuitive because yeah, like everybody can understand that some species, some animals are harder than than others to spot. And so basically

speaker-1: Another a a different way that I like to think about it now is that an occupancy model has a zero inflation component, which happens at the level of the site. So the feature of occupancy models is that every site is multi is surveyed multiple times. And what we're estimating, the first part of this model is is this site occupied yes or no? Which is a latent variable because we we know it's occupied when we found it, but we don't know ⁓ it's not occupied when we didn't find the species. So we're estimating this zero inflation parameter at the level of the sites, but conditioned on that ⁓ occupancy state, we have this detection history, which lets us estimate the detection probabilities or detection rates at those sites. And that's kind of what I meant before about ecological models are complicated, but they're also not. If you can sort of wrap your head around that the zero inflation is not happening at the level of the observation as it does in a typical ⁓ zero inflated model, but it happens at the level of the site, then you're kind of back at just this, you know, model that's in our our toolkit.

speaker-0: Yeah yeah yeah. and added to that you also work out on automated recording units. So what which is ARU and so what's what's that? I surely don't know anything about that. What what is special about that kind of data and why do traditional occupancy models struggle with it if they do?

speaker-1: So automated recording units are things like camera traps and ⁓ microphones that are hanging up in the forest basically to detect birds or vocalizing mammals and it's in the name. They're just these devices that are essentially monitoring continuously over a deployment period. So a camera battery or an S D card might last like two weeks, but d during those two weeks it's running the whole time. But there's ways to extend it as well. Like you could have battery packs, or you could only have the camera be on for one hour every day, could be on every other day. ⁓ so the key feature of it is that it's this thing that's actually monitoring continuously in time versus the traditional types of surveys that we would do where a human would go and visit a site for half an hour, six hours, two days. And look for the species and then record in our notebook if we detected the species or how many we detected. The key feature is that the surveys are not discrete, really. They're occurring in continuous time, ⁓ which is different from how the traditional ⁓ survey methods work. And I wouldn't say that occupancy models struggle with the kind of data. It's just that the models were. kind of designed for discrete surveys as opposed to stuff that's going on in continuous time like these kind of d devices.

speaker-0: Mm-hmm. Okay. Yeah, thanks. That's very clear and that actually illuminates some of the modeling choices we're gonna we're gonna talk about next because ⁓ you created an R package to just that kind of models that you called Ocaru. ⁓ so occupancy and ARU and actually is an Ocaru something special in your field? Is is that even an animal?

speaker-1: I don't think so. The way the way I thought of pronouncing it was AUK A R U, but maybe this is a good a good signal that the name isn't chosen that well.

speaker-0: Yeah, I mean I like okaroo. It it it sounds ⁓ it sounds like an animal, you know. Look at that okaroo. Yeah.

speaker-1: Yeah, you're you're welcome you're welcome to call it you're welcome to call it that.

speaker-0: So actually the most important part is what is it about? You know, ⁓ why did you create that package? ⁓ and what's the ⁓ elevator pitch for it?

speaker-1: Well, elevator pitch might be hard because I feel like I always end up talking a long time when someone asks me this question, but I'll I'll try and keep it brief. ⁓ there's kind of two motivations here. The first one is that in a traditional occupancy model, we normally record whether or not we detected the species, yes or no, during the survey. So we collapse our data to a detection or non detection event, aka ⁓ one or a zero. ⁓ and a lot of analyses that are done even with ARU, so camera traps, people will end up, for instance, binning their survey periods to one day or one week of camera trapping, and then collapse any detections of the target species to a one or a zero, as opposed to modeling the counts directly. And my view is that this kind of happened because, well, when occupancy models were proposed. the what they called the observation model. So the detection side of things was a Bernoulli or logistic regression where you modeled the detection or non-detection during a given survey. In my view, this collapses the actual raw data, which is a count, to a one and a zero or a one or a zero before doing any modeling at all. And you are throwing away information by definition when you collapse a count to a one or a zero. And the my main motivation here was to be like, we shouldn't do that. Or if you're gonna do it, you need a justification as opposed to doing it by default. So there are issues, of course, like ⁓ autocorrelation between detections, like you might have a group of 20 kangaroos that's standing in front of a camera for two hours, which will really ramp up your detections. So I put some features into the package which lets you thin detections. So for instance, you could have a 30 minute thinning window where it retains the image with the highest count during that window. So going back to the kangaroos, let's say they are there for two hours. ⁓ that would be reduced to four detections with the maximum number of kangaroos in an image during each 30 minute window. So it might be like the group ⁓ in the camera might range from like seven to seventeen individuals. So it would retain the highest count, because that's the highest that we know was there at any given time during that window. But regardless of the thinning, the idea is to just not just go to ones and zeros, but model how many did we see because all else being equal, ⁓ if you're seeing more kangaroos in one month than the next, there were Probably more kangaroos in the system in that time. And I'm treading carefully here because this is a highly contested issue in ecological modeling where the holy grail of ecological modeling is estimating abundance, ⁓ for which the gold standard is marked recapture analysis, where you have actually tagged individuals, recapture them over time. And you're essentially modeling these time series of individuals surviving in the wild. But the key the key part of those models is that we have individual identification, which means we can track them over time. There are methods that exist where you use the counts directly to estimate an abundance, but those models are quite sensitive to model assumptions. that are very hard to meet in practice. And one of them is, for instance, that you need to have independent count or counts of independent individuals. If you're double counting individuals, that throws out the abundance estimates. And that gets very hard to do with ⁓ animals that are hard to individually identify, like kangaroos or, you know, plug in any any mammal that doesn't have clear patterns like stripes or spots or whatever. ⁓ and it gets very hard to do. But nevertheless, I still don't think that that justifies collapsing counts to just ones and zeros, because how you choose to interpret the information that's in those counts is up to you. But you might as well model them before just collapsing it to a one and a zero before you even start your analysis.

speaker-0: Yeah, yeah, yeah, yeah. Yeah. And so actually that was one of my following up questions was about that choice that you made in in Ocaru, which is to model the counts directly, ⁓ rather than treating the detection as a nuisance parameter, which is exactly what you ⁓ which you just talked about. ⁓ yeah, I can see why why that's definitely a meaningful shift conceptually in how that'd be very useful. And yeah, basically instead of throwing away that information, you you do want to have it as much as possible in the model. And then if you actually want to throw it away afterwards, well maybe but but you don't have to. And

speaker-1: Exactly. You you can just you can just fit the count model and then be like, I'm not interested in in in ⁓ interpreting this at all. And that's another thing actually where in occupancy models, usually the focus is on the occupancy of any given site, which is just to say, does this species occur here at this site? Yes or no. That's all occupancy is. It's just a one or a zero for if the model thinks that the species was there. My another argument that I'm putting forth is that when you have these continuous deployments of camera traps, simply focusing on this binary presence absence of species is almost kind of boring. Like we have much richer data, right? Like we potentially have these cameras running for months at a time, ⁓ potentially repeated over many seasons, and to just basically have your target of inference be if the species or at which site the species occurs, ⁓ that still comes out of the package because that's just what occupancy is. But I think it's much more interesting to focus on variation in the detection rates ⁓ over time, especially when you get longer deployments. Like I've had some data sets that are many years and you get to see these crazy fluctuations in underlying detection rates of all your species. And it just gives you a lot of insight about what's in your data that's beyond just occupancy. So what I'm saying is this package still gives you occupancy like any traditional package might do, but it just also lets you focus on those detection rates. So you don't lose anything, you only gain stuff. And it all comes from the fact of just basically not throwing away information in your data before you even begin modeling.

speaker-0: Yes. Completely agree. ⁓ and you have also some very cool stuff in there where you have the options for hierarchical, ⁓ special spatial and temporal Gaussian processes to help for multi species modeling. Obviously I I was very ⁓ very charmed by that option because I'm ⁓ I'm a fan of GPs, so ⁓ can you can you unpack that for us? Can you Does what that's doing and what scientific questions it lets ecologists ask that they could not before.

speaker-1: Yeah, sure. I'll I'll start with just the justification for having these spatial and survey effects in there to begin with. And just from first principles, when we essentially have repeated measures at two levels, at the level of the site, because every site is surveyed over potentially many days, many weeks, many months. And also at the survey, because usually we have a grid of cameras, so every survey period, let's say every day or every week that these camera traps are active, we're actually getting measurements from multiple camera traps. So we have repeated measures at the level of the site and survey. And there's plenty of literature that talks about what or how you should handle that. But the the key thing is that you ⁓ essentially have pseudo replication in your data if you don't account for that somehow. So the way I account for that is a random site effect and a random survey effect. And then it's just a question of how do we actually model these random effects. So typically you might have a normal distribution ⁓ where all of these random effects are coming from, where you're estimating some standard deviation. That's just your run-of-the-mill random effect as implemented in any R package. With our site and survey effects. Those sites and surveys aren't necessarily independent because you would think that sites that are closer together are going to be more similar, and you'd also expect surveys that are close together in time would also be more similar to each other. And once you get into this spatial and temporal sort of correlation structure, you're basically knocking on the door of a Gaussian process where you're constructing a covariance matrix over all the sites and over all the surveys. So The way I try and explain this as an elevator pitch is it's just a random effect that looks at these spatial and temporal locations of where we're collecting our data. Where the multi-species comes in is the package is designed to not just look at w data from one species, but to look at data from all the species that are recorded. And the motivation here is when you look at your raw data. Of a camera trap or a song meter, essentially what you have is a timestamp and a species label of what we found. And that means if we have multiple species, we probably should have a random effect for every species and every site because not every site's going to be good for each species. ⁓ And then you get into the multi-species Gaussian process thing where you still have the same Gaussian process kernel. But now you have a realization from this for each species. And there's actually a way to do this to make it quite computationally efficient that you only have to do one Kaleski decomposition, even though you're getting these draws for all of these species. And it ends up being matrix normal, which I hadn't come across matrix normal until about two years ago. But it's basically a matrix of values where you have a row and a column covariance. So if you have a row covariance that governs the correlation between your sites, you have a column covariance that governs the correlation between your species. So what you get is not just these spatially structured random site effects and random survey effects, but you can also directly estimate how species are correlated in these effects. So for instance you can estimate whether two species are likely to co occur in space and time based on this multi species process. ⁓ now I imagine that's quite a lot, but that's that that's what it gets you that I don't think I've seen implemented in other packages. And if I had to if I had to say how the package like clearly differed from other things that are available, it would probably be ⁓ how the multi species hierarchical nature of the model is handled.

speaker-0: Mm-hmm. Yeah. Yeah. And that definitely makes sense. I think it's a it's a great setup for a g for a a GP. So that that makes sense that you have that in there. And and also then the hierarchical nature of your data also marries really well. So that that's super fun to see that in there. How how did you end up implementing that? Like are you using vanilla GPs, are you using approximations? How is the how are you handling the the speed and efficiency side of the GP sampling?

speaker-1: Yeah, there there's a couple of ways. At the moment, the GPs are still exact, and that's simply because so far most of the data sets that I've come across haven't been crazy in the number of sites that they have. So often you'll have like hundred or a few hundred, which I feel still works well for an exact GP. ⁓ also I have built in a feature where Often we'll have different grids of cameras in different locations. So essentially independent grids. And when I have independent grids, it does a separate smaller GP for each region, which just assumes that the regions are independent. They still share the same hyperparameters, but you don't need to construct a full covariance matrix. You do it by region. So if you have a thousand cameras in your data set, That would be an expensive GP to compute. But if those thousand cameras are split between ten different grids, then you've got ten smaller ones, which is actually way easier to manage. And it is definitely ⁓ something I want to add is some form of ⁓ approximation. And someone recently told me about GP tools, which is a a I think it's some sort of package that's implemented in Stan that lets you base I think I want to say it's Fourier transform GPs, ⁓ but some sort of approximation for the inevitable case that you do have a huge data set where you actually need a massive covariance matrix.

speaker-0: Yeah. Yeah, I don't know GP tools per se, but maybe they are using HSGP and Ratherhood most of the time. ⁓ if you're in a maximum 3D case for GP, that's gonna be that's gonna be extremely helpful. and so for spatial temporal data should work because that's three D. ⁓ and that's yeah, the same flavor as if we transform. And that really is ⁓ really is a game changer for big data set because it just becomes almost linear in in compute time with ⁓ with your data data size. So that's

speaker-1: Yeah, no, I I I do have code laying around 'cause I've implemented them in the past for the HSGP stuff. So it's it's honestly it's probably something that will be in the package in the next month or so.

speaker-0: Mm-hmm. Yeah. Yeah. I mean, that's really ⁓ that's really helpful. I mean if most of your data set are not too big, that's awesome. You know, that means you can't use exact GPs. ⁓ this is rare, but ⁓ this is this is great to hear that is that it's possible. ⁓

speaker-1: It's not that rare in camera trapping because just most most of the time people are bottlenecked by by resources and ⁓ I think it would be not uncommon to have a data set of thousands of cameras. I do think it would be unusual for all of them to be part of the same grid as opposed to ⁓ split into multiple smaller grids where you can actually just do them separately within each region. But even there the approximations would ⁓ end up being faster, I'm sure. So ⁓ worth implementing regardless.

speaker-0: I mean yeah, I mean it depends on the approximation because you have some approximations that actually work less well than a vanilla GP if the dataset's not big enough. So, you know, kind of like ⁓ central limit theorem kind of stuff where it's if you're not in the regime where the approximation ⁓ accuracy kicks in, then you actually want the vanilla GP. So use the vanilla GP, just do that.

speaker-1: Okay.

speaker-0: It's like you you have the mathematical guarantee and and that's great. In my case that's been super rare. Like most of the time I cannot use vanilla GPs and I've had to use approximations and thankfully there are some now. And also the compute we can use is just so much better now. ⁓ but yeah, like it I mean it's good to have the different options, I guess.

speaker-1: I have wondered about that with a spatial GP where you have two inputs, a well basically a latitude and a longitude. And yeah, the way I understood the HS GPs, you have a basis function for each dimension, or a number of basis functions for each dimension. So if you set it to the default of like twenty, it ends up being twenty times twenty, I think, because you need the combinatorics. So you end up having to estimate a lot of coefficients, which When you have a a smaller field you could I could see how that would be ⁓ not ideal.

speaker-0: Yes, yes, exactly. ⁓ actually the the choice of hyperparameters for HSGPs is extremely important and and is gonna impact a lot the accuracy of the model and of the predictions. So this is it's always trade offs, right? But basically the approximation is valid for a given universe, if you will. And the hyperparameters partially partially determine that universe. And so you definitely want to pay a lot of attention into how these hyperparameters are chosen. ⁓ in PineC we have that function that B Langles and I developed a few, I think, last year or a year and a half ago, something like that, where you pass a bunch of Constraints that you want on your covariance. And and then PyMC will do some optimization under the hood to find the optimized hyperparameters for you for the model. It's basically the same thing as what you would do by hand, but just automated and it's a bit more user-friendly. A bit like the defined constraint prior function that I worked on that I merged like a now a few years ago in PymeC, where you could just tell. PincyO ⁓ hey, I'm looking for a gamma function with ninety-five percent of the mass between this and that. What are the good parameters to do that? Now it's in prelise and it's so much better, and Oswald Martin and and colleagues made that a g awesome package that's so much better than what I coded at the time. ⁓ but that's kind of the same idea. Yep.

speaker-1: That's good. I should probably look at that because I wanted to add a some functionality to automatically determine the length scales based on the input ⁓ of the GPs because you can set the priors. Well that you can kind of set the priors. I in the model I've basically decided on the prior distribution for each parameter, but you can set the hyperparameters yourself. ⁓ and so it has some sensible defaults, but it would be nice to input, for instance, the distance matrix for your sites and then determine a suitable length scale based on the input. ⁓

speaker-0: Yeah, especially since I mean in your case it's a specialized package, so that's that's actually very good because most of the people are gonna use the package for the same kind of models. The difficulty we have in Prime C is that it's a generic package, so tons of people use it for very different use cases, and so we cannot make a lot of very informed default like that. So the the function we have to find the hyperparameters of the the HSGP, you need to provide a length guard, you need to provide A covariance function, the length scale, and like a couple other parameters that I'm forgetting. but you could go one level of abstraction ⁓ because most of your users are gonna use it for occupancy ARU models, and so you can have very very good informed defaults on that already on the length scales. ⁓ so yeah, exactly. I think you can just go on the Pym C source code. Find that function, I think it's called find H S C P hyperparameters or something like that. and yeah, like start from here and and add a level on of abstraction on that. That'd be that'd be very good from a user friendly perspective.

speaker-1: Yeah.

speaker-0: ⁓ yeah, but that that's super exciting. ⁓ I I think it's awesome that you're making access to GP like that much easier. Because as I always say on the show, GPs are amazing. ⁓ they are in itself, you know, even though they are a beautiful mathematical object, ⁓ they are in itself an approximation. Like no no real process in nature is really a GP. ⁓ you're using a GP because you don't know the actual mathematical process underlying the the phenomenon. Otherwise you would just use the mathematical formula, mathematical formula for it. But they are super powerful because they are so flexible. And and I mean there are there are neural networks. And we all know how good neural networks are. And so GPs are just ⁓ infinitely neural networks with an infinite number of hidden layers. And and that's what neural networks converge to and and they are just amazing because of that. But they are a bit complicated to understand from a theoretical standpoint. And there is a high barrier to entry. And and then also the barrier to entry tends to be high for practitioners. Because you're like it's you don't have a lot of packages that make access to it very easy. And also make them available in packages that people are already using for other kind of models. So often they are specialized packages just for GPs, maybe for one kind of GP, you know, the which is like if you don't know the research, you're like, whatever needs that GP? You know. ⁓ so yeah, I think making GPs more accessible, especially when it makes a ton of sense here, that's that that can really be a game changer for practitioners.

speaker-1: Yeah, and initially I didn't even have the option that you could turn them off because the package started off a bit more opinionated, basically, about the fact that we need to accommodate for these repeated measures ⁓ through random effects. And as complicated as GPs are, the way I explain them to people just getting used to them is it's essentially a a random effect where the realizations aren't independent necessarily, but they're dependent on their neighbors. So you're just adding that proximity component of surveys that are closer together in time and sites that are closer together in space, ⁓ maybe more similar to each other than sites that are farther away, which I think is quite easy to wrap your head around that that is useful and relevant. And yeah, I think at at this point you can turn them off, but by default they're turned on. just just out of principle, basically that to ac to account for these repeated measures that we w we should have them in.

speaker-0: Yeah, yeah, yeah. Bold choice to have that by default. I love it. and so something also another thing you have that's that's also very very useful are global local shrinkage priors for variance decomposition of the occupancy and detection ⁓ linear predictors. So can you tell us what this is about? What is variance decomposition in this context? And why is shrinkage an interesting tool here?

speaker-1: So the the intuition around it is that we have a linear predictor with a certain amount of variance and different model components are in a way competing for that variance or at least contribute to that variance. And you might have a bunch of predictors. You could have two predictors, you could have two hundred predictors. ⁓ and then you have random effects and then you have species level deviations because every species has its own response to the predictor and species level random effects. So we potentially have a large number of variance components. And it's actually traditionally you would just put independent priors on all of the coefficients on the species well the species coefficients might be modeled hierarchically with a shared ⁓ scale. ⁓ but the coefficients would probably be independent. ⁓ the variance the variance components of the random effects would be independent. And they've pretty well shown that the prior predictive distribution of these models, as you add more and more predictors and more and more random effects, the prior predictive distribution gets a bit silly. Where you have a huge amount of variance that you are expecting a priori to be captured by your model components. And what the global locate shrinkage priors does is it kind of turns it on its head by saying, I'm going to put one prior on the overall variance of my linear predictor. And the prior distribution for that variance, I've chosen a student T, just because it has a nice thick tail that if you have. A lot of variance in your linear predictor that can easily be accommodated. The key intuition is that you now partition this global variance to your different model components via a simplex. So at this point, you have a global ⁓ you have a global variance and a simplex that you're estimating from your data. And it basically neatly partitions your available, let's call it, variance budget. Across an arbitrary number of model components. So, like I said, you could have literally hundreds of predictors in there, and it should be fine for the model because your total variance doesn't really change. It just gets partitioned between these different model components. And the simplex has an underlying parameter, or it could be multiple parameters, but in this package it just has one. And that parameter is determines the ⁓ sparsity of the simplex. And the sparsity of the simplex basically determines how evenly our simplex is distributed. So when you have a sparse simplex, there's a few components, a few model components that are explaining most of the variance. So let's say you have 20 predictors in there, none of them are important except for one. You expect the simplex to be very concentrated around that single model component. ⁓ which means the simplex is sparse because most of the simplex values are zero except that one which is very high. ⁓ and if you have an even simplex, ⁓ you expect a more even distribution. So this predictor has like 10% and that predictor has 7% and this one maybe 15% of the simplex. That's actually a hyperparameter that's included in the model as well, ⁓ which basically means that the there's a prior over the sparsity of the simplex. Which means that a huge amount of model configurations essentially are possible using the exact same hyperpriors. So theoretically, ⁓ this model should basically run whether you have one species or a hundred, ⁓ zero predictors or fifty, random effects everywhere or none at all. That should all just be handled by having these two priors on the overall variance and the simplexity of composition.

speaker-0: Okay. Yeah, I love it. This is this is super interesting. ⁓ and so to do that what are you what are you using under the hood from the for I'm I'm guessing it's from Stan?

speaker-1: Yeah, so actually I'll just backtrack and give credit where it's due. ⁓ the this basically comes from the R2 D2 literature where they decompose R squared because you can the you can get the explained variance as a deterministic transformation of your R squared. And so what they do is they set a prior on R squared, deterministically transform that to your explained variance, and then simplex decompose that. I skip the R squared because the R squared gets a bit hard to interpret when you have count variables and also logistic. So the occupancy side remains a Bernoulli, so a one or a zero, and the detection rate is a count, so it's Poisson or negative binomial. R squared doesn't work as nicely intuitively as it does for nice Gaussian models. ⁓ so I skip the R2 stuff. I go straight to the variance prior. And under the hood, I just produce a an R a simplex of arbitrary length. And I have to give thanks to the stand developers because they have recently put in functions that lets you simplex transform just a free vector with one less element than your desired simplex, which means that under the hood in this model, it looks at all the variance components. So How many species, ⁓ how many predictors, how many random effects. Based on that, you know how long your simplex has to be. The simplex is then produced under the hood. And yeah, what I use in Stan is a a variance parameter with a student T prior and a simplex with either a deer schleigh or a logistic normal prior for the simplex.

speaker-0: Mm-hmm. Yeah. Yeah, so you're removing one degree of freedom because well, you have to. Yeah. Sounds it sounds a lot like the the zero subnormal ⁓ distribution that that we worked on on on the Pine C sub with Adrian Seabolt.

speaker-1: Exactly. It is, and there's a lot of overlap between zero sum normal and the simplex and both of those constraints ⁓ are getting a lot of use in the package. ⁓ because the package accommodates, for instance, categorical predictors and instead of setting a reference category, I just model those as zero sum normal based on how many categories there are. But every categorical predictor can have an arbitrary number of categories. So you might have habitat type like forest or swamp or desert as a categorical predictors with three levels, but you could also have one with six or seven levels. ⁓ so based on that, it just produces ⁓ a zero sum vector of the appropriate length for every categorical predictor. ⁓ so actually a lot of the stuff that's in this package now has been has implemented really recent additions to the stand language. So the timing's been really good.

speaker-0: Yeah, I'm guessing that that made your life much easier. Because before that was that was much harder to do.

speaker-1: Yeah. Yeah, and it and and it was just really messy. There's the there's still there's still that you have to have a bunch of unconstrained parameters in your parameters block and then use the transform parameters to actually constrain them. So one downside is still that your final fitted object is going to be bigger because it has a bunch of these unconstrained parameters that you need to produce the transform parameters. But ⁓ I mean that's a small price to pay because It just makes it just makes modelling these types of things so flexible.

speaker-0: Yeah. Yeah. No, exactly. That that that's really that's really being a game changer to have the zero subnom constraint in Pine C and in Stan, so definitely encourage people to to use that more. It's great that you're using that under the hood because again, these are complex mathematical objects, in the sense that they're not very intuitive at the beginning. And so being able to just have that under the hood for the users, they don't even have to care. that it's actually using that. It's just like just use that because it's it's better. Like it's gonna make your model better. This is the only thing you need to know. And actually you don't even need to know about it. It's just like that's that's that's what you need to use. Just just do what the doctor said. And and that's gonna make their model better. ⁓ so we're going fast here and for people who wanna have ⁓ background on that, ⁓ the ⁓ at least three episodes they can Listen to ⁓ seventy-four with Adrian Zabold, who is the the father of the Zero Subnormal. Episode 133 with Sean Pinkney, who is one of the ⁓ persons who developed one of the people who developed Zero Subnormal constraint on the stand side. So if you guys want the the background on that and all the cool thing you can do with that, ⁓ definitely listen to that one. Adrian also make a a guest star appearance in that one. And finally episode one thirty four with David Cohns, because in this one ⁓ David goes through and actually demos live the ARR squared prior, which is the R two D two prior, but for time series, for especially auto regressive time series, and it goes into a lot of the concepts that you just ⁓ you just ⁓ mentioned. Matt. So so if people are interested in that, definitely check that one out too.

speaker-1: Just ⁓ a another nice thing with the zero sum stuff is that in this model, for instance, I have species level random intercepts that sort of capture the overall detection rate of a given species in your data set. But then if you have, let's say, a categorical predictor and you model that hierarchically, you have an overall habitat effect to say this habitat has this influence on just your species and a species level offset. realization of that. So there you get into actually the zero sum matrix where it has to sum to zero over the categories, but also across your species, so that your overall effect keeps that nice, well, that zero sum interpretation. And it means that the intercepts in these models are actually the the averages over the species over all of your categorical predictors. ⁓ whereas if you don't zero sum anything ⁓ you know, you get those wider uncertainties you actually have to go and derive those means afterwards in generated quantities. And gets really hard to keep track of. So yeah, the zero sums have just been amazing because of well, identifiability of your parameters in the strict sense. ⁓ but also sorry, I just actually blanked on what I was saying.

speaker-0: That's normal. These are ⁓ complicated objects as I was saying. ⁓ but yeah, basically you are you're talking about the fact that e it can get intricate because these are like Russian dolls where you have basically you might need zero sum normals all the way down because if you're dealing with categories and on top of that simplexes, then you have a a bunch of Some to zero constraints that you need to avoid other parameterization in your model. And so basically what you're doing is getting the headache out for the users, which I'm sure everybody will appreciate. And actually another another concept I wanted to touch on is within chain parallelization. Because that's something ⁓ I think a lot of listeners have heard of on the show, but never necessarily actually used because it's quite new. So in your experience, when does it actually help and when is it not worth the complexity?

speaker-1: Yeah, that's a that's a a funny one. It's it's not really intuitive, ⁓ I find. So just as a background, with within chain parallelization, ⁓ let's say you still want to run four chains. Normally you just run one chain per core. But if you have more than four cores available on your computer, you can actually split the work of the individual chains across different cores. So for instance, if you had forty cores and wanted to run four chains, you could run the work of each of the chains across ten different cores. ⁓ there's a lot of overhead associated with within chain parallelization as I understand it. A lot of copies to go across the different cores. And I've actually found that it often does not speed up the model. I think normally you cut the likelihood ⁓ computations over your sampling units. And ultimately in these types of models or in occupancy models, the sampling unit is the site. So you compute all these likelihoods for all of the species at each different site. And what I have found is the more ⁓ expensive the per site likelihood computation gets, the more benefit you get from it. ⁓ I have not really figured it out. I end up testing it out on the different ⁓ models that I have. Within chain parallelization is optimal. You can just turn it you can just turn it off because a lot of the time it won't speed up. ⁓ Yeah, I'm also not sure how ⁓ what to exactly parallelize because you could just parallelize the likelihood function and do everything else ⁓ inside, for instance, the model block or transform parameters. But I've actually put all of this stuff associated with the sites inside of the parallelization function, which I think means that there's ends up being a lot of copying. So a lot of data structures have to get copied. There's like 30 arguments to the ⁓ parallelization function. So yeah, maybe at some point I will try and figure out if there's a a better way to do it. But ultimately it does come down to trial and error really and how many cores you have. ⁓ I'm not quite sure. I think it is a from from what I've read on the stand page, it's actually it's there is no ⁓ cut and dry sort of outcome with parallelisation.

speaker-0: Yeah, it's also my experience, but this is not really my area of expertise. I'm sure Adrian, for instance, if he were on the show, ⁓ and he will be soon again on the show, by the way, ⁓ maybe he'll have he'll have ⁓ things to add. but another thing, another concept you're using also in your package is Monte Carlo integration to deal with high perito K values. when you have observation level random effect. So listeners may have seen that in their models when they are using RVs, they may have a warning that there is a high perturkey value in there that makes the the comparison with Lee Wano cross-validation not trustworthy. ⁓ maybe they don't even know what that means. And so this is related to that So can you can you walk us through what that problem looks like in practice and how the Monte Carlo integration fix works?

speaker-1: Yeah, I'm I'm actually with your hypothetical users because I do not have a good understanding of what Leave One Out Cross Validation is actually doing and don't have an intuition around the hyperitocays. ⁓ I do understand when they come up. And in an occupancy model, like I said before, the lowest level where the likelihood factorizes, which is the level where you want to compute. the low likelihood for leave one out cross validation is the site. And as I mentioned before, ⁓ we have these site level random effects because, well, we have the repeated measures at this unit as well. So I feel like they should be there, but it means that when you have a ⁓ random effect or a parameter associated with the level of the observation, you end up getting these high Perito K values. And it's something that I've been ⁓ encountering or struggling with for a lot of ecological models because well we tend to have random effects at the level of the site when you want to capture ⁓ variation between the sites. And I got the tip from Aki ⁓ Vitari or Vetari, I don't know how to pronounce his last name, ⁓ on the Stan Discourse and he said to use Monte Carlo integration, which is essentially ⁓ In generated quantities for each post or for each MCMC iteration, within each draw, you would do a bunch of draws, or you would compute the likelihood many times for what is essentially a posterior predictive random effect. So let's say you have 50 Monte Carlo draws within each MCMC draw, you do you compute the likelihood 50 times. with a different value for that random effect, and then you average it out within each draw. ⁓ and this is still fast because it's happening in generated quantities. So ⁓ it's well stuff that happens in that block is fast. So it doesn't add that much to the runtime. ⁓ but I found it very much reduces the pareto case. So ⁓ I've gone usually I'll get somewhere of I don't know sixty or seventy percent of my observations. have high parido K values without the Monte Carlo integration. And it might go down to 10%. And those probably genuinely reflect ⁓ poorly fitting observations or poorly fitting sites. ⁓ so it's just one I don't know if it's the most principled way. I think it's kind of a brute force way to go about it. But ⁓ I'll take anything to ⁓ yeah reduce those high paridoca's ⁓ and make my model comparison a bit more trustworthy. Because yeah, the the key thing is when you have observation level parameters and you get those high peridocays, ⁓ it's not necessarily that your model is bad. It's just that ⁓ leave one out cross validation isn't working well.

speaker-0: Yeah, yeah, exactly. and in in and in I mean a lot of time, especially with count models, you will see these warnings and that's because you have some highly influential observations. And so that means that if you re if you remove that observation, then the landscape completely changes because that observation is is very highly influential and it determines a lot of what the model is saying about the phenomenon. But if you remove it and if you have a very low sample size, which can happen with count data because it's it's discrete, if you just have, you know, like one count of that data, as I understood, then that's a highly influential observation. And if you remove it, that means leave one out cross-validation is gonna change a lot depending on which observation you remove. And that can change that can change the results a lot and that's not what you want because that mean the the number you see is not necessarily very reliable and reproducible. That's also yeah.

speaker-1: Have you found have you found this in your experience to be ⁓ common with count models more so than co continuous data?

speaker-0: Yes, but that's ac anecdotal evidence. Like I'm sure Aki would be would be a much better interlocutor than myself to tell you the the actual scientific truths about that. But my experience since I work also a lot with count models is that often part okay is is problematic for for count models. I think I remember

speaker-1: Also for ⁓ negative binomial? Because I feel like the Poisson is just like completely restrictive and unrealistic, especially for ⁓ for ecological data. ⁓ have you found that high freedocays tend to persist with negative binomial as well?

speaker-0: Honestly that I don't know. Like I don't remember. I would need I didn't pay enough attention, I think. ⁓ what I know is that yeah, all I almost know and I almost never can use a pure Poisson distribution. Like it's always almost always too restrictive because of the mean equals variance constraint. It's just like the the likelihood's just too constrained and not variable enough and and then I have to use a gamma Poisson.

speaker-1: Or I actually ⁓ I encountered it recently in a wild data set where the Poisson was preferred to negative binomial and ⁓ it was an incredi it was an incredible feeling.

speaker-0: Yeah Yeah, that never happens. Yeah yeah.

speaker-1: No. Yeah. I don't expect it will happen again.

speaker-0: No, yeah, yeah. Almost always I mean I always try the Poisson because like it's the scientific process and you know, I have to justify ⁓ why I'm using ⁓ the the gamma poisson. But like almost always I I can't even remember a case, but that doesn't mean never happened, where I had to I could use a Poisson. Like it it doesn't pop into my mind.

speaker-1: That's where I quite like the I don't know if you've come across the PC that penalized complexity priors. where if you consider the Poisson to just be a special version of the negative binomial, you can set a prior that essentially penalizes how far it differs from a Poisson model. So you can kind of have your cake and eat it too. You can if there's no evidence in your data to suggest that over dispersion is present. you can just fit a negative binomial model that basically looks like Poisson. ⁓ but if there is evidence, then there's plenty of space in the prior for that to occur. ⁓ and that's actually the default prior that's in the negative binomial over dispersion parameter. ⁓ because the the R package the the R package allows you to ⁓ have over dispersed data. So you instead of a Poisson you fit a negative binomial. And I think it was Aki that ⁓ did some work to look at which values for a an inverse gamma on the over dispersion parameters most closely mimicked the PC prior ⁓ on it. And that's the prior that's now default in BRMS for negative binomial ⁓ models. So it's an inverse gamma ⁓ zero point four, zero point three. So That's what's ⁓ that's what implemented by default.

speaker-0: Wow, that's really cool. That is the default. yeah, so penalized complexity priors, I love them. I use them all the time for GPs, actually. ⁓ especially for amplitudes and length scales. ⁓ a bit more for amplitudes because most of the time when I'm using GPs, I'm using HSGP and the length scale is actually fine with HSGP. Because if you're using SJGP, most of the time you have a lot of data. And so you have a lot of data which can inform the length scale, and so you're fine. ⁓ the amplitude is a bit more complicated, and so I always use a penalized complexity prior on the amplitude, which is awesome because it's very intuitive to define. and the on the length scale I've tried it, but it's actually harmful. It's much better to have like a more diffuse prior and then have the data move the the prior on the length scale if you have a lot of it. And so yeah, like and and then the other PC Priors we need to add them in Prime C. I know Bill Engels has been doing a lot of work on that and he's actually being the one opening my eyes to ⁓ the awesomeness of PC priors. That's thanks to him that I use them in GPs. And ⁓ we have a project to add them into Py MC, but we never we never get to it. if anybody wants to work on that with us and actually get that across the finish line, please contact me. ⁓ but yeah, we need we need to add that to Pine C because it's extremely helpful. And and I know Bill has been working on that. He's been working on some GP stuff also, so he'll come on the show very soon. I've been working on on Having him on the show, so that will happen very soon. ⁓ but in the meantime, and and while we're working in in adding PC Priors to Prime C out of the box, you can already do that. But ⁓ you can already listen to episode one thirty-six ⁓ with Janet Van Niekerk and Havard Rue. They work on Inla, so Nestidla Place approximation. And in the episode we talk exactly about PC Priors. And so if you just wanna hear about that, ⁓ I think it's in the in the chapters, you can just jump into that. But it's a it's a fun episode to listen to. And yeah, they are they already have a bunch of PC Priors coded into the in lab package that's both for R and for Python. So Yeah, for ⁓ for people interested in that, feel free to to look at that. I'm sure you'll you'll find some very interesting interesting nuggets of ⁓ codes and concepts in there.

speaker-1: I think the PC priors are actually a bit orthogonal to the variance decomposition priors because the scale or the amplitude of the GPs is just one of the variance partitions that's in s that lives on this simplex. So I don't even know how I don't think it would be possible to even put a PC prior on there because that variance partition is just absorbed into absorb isn't the right word there is just one of the components that lives on that simplex that decomposes the total variance.

speaker-0: Yeah. D that I mean I've never had to work on that case, but I agree with your with your ⁓ prior that ⁓ yeah, it's probably mutually exclusive. Yeah, it's a trade off you have to take depending on what you're mostly interested in in in that case.

speaker-1: I think so, yeah. I do enjoy or I enjoyed seeing that when I ⁓ wrote a function that lets you set the hyperparameters for your priors, that you have a potentially very complicated model with predictors on occupancy, predictors on detection, random effects, et cetera. But the actual list of parameters that you need to set priors on is quite short. Whereas if you didn't have these decomposition priors, you'd potentially have to set priors for every coefficient, for instance. or the variance of every random coefficient slope type thing. and if there's one thing that this package has sort of writing this package has sort of taught me, it's just the the choices you have to make when you write a package that's sort of general purpose and lets you fit arbitrary number of species, arbitrary number of seasons and predictors and stuff like that, like the the sort of Classical approaches of just having a well a regression model with maybe a couple of coefficients is really such a narrow version of what's possible with probabilistic programming. ⁓ and a more complex model really forces you to think about how to set principal priors that are essentially generalizable to a huge number of model configurations.

speaker-0: Yeah, yeah, exactly. Definitely. It's ⁓ it's incredible to see how models are actually incredibly customized, even though you think you're working on a on a classic ⁓ classic simple model. in so I'm gonna start to have ⁓ playing us out here, Matt, because I don't wanna take ⁓ you know your whole ⁓ lunchtime. But we've we've been

speaker-1: Sounds good.

speaker-0: Yeah, I'm very happy because we've gone through a lot of the technical topics I wanted to touch on with you. But a question I always love to ask people who work on packages is when does Okaru break down or underperform? You know, what kinds of data sets or questions is it not the right tool for? And what would you tell someone before they reach for it?

speaker-1: I think there's a I wouldn't say there's a huge learning curve. It's a complicated model that I've tried to make as accessible as possible with just a couple of functions. there is a big part is the chosen survey interval, which is basically how do you bin your data? Do you group it by day or week or month type thing? So if you have years and years of data and you your observation level is how many count how many individuals or detections you got of each species per day, that model will be incredibly slow to fit. Whereas if you have longer like multi-year data, you might group by a week or two weeks. So that's a bit of a ⁓ where big data can still be accommodated by choosing that survey interval. But overall, if you have lots of data, lots of sites, this will probably end up being quite slow. And I I work on a MacBook and I've got a Apple silicon and it's not super slow or anything, but I've seen people try and use it on a not as fast a computer and it's already quite slow. I mean, these people still have a lot of these data sets will have tens of thousands of observations even after a thinning interval. ⁓ so that's ultimately still always going to be a lot of data, especially when you have these random effects. So yeah, I think The thing is I don't know what would otherwise suffice and be a lot faster. because even if you do collapse it down to ones and zeros, like you still end up having a lot of sites and a lot of surveys. You just have a different observation model and maybe not as many random effects. ⁓ so yeah, to I'd say the biggest bottleneck will be really large data sets ⁓ and also potentially not ⁓ ⁓ reducing the resolution enough of the data. So if you if your model runs really slowly, you might even group it by month and just have how many counts per month did we get? And you know, that it actually doesn't make any difference at all if we're getting occupancy. So ⁓ if you're only interested in occupancy, you can just group you can have huge grouping periods and it still gets you occupancy. So ⁓ yeah. Probably speed, I would have to say, given that it's ⁓ stand and potentially it can be quite complicated with with many sites and slow to fit, I would say.

speaker-0: Hm. Yeah. Yeah, that makes sense. ⁓ anyway, grouping and grouping always your friend. I had the same could have the same issue when ⁓ when I was working on electrol for testing models. So I definitely understand what you're saying. And ⁓ actually for for an ecologist or an applied researcher who've been listening to us right now and who wants to start using Okaroo on their own.

speaker-1: Yeah.

speaker-0: Where do you point them? And and w also what's next for you for the the roadmap on ⁓ for the package?

speaker-1: So the the you can download the package on ⁓ GitHub. and the package I don't think will be on CRAN because it uses command stan R in the background, which has the latest version of Stan. So for instance, all this zero sum and simplex constraining ⁓ wouldn't be available on R Stan, but command stand r isn't on CRAN and there's no plans to put it on CRAN. ⁓ Roadmap for the package is I'm nearly finished with a the multi-season version. So that's where you start to get into like dynamic how occupancy can change from one season to the next as a Markovian process. So that's where you get it really in the hidden Markov side of things. ⁓ once I have that down, I have a lot of work to do myself using this package because ⁓ for my my postdoc, I mean my main motivation for developing it was that I didn't feel that the current methods were able to fully exploit the information that's actually present in these data. And when people want to know like, hey, how good is our camera trap data, I can't really answer that if I have methods that only touch on a fragment of the information that's actually contained in them. ⁓ So I hope to be ironing out a lot of kinks as they come up. Like it's unbelievably easy to introduce bugs into code. ⁓ I've come to discover making this R package. ⁓ yeah, I'd like to get smooth sailing and hopefully get people using it and tell them what they need, what they're looking for, and obviously bugs that they're finding. ⁓ and I'd just like to keep improving it to expand on it and making it more ⁓ more customizable as well.

speaker-0: Yeah. Yeah. Sounds ⁓ sounds like a great roadmap, ⁓ very very exciting and ⁓ and hope you'll you'll manage to to do that in the in the coming month. You seem ⁓ you seem to be ⁓ on a good cruising speed so so really rooting for you. and well let's call it a show, Matt. ⁓ I think it's ⁓ it's a good stopping time. We've been ⁓ we've been covering a ton of ⁓ of topics. I'm very happy about that. Before letting you go, of course, ⁓ you listen to the show, so you know I have to ask you the last two questions. I ask every guest at the end of the show. If you had unlimited time and resources, which problem would you try to solve?

speaker-1: ⁓ yeah, so firstly, thanks for having me and ⁓ yeah, it was ⁓ it was great. The problem I'd like to solve is ⁓ habitat loss and just the destruction of the natural world around us is one hundred percent the thing that I would like to that I would like to solve the most. ⁓ I don't know if there is a clear solution. ⁓ I'm fairly cynical or pessimistic about it, but yeah, if I if I had infinite time and resources that would be ⁓ what I would like to work on.

speaker-0: Sounds perfect. ⁓ sounds ⁓ yeah, like a like a very, very important problem indeed. And second question, if you could have dinner with any great scientific mind, dead, alive, or fictional, who'd it be?

speaker-1: Yeah, I had to think about this because I knew I knew the question was coming. And I think it would have to be Jeffrey West, who's a physicist that wrote a book called Scale, ⁓ which is probably my favorite book of all time and changed the way I look at the world. So I think ⁓ having dinner with him would be awesome and I actually think it's worth ⁓ worth people listening or reading his book. It's ⁓ it's a fascinating fascinating story.

speaker-0: Okay, yeah. ⁓ I I didn't know about him but ⁓ that sounds that sounds interesting, so thanks for educating me, Matt. And ⁓ well on that note let's call it a show. Thank you so much again for coming in the show, Matt, everybody. There will be a lot of of show notes for this one, so make sure to check them out if you wanna dig deeper and Matt. Thank you again for taking the time and being on this show.

speaker-1: Thank you.

speaker-0: This has been another episode of Learning Bayesian Statistics. Be sure to rate, review, and follow the show on your favorite podcatcher and visit learnbasetats.com for more resources about today's topics, as well as access to more episodes to help you reach true Bayesian state of mind. That's LearnBaseTats.com. Our theme music is GoodBayan by Beba Brinkman. Fit MC Las and Mega Rand. Check out his awesome work at bebabrinkman.com. I'm your host. Alexandora. You can follow me on Twitter at Alex underscore Andora like the country. You can support the show and unlock exclusive benefits by visiting patreon.com slash learnbase tax. Thank you so much for listening and for your support. You're truly a good baby and change your predictions after taking information. And if you're thinking I'll be less than amazing, let's adjust those expectations. We show you how to be a good daisy. Change calculations after taking fresh data. Those predictions that your brain is making. Let's get them on a solid foundation.