Hi there!
I'm stuck in building a CMR multistate model with a spatial component and I would be very grateful if somebody could provide some help.
My first aim is to generate dispersal distances for all individuals in my study split by sex on the basis of a gamma distribution. Some individuals have observed dispersal distances, but it has not been observed for some others. Therefore, the "Dist" vector is a combination of observed values and NAs.
for(i in 1:nind){
Dist[i] ~ dgamma(sh[sexvec[i]], ra[sexvec[i]]) #Sexvec: vector indicating the sex of #each individual
}
for (g in 1:2){ #Priors
sh[g] ~ dgamma(1,1)
ra[g] ~ dgamma(1,1)
}
The second step of the model would be to compare the dispersal distances of each individual with available recruitment points (which are known with certainty) and allocate every individual to the closest point to each dispersal distance. To do so, I entered an object to the model indicating the relative distances between each birth location and each potential recruitment point (Distmat), and a second object indicating the ID of the birth point of each individual (birthterr). Then, this part of the model must assign every individual to the closest recruitment point and return the ID of the recruitment point of each individual.
Perry de Valpine kindly provided some feedback of this approach and suggested a function a few weeks ago:
sortbyterr <- nimbleFunction(
run = function(Distmat = double(2), birthterr = double(1), Dist = double(1), nterr = double()){
returnType(double(1))
n1 <- length(Dist)
recterr <- numeric(length=n1, init=FALSE)
for(i in 1:n1) {
min_val <- Inf
min_ind <- 1
for(j in 1:nterr) { #nterr just indicates the total number of available territories
this_val <- abs(Distmat[birthterr[i], j] - Dist[i])
if(this_val < min_val) {
min_val <- this_val
min_ind <- j
}
}
recterr[i] <- min_ind
}
return(recterr)
})
The function is incorporated into the model, and a second kernel (newkernel) is later calculated with the distances between birth and recruitment points.
recterr[1:nterr] <- sortbyterr(Distmat[1:51, 1:nterr], birthterr[1:nind], Dist[1:nind], nterr)
for(i in 1:nind){
newkernel[i] <- Distmat[birthterr[i],recterr[i]]
}
However, it seems there might be a problem in which the output of the function ("recterr") might not be created successfully. When compiling, I am warned that "dynamic index out of bounds: Distmat[birthterr_oBi_cB, recterr[i]]".
If that helps, Distmat is a 51*132 matrix and birthterr is length 461 with values between 1 and 51. 461 is the number of individuals in my study, 132 the number of points of potential recruitment, 51 is the number of points of birth.
Further errors are reported when trying to run the model.
Error in if(this_val < min_val){:
missing value where TRUE/FALSE needed
Note: cannot calculate logProb for node recterr[1:132]
Error in mapCopy. Sizes don't match: 132 != 461.
The model will ultimately crash. I've been trying with other equivalent structures that incorporate the resulting "recterr" object and all provide errors I cannot identify as well.
An alternative approach I've been trying to get this recterr object is to work with different built-in functions outside nimbleFunctions. Once the Dist object is calculated exactly identically to the first time, I1) calculate the difference between all distances in Distmat and the Dist object, 2) save the minimum value of each row, 3) use the built-in function "equals" to assign a value of 1 to the row with the minimum value, and 0s elsewhere, and 4) use the built-in function "inprod" to obtain the territory number of the closest territory to the simulated Dist value.
for(i in 1:nind){
for(j in 1:nterr){
diffmat[i,j] <- abs(Distmat[birthterr[i],j] - Dist[i])
}
}
for(i in 1:nind){
minvec[i] <- min(diffmat[i,1:nterr])
}
for(i in 1:nind){
for(j in 1:nterr){
e[i,j] <- equals(diffmat[i,j], minvec[i])
}
}
for(i in 1:nind){
recterr[i] <- inprod(e[i,1:nterr],1:nterr)
}
Here my model complains about e[i,j] calculated by function equals, with several warnings indicating that:
Error in if (x1 == x2) 1 else 0 :
missing value where TRUE/FALSE needed
Note: cannot calculate logProb for node e[(...), (...)].
The model will crash as well at some point following this structure. Is it something wrong with my usage of the equals function?
My plan is to structure survival and recapture by point of recruitment in further steps, so I'll need to use this object as a dynamic index many times too.
Does anybody have any idea why either of the approaches fail? Any guess would be very appreciated. I tried to be as informative as possible but please let me know if you need any more information or remarks.
Thank you guys and have a nice day!
Jaume A.