[Note] Infinite values were detected in model variables: logProb

[emilygs...@gmail.com]

Hi everyone,

I have been having some issues with a survival model I am trying to run in nimble. 

Specifically, I am finding that some of my variables are producing logProb of infinity and whole model likelihood calculates at -Inf. Does anyone have any ideas why this can happen? 

I have attached my code and the data I am using here incase anyone can give any pointers. I have tried to code age by having:

juvenile = 1, adult = 2, dead = 3. If there are better ways please do suggest some. 

Any help is greatly appreciated. 

Emily

[Daniel Turek]

Emily, thanks for your use and interest in NIMBLE, and also your patience with your question.

There are a few things with your code that needed addressing, which were preventing your model from working correctly.  Please take a careful look at the descriptions below, and the updated code.  Assuming these changes are all consistent with your *desired* model structure, then this should get you going.

I'll also point out that the dcat-dcat structure you have with the latent "surv_state" variable and the "surv_obs" data might benefit from the marginalizing distributions available in the nimbleEcology package - specifically in this case, I believe the dHMM distribution.  But before you consider that, please spend some time reading and thinking about the changes and comments below.

(0. I commented out the line "input_data <- output_data %>% filter(Year < 10)", which I don't think was intended to be in this script, where data is loaded from the file you provided)

1. Defined a new constant "first", giving the temporal location of the first-non-3 (alive, juvenile or adult) observation of each individual.  This is a commonly needed value in mark-recapture models of repeated observations of individuals, telling the model when to first use the dcat(initial_survival[1:3]) distribution, and subsequently, when to use the transition matrix.

2. Providing the "first" variable in the "constants" list provided for the model.

3. Changing the initial distribution for the first observation to begin at surv_state[i, first[i]], and the following loop to begin at (first[i]+1) .. (see code file)

4. Changing the initial values for "surv_site", specifically "surv_state_init" to be plausible.  Equal to 1 or 2 on the first observation (the value which is consistent with capture_history), then 2 thereafter.  Previously, your "surv_state_init" initial values were all 1's, for all individuals and all time periods, which given your observation and transition matrices, was generally impossible values given the capture histories.

5. Changing the initial values provided for "offspring_state", rather than to be randomly Poisson-distributed, which often times gave values that were 0 (and thus could not be used for the Poisson rate parameter for non-zero data) to instead use the values of "offspring_obs".  If you want something random here, you could still use something like what you had, but it has to be positive, so perhaps rpois(length(offspring_obs),offspring_obs)+1, adding the "+1" at the end.

These represent fundamental changes to your model structure - in particular, at what time periods the model starts to consider the latent state values and initial probabilities, and what time periods in the capture histories are considered to be observed data for the likelihood.  If I've made incorrect assumptions about your intended model, then you should not blindly accept these changes I've made - it's entirely possible that I've misunderstood your intended model.

Modified code is below, along with ## CHANGE to indicate where substantive changes have been made.  I hope this helps.

Cheers,

Daniel

--
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/67520f02-8e99-469b-be0b-b5d55289743cn%40googlegroups.com.

[emilygs...@gmail.com]

Hi Daniel,

Thank you so much for your response and for taking so much time to look at my issue. 

I also kept working on this while waiting for a response and luckily have some of the same fixes. Specifically, I did have the epiphany that a major issue was I was giving impossible values for the two state variables (offspring_state and surv_state). I had fixed the offspring_state the same way you suggested but had issues with the surv_state as I couldn't work out how to deal with transitions from 'not yet alive' to 'alive'. The addition of the 'first' constant really fixes this so thank you for that suggestion. 

So, I think the assumptions you made about the model were correct, at least they follow the exact line of fixing I was trying too. 

Following on from that, I had a couple of further questions:

1. In my attempts to fix the model myself, I considered re-coding the 'observations' matrix so that unobserved = 1, observed as juvenile = 2 and observed as adult = 3. I felt this coding was closer to the example in a hmm tutorial here https://oliviergimenez.github.io/banana-book/hmmcapturerecapture.html. However, I note you did not need to do that. Does it make any difference whether I use 1 or 3 for 'unobserved' as long as the probabilities and link to the state values are correct? I'm feeling it may not. I have pasted my revised code below. Not sure if I was just overcomplicating for no gain.

2. I am interested in the marginalising distributions in nimbleEcology but was not sure if they would work for a multiple age structure like this model and with the 'first' constant? If they do, please let me know. I assume I would just use the matrices of observations, transitions, initial_surv etc the same as in any binary example? Would the coding choice mentioned in 1. have any impact here? 

Thank you again, I really appreciate being able to talk this through. 

Emily

MY ATTEMPT AT NEW OBSERVATIONS MATRIX: (obviously I also changed the observed data to be 1 = dead/not seen, 2 = juvenile, 3 = adult but kept the initial state values at 1 = juvenile, 2 = adult, 3 = dead)

# matrix of transitions from juv, to adult, to dead (STATE)
# the probabilities should be the probability of being 1, 2 or 3
# STATES = alive juv 1, alive adult 2, dead 3

transition[1, 1] <- 0 # Pr(juv alive t -> juv alive t+1)
transition[1, 2] <- mean_phi_juv # Pr(juv alive t -> adult alive t+1)
transition[1, 3] <- 1-mean_phi_juv # Pr(juv alive t -> dead t+1)
transition[2, 1] <- 0 # Pr(adult alive t -> juv alive t+1)
transition[2, 2] <- mean_phi_adult # Pr(adult alive t -> adult alive t+1)
transition[2, 3] <- 1 - mean_phi_adult # Pr(adult alive t -> dead t+1)
transition[3, 1] <- 0 # Pr(dead t -> juvenile alive t+1)
transition[3, 2] <- 0 # Pr(dead t -> adult alive t+1)
transition[3, 3] <- 1 # Pr(dead t -> dead t+1)

# observation matrix (captures recapture probability)
# OBS = not detected 1, alive juv 2, alive adult 3
# row 1 = alive juv
# row 2 = alive adult
# row 3 = dead

# column 1 = not observed, column 2 = observed juv, column 3 observed adult

observations[1, 1] <- 1 - mean_p # Pr(alive juv t -> non-detected t)
observations[1, 2] <- mean_p # Pr(alive juv t -> detected juv t)
observations[1, 3] <- 0 # Pr(alive juv t -> detected adult t)
observations[2, 1] <- 1 - mean_p # Pr(alive adult t -> not detected t)
observations[2, 2] <- 0 # Pr(alive adult t -> detected juv t)
observations[2, 3] <- mean_p # Pr(alive adult t -> detected adult t)
observations[3, 1] <- 1 # Pr(dead t -> not detected t)
observations[3, 2] <- 0 # Pr(dead t -> detected juv t)
observations[3, 3] <- 0 # Pr(dead t -> detected adult t)

[Daniel Turek]

No problem, Emily.  Apologies that my response is brief, but:

1. The values you use for coding the states makes no difference whatsoever, aside from perhaps lending intuition (so long as the transition matrix indicies, etc, are all consistent).

2. Yes, the nimbleEcology distributions would work fine for this.  Each full capture history surv_obs[i, first[i]:occasions] would follow such a custom distribution from nimbleEcology.  Please carefully read the documentation and examples for the distributions dHMM, dHMMo, dDHMM, and dDHMMo, all from nimbleEcology.

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/aee996ad-c758-4ee3-bd8f-3be2dbb6910fn%40googlegroups.com.

[emilygs...@gmail.com]

Hi Daniel,

Thank you those answers were really helpful.

Sorry, one more question. I wanted to create a matrix from the previous code, using some of the parameters to then calculate eigen values and vectors. 

The intention would be:

#transition_matrix[1,1] <- exp(alpha + (beta_age*1)) # fecundity for juveniles 
#transition_matrix[1,2] <- exp(alpha + (beta_age*2)) # fecundity for adults                          
#transition_matrix[2,1] <- mean_phi_juv # juvenile survival
#transition_matrix[2,2] <- mean_phi_adult # adult survival

result <- eigen(transition_matrix)

However, I am finding that I get errors relating to dimensionality. To my understanding, each input should be numerical of length 1 and this should create a 2x2 matrix. However, this does not seem to be the case. Do you have any idea why e.g. 'exp(alpha+(beta*1))' would not be length = 1 or why the code above would not give a matrix of 2x2? Is there a best practice to create a matrix from estimated parameters as part of a nimble model and estimate derived elements e.g. eigen values. 

If you need more info just let me know. I really do appreciate the help with this. 

Emily

[emilygs...@gmail.com]

So, I continued working on this and have solved some elements. I can now make the transition_matrix and when exploring the model after the model has been initialised with 'nimbleModel()' it seems to be a matrix with dim = 2 2. However, I still get errors when using the nimEigen function. 

I added indexes to the call transition_matrix like so:

eigen_values <- nimEigen(transition_matrix[1:2,1:2]) 

ERROR: Error in model$eigen_values[1] <<- nimEigen(model$transition_matrix[1:2,  :
  incompatible types (from S4 to double) in subassignment type fix

I am a bit at a loss as to why there is the incompatibility as my interrogation of this node in the model suggests it is the correct class and dimensions for nimEigen to work. 

[emilygs...@gmail.com]

Apologies for the thread chat with myself, but I'm hoping this might in the future help someone else!

I found the solution. 

The issue seemed to be that nimble was expecting the variable 'eigen_values' to have type = double. But the list output but the nimEigen() function did not match this. Therefore, I just extracted the actual value I wanted from the nimEigen() output and assigned that instead. Nimble then seems happy with this. At least, the model runs and seems to perform correctly. 

EXAMPLE:

eigen_values <- nimEigen(transition_matrix[1:2,1:2])$values[1]

[Daniel Turek]

Yes, that looks correct, Emily.  The output from the nimEigen function is a nimble-internal data structure which plays nicely with nimbleFunctions and the nimble compiler, which contains two internal elements:

nimEigen(...)$values

nimEigen(...)$vectors

These give a vector of the eigenvalues, and a matrix of the corresponding eigenvectors, respectively.  So, for extracting the first eigenvalue, it looks like you did exactly the right thing.

No problem about responding to your own emails, I'm glad you were able to trouble-shoot the problem.

Cheers,
Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/50eb5cfc-fa92-4bf0-9848-eb7ff087c2fdn%40googlegroups.com.