Issue in estimation of parameters in a two-sex integrated population model

[Jessica Cachelou]

Dear Nimble users,

My question is quite specific, so a lot of thanks in advance for people who will take time to give an answer. 

I am working on the Bayesian Integrated Population Model developed in Gamelon, Nater et al. (2021) (Ecography, https://doi.org/10.1111/ecog.05738). This IPM has been built for wild boar females and is structured in body mass classes: small, medium and large. My aim is to explore the role of males on wild boar population dynamics.  I have thus included males in the IPM code (please, find it in the attached file). 

 

I first used the same datasets (e.g. CMR data, hunting data) for both males and females. I expected to obtain the same estimates of mortality probabilities and population sizes for the two sexes. We obtained estimates for females (population size, stage-specific mortalities) perfectly comparable to the ones published in the paper (Gamelon et al. 2021) last year. But for males, estimates are different from females. The temporal variation in all parameters seems the same for both sexes, but the mean estimates are not correct for males. For example, population size is underestimated for males (especially for medium and even more for large individuals), moreover, the probability to grow from medium to large (gPm[2]) and the probability to grow from small to large (gSLm) were lower in males than in females. In the same way, the probability of hunting mortality (sHm) was higher for medium and large males than for medium and large females. The probability of natural mortality (sNm) was higher for small and medium males than for females. However, population size for offspring and small are similar for males and females.

 

We tested a lot of things to understand why we obtained such a discrepancy between males and females:

-          Change seed: we used 100 instead of 0;

-          One matrix 6*6 for the two sexes: instead of having two different growth matrices (dim = 3*3), one for males and one for females, we wrote one single matrix of dimensions 6*6;

-          A lognormal distribution for population size, instead of the truncated normal distribution;

-          A sex-specific postnatal survival, so we added a parameter for postnatal survival of males, instead of having one common parameter for the two sexes; 

-          We tested the IPM with males only, everything worked perfectly; we tested the IPM for females only; everything worked. But when both sexes were included, males’ estimates became wrong; 

 

After talking with Chloé Nater (many thanks for your help), we tested additional changes:

• The bias is still present when running multiple chains and it is not a matter of convergence to alternative solutions;

• Running the CMRR model separately gives parameters that are identical between males and females, the issue is not in the CMRR data/likelihood;

• Unlinking the sexes (i.e. treating M and F as two independent female populations) also does not eliminate the discrepancies;  it is not the link of the two populations through reproduction that causes the problem

 

We do not have any ideas of the issue; however, we suppose we are doing something wrong with NIMBLE. If you have a bit of time, I would be more than grateful if you could have a quick look at our code (in attachment). There is probably an obvious mistake that we cannot see.

However, I can not publish the data, so if you need to run the code, ask me and I’ll give you a part of this dataset and the code for the initial values.

Many thanks in advance for your valuable help.

Cheers,

Jessica Cachelou 

[rolek...@peregrinefund.org]

Hi Jessica,

I recently had similar issues with an IPM. I resolved this by using posterior predictive checks to check the goodness-of-fit of my count and fecundity models and ensure they were providing reasonable estimates. I needed to switch from normal and log-normal distributions, to Poisson, zero-inflated Poisson, or negative binomial. To speed things up, I checked the fit of count distributions outside of the IPM first, and this translated well within the IPM. 

So you might give this a try,. You could check out the code from my recent pub for an example in JAGS code. Or more complete treatment is contained in the book Integrated Population Models by Kery and Schaub Chapter 7 (code freely downloadable on the website). If you have trouble finding examples for all the different distributions, just email me and I can pass along the code that was not published with the IPM. Lastly, if you hit more road blocks the hmecology listserv frequently has the authors of the IPM book on it, and I'm sure you could ask them there. 

Cheers,

Brian Rolek

[Perry de Valpine]

Thanks Brian for posting suggestions.  Those are great to see.

Jessica, thanks for posting the question.  This is a fairly involved model and a question that may be hard to pin down.  To follow up on your offer, yes it would be helpful if you can make a fully reproducible example, with a subset of data or simulated data or whatever you can provide.  It would also help if you can try in any way to make it a minimal reproducible example, if there are any parts of the model that can be cut out or simplified while still displaying the problem behavior you are trying to understand.

Are you also providing a complete set of initial values for each sex and setting those to be identical?

Is there anywhere at all in the model definition where the two sexes share information?

Perry

 

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/17114f7f-e8a2-439d-8302-1033d3a70814n%40googlegroups.com.

[Jessica Cachelou]

Thanks Brian and Perry for your suggestions. 

Firstly, I send in attachment the R code for the IPM and for the initial values, and a part of my dataset. 

Brian, I have already used posterior predictive checks to check the goodness-of-fit of each part of my model, and everything was estimated reasonably. So, I do not think it is that. Moreover, I also switched from normal distribution to log-normal distribution for the population count and nothing changed (just it was faster). But I can test the switch to Poisson / zero-inflated Poisson / negative binomial !

Thank you so much for all the references ! I'll read them carefully. 

Perry, to answer your questions, I am providing a complete set of initial values (code in attachment) that set identical for both sex. 

Moreover, there is no sharing information between two sexes. We just have a postnatal survival that is common to females and males, but, we tested a sex-specific postnatal survival (in order to totally separate females and males) and nothing changed. 

Cheers,

 Jessica 

You received this message because you are subscribed to a topic in the Google Groups "nimble-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/nimble-users/olonWSd1Ayc/unsubscribe.
To unsubscribe from this group and all its topics, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/CABeTmqgG1MA2w372xE89mk2qdC33TZNN8ZMZvk1rE7cvHedKkA%40mail.gmail.com.

[Perry de Valpine]

Hi Jessica,

A clarification question: It looks to me like the dynamics of the two sexes are not entirely separated in the model but rather intersect at reproduction, where only females produce offspring that are then split into female and male offspring.  This is of course realistic.  But in terms of the statistical structure of the model, wouldn't it be a reason that parameters for males and females could differ?  For example, the number of females could be pulled up or down by the data (e.g. via sex-specific survival or growth rates) so that the reproduction step helps the model match the data, and this could drive differences in male vs female parameters?  Since I think you're trying to determine if the issue is statistical vs a purely software issue, would you want to make the males and females completely not share information in some artificial way?  I thought that's what you were saying but I am not clear that is what is happening in the model due to information-sharing via reproduction.  It is possible I'm not following your question perfectly.

Perry

[Jessica Cachelou]

Hi Perry,

Yeah, sorry I did not understand what you meant in the previous mail ! Indeed, females and males are linked in the production of the offspring. 

But, I agree with what you said, and it was one of our suggestions, so we also tested an IPM where females and males were treated as two independent populations. We used the following code where s is the sex (1 for female and 2 for male), t the time and z the body-size class:

# Unlinked model (treats m and f as two independent f populations)
   for(s in 1:2){
     for(t in 1:Tmax){
       for(z in 1:Z){
        
         # Number of offpring produced by parents in size class z
         Off_tot[z,t,s] ~ dpois(marN[z,t,s]*pB[z,t]*nF[z,t]*0.5)    
       }
       
       # Number of produced offspring recruiting into the population
       YOY[t,s] ~ dbin(s0[t,s], sum(Off_tot[1:Z,t,s]))
     }
   }

And we did not see any difference in the results. That is why we supposed a potentially software issue. 

Jessica 

[Perry de Valpine]

Hi Jessica,

Thanks.  I replaced the model code labeled as "# Linked model" with the "# Unlinked model" you've provided and tried to run your code.

Your code doesn't run for me.  I get:

Defining model
Error in addMissingIndexingRecurse(code[[i]], dimensionsList) :
  inconsistent dimensionality provided for node 'Off_tot'

I also tried the build the model directly, outside of nimbleMCMC:

m <- nimbleModel(WB.code, constants = WB.constants, data = WB.data, inits = Inits[[1]])

and got the same error.

So I am wondering if you can check that what you sent as a reproducible example does work on your system?

I inspected possible mismatches in size and/or dimensionality between the model definition and the Inits as follows:

m <- nimbleModel(WB.code, constants = WB.constants, data = WB.data, calculate = FALSE)
for(v in names(Inits[[1]])) {
  dim_or_length <- function(x) if(!is.null(dim(x))) dim(x) else length(x)
  if(v %in% m$getVarNames())
    cat(v, ':', dim_or_length(Inits[[1]][[v]]), ';', dim_or_length(m[[v]]),'\n')
}

The output shows first the dims in the Inits[[1]] and next the dims of the variable in the model.  I mark the ones that don't match:

marN : 3 12 2 ; 3 12 2
initN : 3 2 ; 3 2
marN_g : 3 3 11 2 ; 3 3 11 2
marN1_plus : 12 2 ; 11 2
octN_g : 3 11 2 ; 3 11 2
H : 3 12 2 ; 3 12 2
Off_tot : 3 11 ; 3 11 2. ****** Missing a 3rd dimension in Inits
YOY : 11 2 ; 11 2
Mu.ll : 2 ; 2
sigma.ll : 2 ; 2
epsilon.ll : 11 2 ; 11 2
Mu.pp : 3 2 ; 3 2
sigma.pp : 2 ; 2
epsilon.pp : 11 2 ; 11 2
Mu.mH : 3 2 ; 3 2
sigma.mH : 2 ; 2
epsilon.mH : 11 2 ; 11 2
Mu.pB : 3 3 ; 3 3
sigma.pB : 1 ; 1
epsilon.pB : 11 ; 11
Mu.nF : 3 ; 3
sigma.nF : 1 ; 1
epsilon.nF : 11 ; 11
Mu.m0 : 2 ; 2
sigma.m0 : 2 ; 2
epsilon.m0 : 11 2 ; 11 2
Mu.mN : 3 2 ; 3 2
sigma.mN : 2 ; 2
epsilon.mN : 11 2 ; 11 2
Mu.gP : 2 2 ; 2 2
sigma.gP : 2 2 ; 2 2
epsilon.gP : 2 12 2 ; 2 11 2 ************* different sizes
Mu.gSL : 2 ; 2
sigma.gSL : 2 ; 2
epsilon.gSL : 12 2 ; 11 2 ********** different sizes
x : 870 11 ; 870 11 

I'm assuming you had it running, and if you send updated code to fix these issues I will keep trying to reproduce the problem you're encountering.

Perry

[Jessica Cachelou]

Hi Perry,

Sorry, when I sent you the part of the code of  # Unliked model, I forget to sent you these lines to reformat the initial values :

# Reformat initial values for unlinked model
for(i in 1:nc){
  Inits[[i]]$Off_tot <- abind::abind(Inits[[i]]$Off_tot/2, Inits[[i]]$Off_tot/2, along = 3)
  Inits[[i]]$Off <- NULL

}

You need to add these lines after the "Inits", line 728. You can find in attachment the updated version of code if you prefer.  

However, for the dimensions of epsilon.gP and epsilon.gSL, I have the same dimension for the two chains in my system... So I don't understand why it's not right for you.  

Jessica