| Question | Answer(s) |
|---|---|
| Gilbert is still doing the fieldwork every year :exploding_head: | Field work is life! Good for him! |
| Are there still confounding parameter estimates in a fully time dependent model in a Bayesian framework for CJS? | Yes, this is still the same. At least as long as time-variation in parameters is not constrained in any way. |
| NA | will talk about param redundancy later on :slightly_smiling_face: |
| Are there tips for reusing as much code as possible to limit redundancy when writing many versions of models <http://phi.pt|phi.pt> phi.p. Phitpt etc and all the other covariates we may want to test?? | Yes, you could define gamma and omega individual- and time-specific, and do the same in the likelihood. By doing so, you'd have to change just the parameters as you see fit. The thing is that by doing that, you have to keep big arrays along the way, which makes calculations slower. has probablement comments on that. |
| NA | says |
| is the prior for survival always between 0 and 1 or do you sometimes bound it based on previous research for the species i.e between 0.6 and 0.9? | Yes, we'll see an example later on. This is a situation where the beta distribution is useful. |
| In nimble we support model-definition-time if-then-else statements. What that means is that you can use if-then-else where the if condition uses a variable (or expression) evaluated in your *R environment* to decide which piece of code to include in the model definition. It does not create if-then-else in the model, but rather creates different models depending on the if-then-else condition. It would be cumbersome for many model variants but can be a handy way to support a few variants in the same model source code. | Oh sorry, this was a reply to but I forgot to click "reply in thread". |
| How's model selection when we model with random effects? Does model selection lose sense at all then? | I think it's a complicated topic. From a Bayesian perspective, one can do model selection at different levels of the model. Model selection conditioning on random effects may not necessarily make sense for choosing variables at a higher level of the model, and vice-versa, for example. Here is one article that may be of interest: |
| NA | Thanks, Perry! |
| Not sure if this will come up later but how would you code for variation in Phi/p for specific occasions? e.g. for an oil spill year where you know there will be an effect on survival? | For a case like this, the best would be to make a dummy variable that takes 0 in all years without oil spill, and 1 in years with oil spill. You can then include this like any other covariate. |
| NA | If you need more than 2 levels (i.e. on/off), you can still make a variable with integer entries representing different levels, but could then estimate intercepts specific to the different index levels. |
| NA | Thank you - I'll have a go! Sorry if this is a silly question but just to confirm you using the phrase dummy variable which do you specifically mean? You mention on/off is this related to detected/non-detected or another matrix (e.g omega)? |
| NA | I'll write it as an example. Check back later :wink: |
| NA | Thank you! Sorry I've been working with a complicated dataset in Mark so know how to do things there but trying to make sure I get how Nimble accepts it! |
| NA | Binary example (on/off covariate, e.g. oil spill vs. no oil spill): dummy <- c(1, 0, 0, 1, 0) In your model: # Constraint for(t in 1:Tmax){ logit(phi[t]) <- logit(Mu.phi) + beta*dummy[t] } # Priors Mu.phi ~ dunif(0, 1) beta ~ dunif(-5, 5) |
| NA | 3+ level example (categorical covariate, e.g. bad oil spill, moderate oil spill, no oil spill): dummy <- c(2, 1, 1, 3, 1) In your model: # Constraint for(t in 1:Tmax){ phi[t] <- Mu.phi[dummy[t]] } # Priors for(x in 1:3){ Mu.phi[x] ~ dunif(0, 1) } |
| NA | I recommend checking out "Bayesian Analysis using WinBUGS" to get familiar with the basics of implementing different covariate effects. |
| NA | Fab - thank you again, I'll have a go with testing that out and yes, actually get a copy of the book rather than sneak a peak of my colleagues copy! |
| NA | You won't regret getting that book :wink: |
| NA | Hi can I check on the author of "Bayesian Analysis using WinBUGS" that you mentioned? The book I was thinking of yesterday has a different title and Google didn't direct me neatly to your one? |
| Checking, out of curiosity - I guess I should know this having heard about this data set so many times but hey - in the dipper data wing length is measured when (since wing length would change over the life time of the animals?) | Haha, I simulated this covariate for the purpose of illustration. Otherwise, you're right and we get time-varying individual covariate for which we need to do something for individuals that are not captured. |
| NA | aha - many thanks! |
| NA | More soon on time-varying individual covariate. |
| NA | Not the case here if has simulated the values but in real data wing length could be used to sex this species. How much bird wing length changes over an individual's life varies quite a lot between species and isn't well studied as far as I know... |
| Is overdispersion (the ominous c-hat) an issue with Bayesian estimates? | It's a general issue I guess, often caused by non-independence among individuals, which you can address with appropriate distribution (e.g. negative bin vs poisson) or random effects. We'll briefly talk about goodness-of-fit tests tomorrow. |
| NA | Cool, thanks! |
| On a similar note. Should we look at goodness of fit of the dataset prior to the analysis? | Yes, more about that tomorrow. |
| Does it matter if the population is open or closed? | Yes, here we assume open populations, hence the estimation of survival. If your pop is closed, then you may assume survival is 1 :slightly_smiling_face: |
| Ironically just got briefly distracted by the capture and marking part of capture-mark-recapture/resighting! Woodpigeon (Columba palumbus) freshly ringed by my partner with colour ring white CH added for a BTO retrapping adults for survival (RAS) study. As its easier to add the colour ring with 2 people, there are now feathers in my home office... | jealous :rolling_on_the_floor_laughing: |
| NA | Tempted to use this picture as the cover of our website. Meta-capture-recapture :slightly_smiling_face: |
| NA | Is: from field work to real data to theory and statistics, perfect illustration! I have never introduced a computer or a model to a kittiwake. Computers lack salt (according to seabirds). |
| NA | Ha I would have got a better picture but I was still trying to pay attention to what was saying! |
| NA | I love the picture thanks for sharing ! Theory and practice finally re-united :slightly_smiling_face: |
| NA | Some are lucky (you) and others are less (me - 1 day of fieldwork in the last 3 years...) |
| In this kind of models, we do not need to re-sample the same spots (camera-sttations, etc) over the years/seasons, right? | Yes, resampling is key to separate non-detection from actual death. |
| NA | I mean, do i have to re-sample the same specific habitat/area/sampling station? or just if i have a re-capture (no matter if i re-sampled different spots next season) in a different spots with different sampling effort/sampling design? |
| NA | Oops, sorry. Yes, you need to resample the same station. If you change the sampling design (increasing the sampling area for example) then things you get a different population under sampling. Different sampling effort (less people doing the recaptures some years for example) can be accounted for in the detection process by having time variation, and possibly a covariate to explain time variation. I hope it answers better your question :slightly_smiling_face: |
| NA | Technically, you do not need to re-sample the exactly same location. However, if your sampling locations and effort changed, you have to incorporate this into detection/recapture probability (i.e. you need time-dependence). Additionally, if you changed locations, think carefully about what this means for your survival and detection estimates. Survival is always confounded with emigration (unless you model them separately with appropriate data and study/model design), and if you change locations, that may introduce bias. |
| NA | Thanks for both :slightly_smiling_face: |
| NA | Chloé, now that you have finished this lecture (transition stimation)... I have 7 years data of jaguars from an specific protected area; but the first years, we sampled in waterwholes (because the projects aims), and the last years, we sampled both waterwholes and roads/trails. From the 7 years of data, I have very nice sex and identity data from all jaguar individuals. Of course the sampling design and effort is different every year (dry season), but always from the same area (number of waterwholes, number of trail-stations, etc., but always during dry seasons every year). May I use these models (survival and transition estimations) to compare between sexes? How do I may incorporate this data into detection/recapture probs? |
| what should be the minimum dataset size to correctly run such analyses ? (number of individual histories) | This depends on the model and on the number of years you surveyed for (in addition to number of individuals). There is no simple guideline here. You need to try it out. If you are fitting a way too complex model for your amount of data, you will not get convergence/good estimates. That then means you need to simplify your model. In general, it's usually best to start with a simple model, and build complexity stepwise. |
| NA | :+1: |
| NA | You may use simulations to investigate and get an idea of statistical power to detect some effects or the precision of your estimates. |
| Is WAIC also a recommended procedure for deciding to include covariates (e.g. sex, age, other individual characteristics) in survival models? | Yes, I think variable selection is viewed as a particular kind of model selection. Note that for Bayesian variable selection, the so-called indicator variable method is an option, and nimble provides "reversible jump MCMC" samplers, which is a fancy way to say it can handle that approach well. |
| NA | It's important to keep in mind that WAIC (just like any other whatever-IC, p-value, etc.) is not a panacea. Making mindless binary decisions based on one metric or another is always a bad idea. Before fitting any covariates, think carefully about what you would expect to have an effect (based on biology, sampling regimes, etc.). Fit that model and check the posteriors for your effect. The degree to which your effect posterior overlaps with 0 will tell you how much evidence there is for this effect. |
| As I have experts at hand...I have a question :grin: basic CMR makes a closure assumption. But what if the population is exposed to immigration only (e.g. a park where animals can be introduced). As animals enter the analysis at first encounter, does this mean that CMR estimates of survival and detection prob are unbiased? Can this question enters the competition for the best stupid question of the day? | Good morning! Assuming that there is only immigration, no emigration, CMR estimates of survival and recapture probability will not be affected. That is because CMR models typically condition on first capture. |
| NA | Too early for me :sleeping: |
| NA | In your case, it may be worth thinking about whether number of immigrants and/or population density could affect survival or catchability (biologically) though :wink: |
| NA | Ouah, a reply in less that one minute (from awake-experts)!!! Thanks Chloé. And good point about true biological effects of increased density. |
| Good morning :slightly_smiling_face: So, Sarah is showing us to model a parameter (phi) as a function of covariates. There is no error parameter in the regression, what about that? | Good morning :sleeping: You're correct, we're assuming that temporal variation in survival is fully explained by water flow. Easy though to add a time random effect and account for non explained variation, wait for it :wink: |
| *"More of a comment than a question":* Personally, I always specify time-dependence on e.g. survival using the link function (logit) twice: logit(phi[t]) <- *logit(Mu.phi)* + beta * flow [t] The advantage is that Mu.phi then is the "average survival probability" and giving it a Uniform(0, 1) prior as we did so far is still possible. |
Thanks, I was going to ask something along those lines having used jags to make similar models! |
| NA | Yes I always do it this way, not just with survival probabilities. Makes book-keeping much easier, especially with large models :wink: |
| in the line logit(phi(t)) etc is there possibility to include an offset in the right part, as in R GLM ? | Technically, you can do whatever you like on the right side. I do not really see the sense in adding some arbitrary offset, however. This is a different story of course if you include intercept-offsets based on categorical covariates. That can make sense. |
| NA | I was thinking to an "explanatory" factor like sampling effort. We could not be interested by estimating the coeff but we have to account for it in the model. No ? |
| NA | So this depends on whether or not you have data on sampling effort. If you have not measured it but have to assume it varied over time (= a nuisance parameter), then you include it as a random effect as has been explaining just now. |
| NA | If you have measured sampling effort, you would still include it as either a continuous or categorical covariate. You will estimate a coefficient, even if you are not interested in it. |
| NA | No discrimination of coefficients!!!:stuck_out_tongue_winking_eye::no_entry_sign: |
| once the models are fitted, how can we make predictions? | You can make predictions using the parameter estimates you obtain and the equations as described in your model. For example, let's say you estimates survival in your model as in the example: logit(phi[t]) <- beta[1] + beta[2]*flow[t] You get posterior distributions for betat[1] and beta[2]. So in R, you can then simulate whatever flow values you want to make predictions for. Let's call this vector flow.sim. It has X elements. Prediction then goes analogous to the equation in your model: for(x in 1:X){ phi[x] <- plogis(beta[1] + beta[2]*flow.sim[x]) } (plogis() is the inverse of the logit) |
| NA | Ideally, you are going to make the prediction above for each sample in your posterior distribution. That way, you also get a posterior distribution for your prediction (= uncertainty). |
| NA | presumably you can do this directly in NIMBLE (analogous to JAGS), by feeding your flow.sim vector (for example) in as data to the model? Thus you get the posterior distributions as model output |
| NA | you can do this if you really want to. However, for efficiency/speed/memory, it's best to leave any additional calculations out of the MCMC. Predictions can be easily done using the joint posterior distributions, so there is no reason to burden your MCMC with this. (Also, it would suck if you later want to do your simulation slightly differently, and have to re-run the whole MCMC for it...) :wink: |
| NA | Thanks a lot |
| In the initial values code: initial.values <- *function*() list(beta = rnorm(2,0,1) We draw 2 values because we need one value for each beta? And if we had defined a model also for detection rate, and therefore had more betas, would we draw more values here? |
Yes, you will need as many values as beta has index values inside your model. |
| For the random effect using Nimble, is it still correct to specifify it using "tau" (instead of directly "sd") ? | Yes, that still works. The default is still tau, but you can do standard deviation instead by writing "sd = ...." |
| NA | ok, thanks! |
| NA | Besides being more natural imho, the "sd = " parametrization is quite useful to use non-centering a trick to improve convergence, more in the live demo to follow. |
| NA | Yes, sd is easier to interpret than tau :) |
| I think that I didn't get it. When you refer to a 'random effect', is it in the same way as we refer to random effects on linear mixed models? Or it has a different definition? I think I mixed things around :sweat_smile: | Same random effect as in (G)LMM yes. |
| Stupid question ... there is no "probit" function already defined in Nimble? | Looks like there is. See manual: |
| NA | If there was not, you could easily define it as a nimbleFunction, that you could then include in your CMR model :wink: |
| NA | Thanks :smiley: |
| Right, this is possible a stupid question, but it's a nightmare when running models. What about if you don't have the sex of all individuals? will the model still run or you need to add something like na.action = na.exclude? | For this one, it's important to think a bit practical first. If only a very small, random set of animals are missing sex, the pragmatic solution is to just exclude non-sexed individuals from analyses. |
| NA | If animals in one life stage are often missing sex (one common one is chicks/pups that cannot be sexed before they are adults), then you may want to include a sex ratio in your transition probabilities (more on multistate models in the next session). |
| NA | If a large junk of animals across the dataset are missing sex AND you know that sex has a strong impact on survival and/or recapture rate, then you can treat sex as a "partially observed covariate". If you do that, you pass your sex vector including the NA values and you have to specify a likelihood for sex assignment (likely involving a binomial or multinomial likelihood and a sex ratio). Note that any NA data values need initial values. |
| NA | The problem is that if you are testing several covariates, if you remove all the ones with missing values you end up without data :slightly_frowning_face: |
| NA | Yes, in that case, you may want to look into covariate imputation. See my reply above :slightly_smiling_face: |
| NA | Yes sorry I hit reply before seen that reply |
| NA | I've found its worth running the models with just the subset that are sexed (eg adults only) to see if there is a potential difference - so far I've been lucky and there hasn't been but at some point I won't be and will have to consider the above! |
| NA | listen up now for multi-event models. This may also be relevant for your issue with partially observed sex :wink: |
| Does including the back-transformation lines within the model notably increase running time? Maybe with more than just 2 parameters to back transform | Not notably I think. You could benchmark it :slightly_smiling_face: |
| NA | Ok thanks, I'll give it a try then :)) |
| I can't seem to be able to load the dipper dataset. In which package is it supposed to be? | you can download the dataset here: |
| NA | dipper <- read_csv("dipper.csv") > |
| NA | that doesn't work for me |
| NA | What's the error message? |
| NA | Error in read_csv("dipper.csv") : could not find function "read_csv" |
| NA | where should be the dataset? |
| NA | in nimble? |
| NA | You need to specify the path (say my_path) to your file, and use read_csv("my_path/dipper.csv") |
| NA | library(readr) |
| NA | Oh yes, is right, load the tidyverse library first, with load(tidyverse) |
| NA | Ahaha. I finaly found dipper.csv in the github repository. ..I had trouble with tidyverse too indeed. I do have read_csv now. Thanks! |
| NA | Cool. |
| Slow start with no coffee :exploding_head: and coming late to the party. Can someone remind me what the line below (slide 192 with the model hmm.phiflowp) is doing? It sets the z's per animal - but what are the z's again (state, dead or alive?)? The loop sets the first time first[i] that animal i animal was seen to 2? Consider animal 1, and first[1] is 3, i.e. animal first seen at time point 3. A bit confused where we define the initial value of say z[1,2] and z[1,4]? | States z are alive (z=1) and dead (z=2). All individuals are alive at first capture so z[i, first[i]] = 1 for all i. We use a placeholder with the dcat(delta[1:2]) with delta <- c(1,0) for more complex models to come, but you could use z[i, first[i]] <- 1. |
| NA | aha - I think it was the placeholder that was confusing me. Plus I had forgotten that all animals are alive at the start dah! |
| NA | something still confusing me. Maybe a confusion in R/nimble notation translation... I thought that z is a matrix with as many rows as animals and as many columns as time occasions? Is that correct?We are assuming for z[i,first[i]] is 1. Where do we define what we assume for initial states for z[i,first[i]-1] or z[i,first[i]+1] say? |
| NA | z[i,first[i]+1] is drawn in the dcat distribution for z(t)|z(t-1), z[i,first[i]-1] is not defined as it's not used. Does it make sense? |
| NA | Thanks, sorry for being daft, there's something I am missing here. I think I am getting confused with R notation. When we have dcat(delta[1:2]) this is the same as dcat(c(1,0)), but that gives me an error in R since dcat needs the probabilities as arguments? |
| NA | FFS... ok I get it, you are defining the probabilities, |
| NA | so dcat(c(1,0)) means a categorical dist with 2 states, fist state P=1 and second state P=0 |
| NA | presumably this works because it's inside the "nimbleCode" function and so it treats dcat differently than it would treat it in R? |
| NA | so that line really just defines the initial state, which as you say is undefined for each animal before it was detected for the first time, and then the Bayes machinery takes care of generating the following states at each iteration |
| NA | am I getting there? |
| NA | Yes, you're definitely getting there. I wonder whether we should forget about dcat for the first lectures... |
| NA | That might be something to consider. I think in general I would avoid placeholders in the sense that they confuse folks if strictly not needed - you feel like "why is this thing here that loos weird when a constant would be all that is required, what am I missing?". Having said that, this has been the clearest course I have attended in the last decade :wink: so I really feel that you have done some much (for free!) for a community already that only compliments are in order! |
| NA | Thank you, means a lot coming from you. There is much room for improvements! I'll think about this placeholder thing. Maybe something to ask in a feedback form. Btw, do you use feedback forms for the distance sampling workshops? If yes, would you be ok to share it w me so I can get some inspiration? |
| NA | We do use feedback forms. Let me check where I can find them. I am in a course right now and I should be paying attention to a person called :rolling_on_the_floor_laughing: |
| This question could fit in any channel, but here we go; it's a basic, but important one. Are there variations of these models we're talking about that are relevant for non-individual-based data? I.e.if you don't know the identity of an animal (e.g. is not marked naturally or artificially), should you use different models for demographic parameters estimation (which ones)? :thinking_face: | Individual identification is indeed key to the application of capture-recapture models to estimate survival, dispersal, recruitment, etc. For unmarked individuals, you may use N-mixture to estimate abundance and their extension to open populations to estimate survival and recruitment. For fecundity, GLMs are great. You may even combine all these models together with integrated pop models :slightly_smiling_face: |
| NA | Thanks! :+1: |
| The dcat distribution still feels a bit magical to me, in the code you have: ```z[i, j] ~ dcat(gamma[...])``` Which if we replace the indeces gives us: `1 ~ dcat([phi, 1-phi])` I can't get my head around to understand how you can link one value on the left with two values on the right. Checking JAGS documentation (see page 19 in <https://onlinelibrary.wiley.com/doi/pdf/10.1002/9781119942412.app1>) did not lead me anywhere. |
Yes. What I find useful is to play around with the nimble::rcat() function to better understand what the categorical distribution is doing. |
| NA | The Wikipedia page is helpfull too I think |
| How should I code in the model two discrete individual variables with nested indexing? If I use beta1[sex[i]] + beta2[age[i]] I get weird trace plots because it has problems with the intercept (reference group) I guess. It works fine having only one discrete variable (e.g. sex). | It should work with two (or any number, really) of discrete covariates. You have to make sure to properly specify the prior for the different levels of intercept and slopes. Pay particular attention to the fact that your differenct discrete covariates may not have the same number of levels. |
| NA | Also, you will only be able to separately estimate beta1 and beta2 if your data contains enough occurrences of all different sex-age combinations. |
| NA | Hi Chloé, Thanks! Sex has two levels and age has 3 levels. So that should be fine. But I will check if the sample size for each sex-age combination. Maybe that is the problem. Thank you! |
| NA | Hi , so indeed there are maybe too limited data for the second age class in females? |
| NA | That's not so bad. It should still work I think. I recommend running two separate models: one with only the Age effect, and one with only the sex effect. Check that both of those work fine. |
| NA | Ok thanks. Just started the models! |
| > mcmc.phip <- nimbleMCMC(code = hmm.phip, + constants = my.constants, + data = my.data, + inits = initial.values, + monitors = <http://parameters.to|parameters.to>.save, + niter = n.iter, + nburnin = n.burnin, + nchains = n.chains) defining model... building model... setting data and initial values... running calculate on model (any error reports that follow may simply reflect missing values in model variables) ... checking model sizes and dimensions... checking model calculations... model building finished. compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details. Error: Failed to create the shared library. Run 'printErrors()' to see the compilation errors. |
You have an install problem I think. Have you installed the latest version of R and RStudio, and Rtools? |
| NA | I am using the R 4.0.2 |
| NA | Cool. What about Rtools? |
| NA | Are you working with Windows? |
| NA | Version 1.3.1073, yes, windows |
| NA | RStudio Version 1.3.1073 |
| NA | What does your .libsPath say? |
| NA | I use the default path |
| NA | Sounds like you don't have a C++ compiler. Did you follow nimble installation instructions? |
| NA | Are you the administrator of your computer? On windows sometimes writing is lock and cause troubles with installation. Otherwise shutdown RStudio, restart it and check installed librairies |
| NA | Yes, I am the administrator of my computer, maybe the problem is the C++ compiler, I try to see if I can fix it by following the nimble installation instructions... |
| NA | thanks all!! |
| NA | Let us know how it goes. Thank you all for helping. |
| NA | Thanks all for the information, it is done!! I run the printErrors() and it says the problem is make.exe, and I reinstall rtools and put rtools on the PATH, then it is perfectly fixed!! |
| NA | defining model... building model... setting data and initial values... running calculate on model (any error reports that follow may simply reflect missing values in model variables) ... checking model sizes and dimensions... checking model calculations... model building finished. compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details. compilation finished. running chain 1... |-------------|-------------|-------------|-------------| |-------------------------------------------------------| running chain 2... |-------------|-------------|-------------|-------------| |-------------------------------------------------------| |
| NA | Awesome, glad to hear it! |
| Naive question, but there is a way to implement in the models the time effect directly on recapture data ? How we can deal with missing year or a decrease in sampling effort across year for example ? | If you want a direct time-effect (i.e. a trend), you would just have a slope parameter multiplied by the time index, e.g. beta*t |
| NA | For dealing with missing years: you could include that by using a "dummy covariate" that has a value of 1 whenever there was sampling, and 0 when there was none. In your model, you then multiply your recapture probability with that covariate. So your p will become 0 in all years when there was no sampling, and will be estimated "as normal" when there was sampling. |
| NA | Analogously, you could include a dummy covariate with values of 1 and 2, and use it as a nested index. it works similarly if you have multiple categorical levels of sampling effort (i.e. "none", "low", "high"). |
| NA | Finally, continuous measures of sampling effort can just be included as continuous covariates :wink: |
| NA | Many options! |
| NA | Hi , sorry for coming back so late after the course but I have recently used the dummy variable option in my models and one of the reviewer is asking for published papers describing this method. I know I have read it on various occasion before but I cannot find an example. Would you have such papers in mind? Thanks a lot |
| NA | great to hear that this was useful to you! I'm working with a model where I use different (dummy) variables to capture variation in sampling effort, but it's not published yet. I am also sure that there are many other published analyses that use this approach. |
| NA | That being said: using a dummy variable is not really a "method", it's just one way to parameterize a model including variation in sampling effort. If you look at a cake recipe, it will say "add sugar", but it usually will not specify if the sugar is to be added by spoon, by hand, or in a different way. Because at the end of the day, that does not matter :wink: |
| NA | You still have to answer to that reviewer though. My suggestion would be to explain that it's a common approach (maybe plus a reference to one or too examples, and maybe the butterfly book - I do think they introduce dummy variables there too, even if not exactly in this context, check!), then illustrate it to the reviewer with an example. Something like: The equation for recapture probability is p.eff[t] = dummy[t]*p[t] In years when there was sampling (recaptures of individuals), dummy[t] = 1, so p.eff[t] = 1*p[t] = p[t] In years when there was no sampling, dummy[t] = 0, so |
| NA | p.eff[t] = 0*p[t] = 0 |
| NA | I also ping . Maybe he knows something good to cite here! Olivier, please check above :arrow_up: |
| NA | As it is often the case, there is something in the Gentle introduction to Mark |
| could you send the ref on the 27 european passerines? thanks! | |
| NA | thanks! |
| What happens if your prior is wrong. Say for the dipper you use a prior for non-flood years but then estimate survival during a flood year? Is that an issue? | Depends on how much "information" is in your data and what type of prior you use. If you have plenty of data and you use a prior without limits, then the information in your data will "outweigh" your prior info, and your posterior will be adequate. |
| NA | If you have plenty of data and you use a prior with limits (i.e. uniform prior), the same as above will happen but if the data wants to "push" outside the prior, and you will get a posterior which looks "stuck" to one of the boundaries of the prior. |
| NA | great thanks! |
| NA | If you do not have much data and your prior is "wrong", then your parameter will be estimated close to the prior and the "mismatch" with reality will propagate as bias into the other parameters in your model. |
| NA | If you work with informative priors, it's always good to do a prior sensitivity analysis, i.e. run the model with a bunch of different priors and check if the estimates / conclusions are affected :slightly_smiling_face: |
| NA | :+1: |
| NA | Does the WAIC give this information somehow? Or is there another way than simulations to estimates the "weigth" of data over prior estimations (or the other way around...)? |
| NA | You have to do a prior sensitivity analysis. Prior and/or posterior predictive checks can also help. WAIC does not help here. Working with Bayesian models requires a good deal of "common sense", and there are no easy yes-no fixes for things like prior sensitivity :wink: |
| NA | Learning to use simulation tools to investigate things like prior sensitivity is a very handy skill to have, in a lot of different setting. It takes some effort to get going, but it will be worth your time :slightly_smiling_face: |
| NA | "Making mindless binary decisions based on one metric or another is always a bad idea" - Yourself. One of my take home message :wink: |
| NA | Actually working on a commentary paper on this with colleagues :stuck_out_tongue: |
| NA | I look forward to seeing it :slightly_smiling_face: |
| NA | Sounds great! |
| NA | Hi all, interesting discussions! Just wondering - if there are quite a lot of data, and the data kind of outweighting the prior info, does it mean that the results will be pretty close to likelihood results? |
| NA | yes, that should be the case. |
| NA | Thanks Chloe, if this is the case, would there be a concern for losing advantages of Bayesian inference, as the results won’t be too different from likelihood methods? |
| I don't understand why the prior and posterior should not be too confounded... Wouldn't it help the model to estimate better survival if the prior of survival, for example, is already well estimated? I thought that was the case, following what you explained on the european dipper. Thanks in advance for your answer ! :slightly_smiling_face: | You are correct that a good prior will help the model make the estimate. BUT: how do you know your prior is really, truly good? For that reason, when we work with informative priors, we always do a prior sensitivity analysis, i.e. we run the model with a bunch of different informative priors. What we want to see is that the posterior is always similar independent of what prior information we provide. |
| NA | So what we need to ensure is that prior and posterior are not too "confounded" ACROSS a variety of models using different prior information. |
| NA | Thank you Chloé for your answer ! So if I understand well, you are testing a bunch of different priors, and for each of these priors, you check whether the posterior is confounded with the prior, it that correct? And if for each prior, posteriors are confounded with their corresponding priors, that we mean that the posterior are too dependent on the prior values? |
| NA | Yes, exactly! In addition, we also want to check to what degree prior choice affects the estimates of other parameter in our model and - most importantly - our conclusions. |
| NA | Okay, I see !! Thank you very much again !! :smile: |
| is there a rule of thumb on the minimal ratio of marked animals in a given population you should have to have an accurate enough estimate of survival? | I'm afraid we don't have a rule of thumb, I guess the answer is problem-specific, and my advice would be to use simulations to assess the precision of your survival estimate given the biology of your species, the effort you're ready to invest, etc... |
| NA | I will try and provide a few references in a minute or so. |
| NA | Let's start with |
| NA | And this one |
| NA | aaaaaaaand |
| NA | Last one |
| NA | super! last one is probably the closest from my issues |
| NA | Excellent references - my biblio is growing exponentially during this workshop! |
| NA | I love biblio, don't get me started :hugging_face: |
| NA | We used a multistate CJS model to estimate annual marine return rates of Atlantic salmon ( |
| NA | Very nice paper, working w the :fr::baguette_bread: crowd I see :wink: |
| NA | Indeed - Marie and Etienne have been my mentors. They talk very highly of you and your team too :slightly_smiling_face: |
| NA | Working with them would be another good reason to move to Brittany :heart_eyes: |
| First of all, great workshop, thanks! I understand the issue of Bayesian vs. frequentist modelling is not the topic here, but I was wondering: let's say were are in a CMR applied context, unconcerned about using elegant, "fashionable" methods and maximising the probablity of producing a high-impact paper, why would we opt for a Bayesian approach? Particularly if we don't have informative priors? | Let's see how this could work: P(high-impact paper | bayesian analysis) = P(bayesian analysis | high-impact paper) * P(high-impact paper) / P(Bayesian analysis) |
| NA | Jokes aside the Bayesian framework gives you way more flexibility to expand models and add additional needed complexities compared to the frequentists one where the algorithm often cannot deal with the complexities. |
| NA | Thank you for your kind words . Not much to add to answer. On a technical side, i) we often get boundary estimates with max lik estimation (survival estimate hits 1 or 0 for example), and using weakly informative or even vague prior fixes that issue, and ii) random effects are much easier to accommodate in the Bayesian/MCMC framework. Others will probably have something to add |
| NA | :joy: provocative statements always get a response ! Thanks for the replies. |
| NA | Yes indeed :slightly_smiling_face: |
| NA | Yep - definitely using it for the random effetcs part. Also when you want to write your own model (for example if you have a very particular study design) rather than using a pre-definite model, it’s relatively easy to do in a bayesien framework which is very flexible. |
| NA | Beyond the scope of this workshop, but in a Bayesian modelling context it is easier to integrate multiple datasets. Also, it might not increase your chances of a high impact paper, but when you use Bayesian inference you get to use concise, clear, and intuitive language to discuss your results because everything is interpreted as a probability. |
| NA | Yes, I've been sold the benefits of integrated pop modelling... :+1: |
| Can the Eps[t] part of the phi parameter definition be interpreted as the background survival rates, independent of the covariate effect on survival rate (beta2)? | The eps[t] are the levels of a time-random effect. They can be interpreted as residual variation in addition to covariates. So, if you have survival affected by temperature and two model like this: *Model A* logit(phi[t]) <- beta[1] + beta[2]*temperature[t] + eps[t] |
| NA | *Model B* logit(phi[t]) <- beta + eps[t] |
| NA | Then the (variance of) eps[t] will be larger in Model B, because they also include the temperature effect. In Model A, we will have lower "residual variation" (= variance of eps[t]) because we explicitly account for the temperature effect. Does this make sense? |
| NA | Yes, great thanks. I think I was considering background = residual in the case of the example, but it all depends on what else is defined in the model, so it's not quite the same. |
| I am not sure why dcat is better than multinomial. can you explain? | dcat is a multinomial, just that the number of trials is 1 (just like the bernouilli is a binomial with 1 trial). |
| NA | ok, one roll of die? so it might be appropriate for multinomial outcome with one trial? |
| NA | one roll of a die that has different probabilities of landing on each of it's sides (technically, they could be equal probabilities too I guess)! |
| NA | ok, thanks. does it work also for jags or just nimble? |
| NA | Works for Jags as well. |
| NA | Missed Sarah's lecture on priors m oment matching. when is it best to use informative vs uninformative priors vs. vague priors or with very large variances (as seems to be a common practice among many bayesians) |
| NA | can you ask this question in a new thread? |
| NA | yep |
| NA | thanks! |
| First of all, thank you very much for this fantastic workshop! I wonder what are the consequences of assuming that all individuals are alive at any state (zinits)?? Thank you! | The consequences of the initial values for latent states are related to speed to convergence. The zinits tells the MCMC where to start. The closer to the "truth" it starts, the faster it will converge to that truth. Of course, if you work with a short-lived species, it may not be realistic to set "alive" as the initial state 20 years later. So you can also change that at some point (as long as zinits between two observations are always 1). What we do more commonly, actually, is that we only analyse the capture histories of an individual for a "biologically realistic" number of years, to save computation time :wink: |
| NA | Thank you very much! So for example if your monitoring period is 12 years but the lifespan of the species is 4 years you can assume that each individual would be alive during 4 years from first capture and for the rest of years you can assume they are dead? |
| NA | Great, you have just answered my question during the talk :slightly_smiling_face: It can be solved by setting a "last" capture. Thank you! |
| NA | Great that that's clear now :slightly_smiling_face: |
| Missed Sarah's lecture on priors m oment matching. when is it best to use informative vs uninformative priors vs. vague priors or with very large variances (as seems to be a common practice among many bayesians) | I'd say, if you have prior information, just use it, but carry out a sensitivity analysis to better understand how much posterior inference is driven by your prior. |
| NA | sensitivity analysis by changing priors? just the moments, distributions or both? |
| NA | Wrt large variances on the regression parameters (the betas), I tend to avoid them cause you risk getting very informative priors when you back transform to get your prior on survival (bathtub shape). That's why we use dnorm(0, sd=1.5) in the lecture/live demos. My advice here is to simulate data from your priors on betas, then visualise the prior you get for survival. |
| NA | By changing the priors yes. |
| NA | We need an example of a sensitivity analysis in the workshop I think :slightly_smiling_face: |
| NA | ok, thanks - that would be helpful! I have a related question about situation when posterior distribution is weirdly and perhaps unreasonably bimodal. should i ask that question in a different thread? |
| NA | Yes, please. |
| NA | Just to illustrate the bathtub thing |
| We have a situation where we are trying to determine a kind of threshold when density-dependence begins affecting say survival (e.e.g, via intra-specific aggression) in a large carnivore. we were advised by a Bayesian Guru to use flat prior [unif (10, 200)]. Now we get two distinct peaks in the posterior distribution, one around 50 and another around 160. Any suggestions on how to find the right priors that will get us closer to the "truth"? | To find the "truth", you need auxiliary data. Either in terms of expert opinion, literature reports (on similar species), and/or hard data from your system. |
| NA | A bit tricky to answer without knowing more about the model. Did you change your priors and still found two peaks? Do you use some king of mixture model in which case label switching might be an issue? |
| NA | If none of these are available...I'd suggest going with the flat prior and reporting the results as they are, i.e. there may be two solutions. Use the discussion to argue about evidence for/against one or the other. |
| NA | Have you tried simulating data from your model? Could be a way to understand what's going on. |
| NA | Olivier, should I write to separately about this if you think it is not relevant for this group/thread? |
| NA | Sure, shoot :slightly_smiling_face: |
| A bit off topic: Do we have reasons to believe NIMBLE will stay here for long? Or will it be another JAGS/BUGS? | is in his prime age so I'd say yes, we may be confident that it'll be here for long :wink: |
| <@U01UMBR0K38> | we should make a list of all the excellent advices in pictures !! |
| NA | I love this! |
| Are dirichlet priors the priors of choice for multistate param? There is no such things as multi-logit? | Can you ask this question in the questions-6 channel please. It'll be easier for us to organise the FAQ afterwards. Thanks. |
| NA | done - thanks |
| NA | Thank you! |
| I am not sure where I should do this question. Is it possible to perform some model similar to the POPAN approach with Nimble to estimate the (super-)population size? | I am not familiar with that particular approach, but if it's a type of open population CMR, then yes, I am sure it can be done. I like to think that if you can understand a model, you can write and run it in nimble :wink: |
| NA | Maybe knows of a Bayesian application of this particular one? |
| NA | There is Jags code in Kéry and Schaub's book, chapter 10 |
| NA | Thank you so much!:hugging_face: |
| There are identifiability issues with the last survival and recapture probabilities. If we have a situation where we definitely think these two are time-dependent, can we still fit a model where we specify them as such, but then do not report the last survival and recapture probs? (I am aware I would be "dropping" one year of data) | Yes, its fine as long as you mention somewhere that you dont report these estimates because the corresp parameters are redundant. A trick to make this model identifiable would be to use temporal random effect on survival and/or detection probs. |
| me again...:zany_face: In live demo 5 you say "The easiest way to cope with time-varying individual covariate is to discretize and treat levels of the covariates as states. More in the next live demo." and indeed in live demo 6 there are examples of this, where one tries to estimate probabilities of transition between states. But what about age (as classes/states)? The transitions are not probabilistic (is this a word?, well ok they are in a way) and are known (linked to duration between survey). I think I would vaguely remember how to do this in mark (because I did 20 years ago, so how could have I forgotten?), but a nimble example would be great - any example to share? (you see, you provide amazing workshop and material, and people want even more). | You can certainly treat age as a state, but it is easier to have it as a covariate (unless you have uncertainty in age assignation, in which case it is useful, see e.g. |
| NA | well yeah, that's what I would do, I was getting confused by the sentence. But thanks for the examples, that indeed will be useful. :pray: |
| Hello, I don't know if the SLACK is still active but I'm trying... I use a data for which no capture has been done some years. I would like to run a phi(.)p(t) model by fixing the p to 0 for the years without capture. I tried to add a covariate (0 for years without capture and 1 with) but I do something wrong because after multiple attempts it still doesn't work. In fact my code do not really consider time effect but only consider a p for years "without capture" and a p for years "with captures". I'm not sure it's really clear, sorry ... Does anyone have a solution? My nimblecode is : survey<-c(1,1,1,1,1,0,1,1,0,1,1) ##covariate hmm.phipt <- nimbleCode({ delta[1] <- 1 delta[2] <- 0 for (t in 1:(T-1)) { logit(p[t]) <- beta[1] + beta[2] * survey[t] omega[1,1,t]<- 1-p[t] omega[1,2,t] <- p[t] omega[2,1,t] <- 1 omega[2,2,t] <- 0 } phi ~ dunif(0, 1) gamma[1,1] <- phi gamma[1,2] <- 1 - phi gamma[2,1] <- 0 gamma[2,2] <- 1 beta[1] ~ dnorm(0, 1.5) beta[2] ~ dnorm(0, 1.5) # likelihood for (i in 1:N){ z[i,first[i]] ~ dcat(delta[1:2]) for (j in (first[i]+1):T){ z[i,j] ~ dcat(gamma[z[i,j-1], 1:2]) y[i,j] ~ dcat(omega[z[i,j], 1:2, j-1]) } } }) |
The slack is still active yes :wink: and I might have a clue. When `survey[t] = 0` for some `t`, `logit(p[t]) = beta[1]` and the prior on `beta[1]` being `dnorm(0,1.5)`, this won't give you `p[t] = 0`. Nested indexing might help. Try something like `p[t] <- beta[survey[t]]` with `survey[t] = 1` when `t` is a year with survey, and `survey[t] = 2` when `t` is a survey with no survey. You'll need to redefine your covariate survey, 2 - survey where survey is your current covariate will do the job. Then assign `beta[1] ~ dunif(0,1)` and `beta[2] <- 0`. Should work. If not, send me a fully reproducible example offline, and I'll try to debug it, then we'll report here in case it's useful to others. |
| NA | You could also use `p[t] <- 1 / (1 + exp(-beta) * equals(survey[t], 1)` with `beta ~ dnorm(0,sd=1.5)`. When `t` is a year with survey, then `equals(survey[t], 1)` is 1 and `p[t]` is estimated. When `t` is a year without survey, then `equals(survey[t], 1)` is 0 and `p[t]` is 0 too. Note that you could use `step(1 - survey[t])` instead of `equals(survey[t], 1)`, I think. This solution is probably more practical when it comes to have other covariates on detection, with e.g. `p[t] <- 1 / (1 + exp(-(beta[1] + beta[2] * x[t])) * equals(survey[t], 1)`. |
| NA | Thank you ! But when I do the both options a have the same problem that for my first code : I have 10 p bu 8 are the same ( for the years with surveys) and 2 are equal to zero (for the years without surveys)... I tried others solutions |
| NA | Almost there :wink: My bad, I forgot to include time-dependent detection probabilities. With the second option, try `p[t] <- beta[t] * equals(survey[t], 1)` with a loop for `beta[t] ~ dunif(0,1)`. If not, email me your code and data, I'm happy to help. |
| NA | Thank you very much for your help! After some attempts, it perfectly works with "survey<-c(1,1,1,1,0,1,1,0,1,1) hmm.phip.t <- nimbleCode({ delta[1] <- 1 delta[2] <- 0 for (t in 1:(10) { beta[t]~dunif(0,1) p[t]<- (1/(1+exp(-beta[t])))*survey[t] omega[1,1,t] <- 1 - p[t] omega[1,2,t] <- p[t] omega[2,1,t] <- 1 omega[2,2,t] <- 0 } etc..." |
| NA | Ha of course, just use survey as an indicator :wink: Instead of the uniform prior, i would rather use a normal prior. The parameter is not necessarily positive and lower than 1. You may use prior predictive checks to see what prior you imply on pt. E.g. `hist(1/(1+exp(-(rnorm(1000,0,sd=1.5)))))`. |
| NA | Thanks! Last question (I hope). I am not sure about the goal and the interpretation of the histogram (I am begginer on bayesien method sorry) |
| NA | No pb. It tells you that a normal prior with mean = 0 and sd = 1.5 on beta induces a uniform prior between 0 and 1 for `p = 1/(1+exp(-beta))`. |
| NA | Ah yes, perfect ! Thank you very much for your help |
| <@U01UMBR0K38> might have examples of dead-recovery data analyses in Nimble (trouts?). In any case, the translation Jags -> Nimble shouldn't be too difficult, check out some advices here <https://r-nimble.org/quick-guide-for-converting-from-jags-or-bugs-to-nimble>. Happy to help if you bump into problems. As for IPMs, <@U01UQ29DG3C> in her recent paper <https://esajournals.onlinelibrary.wiley.com/doi/abs/10.1002/ecm.1459> has made her Nimble code available <https://github.com/maudqueroue/MultispeciesIPM_SkuaPetrel> :wink: See also <https://esajournals.onlinelibrary.wiley.com/doi/abs/10.1002/ecm.1364> and probably other resources. Hope this helps. | Thanks, I'll check them out - having been influenced by this workshop my plan was to switch to Nimble after finishing this work as it was already well advanced having learnt jags over the last few months but it looks like I need to do it during it! Also btw great EURING talks and and potentially will inform some of the stuff I'm doing! |