Questions multi-session mark-resight models
Vanessa
Hello,
I am using the SECR package in order to estimate the density of wild dogs for a camera-trap study, with mark-resight multi-session models. The data spans over an entire year, divided into 4 sessions. The first 2 sessions correspond to the dry season and the last 2 sessions correspond to the wet season.
When running the multi-session models, I always need to specify autoini on session 2. Specifying autoini on any of the other sessions results in autoini failure and NA results for all estimates. So far, the multi-session models (using autoini on session 2) seem to work, although I sometimes get warning messages of the type ‘possible maximization error: nlm returned code 3’ or ‘possible maximization error: nlm returned code 4. See ?nlm’ or ‘probable maximization error: optim returned convergence 1’ or ‘at least one variance calculation failed’. I have ignored these error messages in some cases (or tried some tricks from the troubleshooting vignette in other cases), but these messages might be indicating more serious problems with my models. Some examples are provided below.
Here are the main issues that I have encountered with my models:
1. I obtain a value of pID = 1.0 for most of the models that I have ran so far. From my understanding, the argument details = list(knownmarks = FALSE) is necessary in my model since we do not know the number of marked animals in the population (we are using natural marks to identify individuals). When I remove the argument details = list(knownmarks = FALSE) from secr.fit, the pID becomes <1.0 but g0 then becomes equal to 1.0. Is this a normal behavior and if not, what does this indicate?
I realised that including the argument method = “Nelder-Mean” seems to result in values for both pID and g0 lower than 1.0. Should I be using this type of method instead? I have read that this method is more robust than the default, but it also provides different results. Unfortunately, the Nelder-Mean method sometimes results in NaNs of SE estimates, as you can see below with the HN function, as an example. (The EX function seems to work when using the method “Nelder-Mean”).
Here are examples of the 3 different scenarios, with the HN function. I provide only results of session 1, as an example:
> fit.HN <- secr.fit(MRCH_MR, model = D~session, mask=clippedmask_MR_cov, detectfn = 'HN', trace = FALSE, details = list(autoini = 2, knownmarks = FALSE))
Warning message:
In secr.fit(MRCH_MR, model = D ~ session, mask = clippedmask_MR_cov, :
possible maximization error: nlm returned code 3. See ?nlm
> predict(fit.HN, all.levels = TRUE)
$`session = 1`
link estimate SE.estimate lcl ucl
D log 1.680162e-03 3.979370e-05 1.603961e-03 1.759984e-03
g0 logit 2.002803e-01 1.339174e-02 1.753150e-01 2.278186e-01
sigma log 1.020530e+03 2.650073e+01 9.698975e+02 1.073806e+03
pID logit 1.000000e+00 7.641795e-18 1.000000e+00 1.000000e+00
> fit.HN_without_knownmarks <- secr.fit(MRCH_MR, model = D~session, mask=clippedmask_MR_cov, detectfn = 'HN', trace = FALSE, details = list(autoini = 2))
Warning message:
In secr.fit(MRCH_MR, model = D ~ session, mask = clippedmask_MR_cov, :
possible maximization error: nlm returned code 4. See ?nlm
> predict(fit.HN_without_knownmarks, all.levels = TRUE)
$`session = 1`
link estimate SE.estimate lcl ucl
D log 6.076282e-04 1.095946e-06 6.054840e-04 6.097801e-04
g0 logit 1.000000e+00 6.018427e-16 1.000000e+00 1.000000e+00
sigma log 1.340144e+03 4.538150e+01 1.254109e+03 1.432081e+03
pID logit 6.803114e-02 6.415837e-03 5.648349e-02 8.173512e-02
> fit.HN_Nelder <- secr.fit(MRCH_MR, model = D~session, method = "Nelder-Mead",mask=clippedmask_MR_cov, detectfn = 'HN', trace = FALSE, details = list(autoini = 2, knownmarks = FALSE))
Warning messages:
1: In secr.fit(MRCH_MR, model = D ~ session, method = "Nelder-Mead", :
probable maximization error: optim returned convergence 1. See ?optim
2: In secr.fit(MRCH_MR, model = D ~ session, method = "Nelder-Mead", :
at least one variance calculation failed
> predict(fit.HN_Nelder, all.levels = TRUE)
$`session = 1`
link estimate SE.estimate lcl ucl
D log 1.628452e-03 3.430370e-05 1.562594e-03 1.697086e-03
g0 logit 2.742999e-01 NaN NaN NaN
sigma log 1.004315e+03 2.707621e+01 9.526338e+02 1.058801e+03
pID logit 8.086541e-01 NaN NaN NaN
I tried calling secr.fit with method = “none” and using the start values from the previous fit, as explained in the vignette ‘Troubleshooting’, but this results in NaNs for all SE estimates (including D and sigma).
2. When I include the session-specific covariate season in the model (dry for the first 2 sessions, followed by wet for the last 2 sessions), the results for the dry and wet seasons repeat themselves for each session (1,2,3,4) as if there was a dry and wet season for each session (which is not the case). Below, I provide only the results of the first 2 sessions as an example. Did I include the session covariate correctly and can I trust my results?
session_cov <- data.frame(season = factor(c('dry','dry','wet','wet')))
fit.density_Season <- secr.fit(MRCH_MR, model = D~season, sessioncov = session_cov, mask=clippedmask_MR, trace = FALSE, details = list(autoini = 2, knownmarks = FALSE))
predict(fit.density_Season, all.levels = TRUE)
$`session = 1, season = dry`
link estimate SE.estimate lcl ucl
D log 1.833007e-03 4.480641e-05 1.747271e-03 1.922950e-03
g0 logit 1.339645e-01 8.042522e-03 1.189702e-01 1.505257e-01
sigma log 1.057070e+03 2.695429e+01 1.005547e+03 1.111232e+03
pID logit 1.000000e+00 7.965126e-20 1.000000e+00 1.000000e+00
$`session = 1, season = wet`
link estimate SE.estimate lcl ucl
D log 1.708045e-03 4.000192e-05 1.631425e-03 1.788264e-03
g0 logit 1.339645e-01 8.042522e-03 1.189702e-01 1.505257e-01
sigma log 1.057070e+03 2.695429e+01 1.005547e+03 1.111232e+03
pID logit 1.000000e+00 7.965126e-20 1.000000e+00 1.000000e+00
$`session = 2, season = dry`
link estimate SE.estimate lcl ucl
D log 1.833007e-03 4.480641e-05 1.747271e-03 1.922950e-03
g0 logit 1.339645e-01 8.042522e-03 1.189702e-01 1.505257e-01
sigma log 1.057070e+03 2.695429e+01 1.005547e+03 1.111232e+03
pID logit 1.000000e+00 7.965126e-20 1.000000e+00 1.000000e+00
$`session = 2, season = wet`
link estimate SE.estimate lcl ucl
D log 1.708045e-03 4.000192e-05 1.631425e-03 1.788264e-03
g0 logit 1.339645e-01 8.042522e-03 1.189702e-01 1.505257e-01
sigma log 1.057070e+03 2.695429e+01 1.005547e+03 1.111232e+03
pID logit 1.000000e+00 7.965126e-20 1.000000e+00 1.000000e+00
3. I am trying to include a station X occasion covariate (using timevaryingcov). For each session, the trap file contains a column for the station ID, two columns for the x and y coordinates, followed by 1 column for each occasion of this time varying covariate. (The first column of the covariate contains the forward slash ‘/’ in front of all values for that first occasion). The following code results in an error message, and I cannot figure out what is my mistake.
MRCH_MR <- read.capthist(captfile, c(trapfile_cov3_s1,trapfile_cov3_s2, trapfile_cov3_s3, trapfile_cov3_s4), detector = "count", markocc = list(c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)), binary.usage = FALSE)
MRCH_MR <- addSightings(MRCH_MR, unmarked = list(UM_s1, UM_s2, UM_s3, UM_s4), nonID = NULL, uncertain = NULL, verify = TRUE)
timevaryingcov(traps(MRCH_MR$'1')) <- list(obstruction = 1:17)
timevaryingcov(traps(MRCH_MR$'2')) <- list(obstruction = 1:15)
timevaryingcov(traps(MRCH_MR$'3')) <- list(obstruction = 1:15)
timevaryingcov(traps(MRCH_MR$'4')) <- list(obstruction = 1:15)
> fit.g0_Obstr <- secr.fit (MRCH_MR, mask = clippedmask_MR_cov, detectfn = 'HN', model = g0 ~ obstruction, trace = FALSE, details = list(autoini = 2, knownmarks = FALSE))
> predict(fit.g0_Obstr)
Error in secr.lpredictor(formula = models[[x]], newdata = newdata, indx = parindices[[x]], :
model covariate obstruction not found in 'newdata'
4. I get slightly different results when I run the exact same models. Here is an example below. Are these differences considered normal?
> fit.HN_MR <- secr.fit (MRCH_MR, mask = clippedmask_MR_cov, trace = FALSE, details = list(autoini = 2,knownmarks = FALSE))
> predict(fit.HN_MR, all.levels = TRUE)
$`session = 1`
link estimate SE.estimate lcl ucl
D log 1.809579e-03 4.119627e-05 1.730621e-03 1.892140e-03
g0 logit 1.254606e-01 7.276727e-03 1.118788e-01 1.404304e-01
sigma log 1.072931e+03 2.648941e+01 1.022257e+03 1.126118e+03
pID logit 9.999998e-01 2.459475e-15 9.999998e-01 9.999998e-01
> fit0.usage_MR2 <- secr.fit(MRCH_MR, mask=clippedmask_MR, trace = FALSE, details = list(autoini = 2, knownmarks = FALSE))
Warning message:
In secr.fit(MRCH_MR, mask = clippedmask_MR, trace = FALSE, details = list(autoini = 2, : possible maximization error: nlm returned code 3. See ?nlm
> predict(fit0.usage_MR2)
$`session = 1`
link estimate SE.estimate lcl ucl
D log 1.809332e-03 4.119854e-05 1.730369e-03 1.891897e-03
g0 logit 1.254171e-01 7.275238e-03 1.118383e-01 1.403841e-01
sigma log 1.073127e+03 2.650127e+01 1.022430e+03 1.126338e+03
pID logit 1.000000e+00 2.804601e-19 1.000000e+00 1.000000e+00
5. Similarly, I tried comparing the following model results, which, in theory, should give the same results:
- the multi-session model, while letting all parameters vary (density, g0, sigma and pID). Examples below for sessions 1 and 2:
> fit.EX_all <- secr.fit(MRCH_MR, , detectfn = 'EX', model = list(D~session, g0~session, sigma~session, pID~session), mask=clippedmask_MR_cov, trace = FALSE, details = list(autoini = 2, knownmarks = FALSE))
Warning message:
In secr.fit(MRCH_MR, , detectfn = "EX", model = list(D ~ session, :
possible maximization error: nlm returned code 4. See ?nlm
> predict(fit.EX_all)
$`session = 1`
link estimate SE.estimate lcl ucl
D log 2.138799e-03 1.270303e-04 1.903964e-03 2.402598e-03
g0 logit 4.525159e-01 4.721200e-02 3.626191e-01 5.456211e-01
sigma log 4.732434e+02 3.080375e+01 4.166178e+02 5.375655e+02
pID logit 9.976084e-01 1.844694e-02 1.093749e-04 1.000000e+00
$`session = 2`
link estimate SE.estimate lcl ucl
D log 1.581688e-03 1.158562e-04 1.370424e-03 0.00182552
g0 logit 2.075733e-01 3.549257e-02 1.464735e-01 0.28563105
sigma log 7.995948e+02 5.996985e+01 6.904291e+02 926.02098939
pID logit 9.994117e-01 7.616414e-03 1.599157e-08 1.00000000
- the models for each session individually. Models from sessions 1, 3 and 4 do not converge (I obtain NA results for all estimates). I therefore use the start argument to fix the value of g0 to the same value as in the results from the multi-session model. Examples below for sessions 1 and 2:
> fit0.EX_s1 <- secr.fit (MRCH_MR_s1, mask = clippedmask_MR_s1, start = list(g0=0.4525159), trace = FALSE, detectfn ='EX', details = list(knownmarks = FALSE))
Warning message:
In autoini(ch, msk, binomN = details$binomN, ignoreusage = details$ignoreusage) :
'autoini' failed to find g0; setting initial g0 = 0.1
> predict(fit0.EX_s1)
link estimate SE.estimate lcl ucl
D log 1.578757e-03 1.059090e-04 1.384451e-03 1.800334e-03
g0 logit 1.000000e+00 6.269059e-16 1.000000e+00 1.000000e+00
sigma log 5.264464e+02 3.446283e+01 4.631177e+02 5.984349e+02
pID logit 3.349430e-01 3.876866e-02 2.636611e-01 4.146443e-01
> fit0.EX_s2 <- secr.fit (MRCH_MR_s2, mask = clippedmask_MR_s2, trace = FALSE, detectfn ='EX', details = list(knownmarks = FALSE))
Warning message:
In secr.fit(MRCH_MR_s2, mask = clippedmask_MR_s2, trace = FALSE, :
possible maximization error: nlm returned code 3. See ?nlm
> predict(fit0.EX_s2)
link estimate SE.estimate lcl ucl
D log 1.579104e-03 1.130677e-04 1.372589e-03 0.00181669
g0 logit 2.082713e-01 3.626389e-02 1.459886e-01 0.28815986
sigma log 8.000058e+02 6.040849e+01 6.900971e+02 927.41916965
pID logit 1.000000e+00 5.911787e-20 1.000000e+00 1.00000000
As you can see, results for session 1 are quite different when comparing both methods and results for session 2 are similar but not exactly the same as I would expect. Is this normal?
6. When using a 2-class mixture model (h2), is it possible to identify which marked individuals were attributed to each group? I would like to explore if sigma varies spatially, i.e. if sigma is different for wild dogs close to communities vs wild dogs far from communities. I do not see how I could explore this idea other than introducing the covariate h2 (but at the moment, although the model seems to work, I cannot tell what these two groups represent). If this is not possible, would there be an alternative way to test this hypothesis?
Any help would be very much appreciated. I am happy to share my data, but by email privately, if possible.
Thank you very much in advance!
All the best,
Vanessa
Murray Efford
Vanessa
Hello Murray,
Thank you so much for your answers. I appreciate your help. I just need further clarifications on two of these issues:
Regarding point #1, I am still unsure if obtaining the pID = 1.0 is indicative of model fitting issues as well? I get these results (pID = 1.0) in most of my models, many of which did not give any warning or error messages. Can I trust my results when I get this value (without error/warning messages), and if not, do you have any suggestions to bypass this problem?
Regarding point #3 (timevaryingcov), I tried using a one-session, capture-recapture model, rather than multi-session mark-resight model, and the same error message pops up each time:
Error in secr.lpredictor(formula = models[[x]], newdata = newdata, indx = parindices[[x]], :
model covariate obstruction not found in 'newdata'
This error message is therefore not related to my mark-resight models being fragile. I am assuming that I am using the wrong approach. Below, I copied my code. Would you have an idea what might be causing these error messages?
These are the two remaining obstacles that are holding me back.
Thank you again for your help,
Vanessa
> timevaryingcov(traps(MRCH_CR_s1)) <- list(obstruction = 1:17)
> fit.CR_obstruction_s1 <- secr.fit(MRCH_CR_s1, mask=clippedmask_CR_cov_s1, model = g0 ~ obstruction, trace = FALSE)
> predict(fit.CR_obstruction_s1)
Error in secr.lpredictor(formula = models[[x]], newdata = newdata, indx = parindices[[x]], :
model covariate obstruction not found in 'newdata'
Based on the steps provided in the ’Datainput’ vignette, my trap file follows the format:
Station01 649864 8789222 /0 0 0.2 0.2 0 0.2 0 0 0 0 0 0 0 0 0 0 0
Station02 648892 8790539 /0 0 0 0 0 0 0 0 0 0 0.4 0 0.5 0.7 0 0 0
Station03 647938 8792038 /0 0 0 0 0 0 0 0 0 0 0 0 0 0.4 0.5 0.2 0
Etc…
(all 17 columns following the forward slash (/) represent the time varying covariate, obstruction, one for each occasion)
Murray Efford
--
You received this message because you are subscribed to the Google Groups "secr" group.
To unsubscribe from this group and stop receiving emails from it, send an email to secrgroup+...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/secrgroup/35cfb4ad-2aee-4a37-9162-2d24e9a66ed6%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
Murray Efford
On Wednesday, July 3, 2019 at 11:06:14 AM UTC+12, Murray Efford wrote:
Point #1I'm really not able to help. I don't use mark-resight models myself. Others may have relevant experience. Sometimes it can help to fix parameters (the 'fixed' argument of secr.fit) to see what effect that has on other estimates and the logLik.Point #3This should be easy to fix - just specify 'obstruction' in the newdata argument of predict. See Examples in ?timevaryingcov. The fact that it is not included by default may be a bug, or maybe there was a good reason... will check.
To unsubscribe from this group and stop receiving emails from it, send an email to secrgroup+unsubscribe@googlegroups.com.
Vanessa
Hi Murray,
Thanks for the clarifications. For point #3, specifying my covariate in newdata easily resolved the problem, so that seems all good! Thanks.
For point #1, I tried to fix the pID parameter to different values. In most of the cases, the model provides an estimated pID equal to 1 (even though I fixed it at lower values, such as 0.01 for example), which results in identical density estimates and loglik. Whenever the estimated pID obtained is lower than 1, the resulting density estimates are different and the loglik (and AIC) becomes a lot larger.
I would be happy to hear any further thoughts that you, or anybody else in this group, may have about this last point. If not, do you know a person who might be able to help me with understanding the behavior of these mark-resight models?
Thanks again,
Vanessa
To unsubscribe from this group and stop receiving emails from it, send an email to secr...@googlegroups.com.
Vanessa
> fit.HN_e2_DSes_g0Vis_2 <- secr.fit(MRCH_e2_2, model = list(D~session, g0 ~ visibility), mask=clippedmask_e2_cov, detectfn = 'HN', trace = FALSE, details = list(autoini = 2, knownmarks = FALSE))
> predict(fit.HN_e2_DSes_g0Vis_2, newdata = data.frame(visibility = 0:1, session = c("S1","S2")))
Error in mat %*% beta[indx] : non-conformable arguments
> predict(fit.HN_e2_DSes_g0Vis_2, newdata = data.frame(visibility = 0:1, session = "S1"))
Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :
contrasts can be applied only to factors with 2 or more levels
Murray Efford
Vanessa
Vanessa
Dear Murray,
Thank you very much for all of your help, it has been very useful.
I just wanted to follow-up on this particular point, from the first messages sent (My question is in gray and your answer is in blue):
6. When using a 2-class mixture model (h2), is it possible to identify which marked individuals were attributed to each group? I would like to explore if sigma varies spatially, i.e. if sigma is different for wild dogs close to communities vs wild dogs far from communities. I do not see how I could explore this idea other than introducing the covariate h2 (but at the moment, although the model seems to work, I cannot tell what these two groups represent). If this is not possible, would there be an alternative way to test this hypothesis?
6. 'secr' doesn't do this (note lower case). I included a function classMembership in openCR that does this for the different set of models in that package. The procedure is to estimate the likelihood components for each animal given membership in each class and assign a proportional probability of membership. Although in principle it could be implemented in 'secr' the coding is very different so it won't happen soon. On the other hand, you should not expect that the classes align with biological groups - they needn't 'represent' anything other than variability in the parameter. One might address your hypothesis of spatial variation in sigma via non-euclidean distance models, but I'm not going to push that as I foresee other complications and the signs are you have limited data..
What are your reasons for suggesting to not use non-euclidean distance models and what complications do you foresee? Are the complications related to the use of multi-sessions, mark-resight models with this approach?
Thanks!
Vanessa
Murray Efford
Vanessa
Hi Murray,
Thanks for your answer.
I am trying to fit a model, in which sigma and density can vary according to a two-class habitat covariate (which I already included into the habitatmask, called ‘Community’). Following the non-euclidean secr vignette, I seem to understand that example #5 corresponds exactly to my situation. Is that correct?
I used the function fn2 from your vignette:
fn2 <- function (xy1, xy2, mask) {
if (missing(xy1)) return("noneuc")
xy1 <- addCovariates(xy1, mask)
xy2 <- addCovariates(xy2, mask)
sig1 <- as.numeric(covariates(xy1)$noneuc) # sigma(x,y) at detectors
sig2 <- as.numeric(covariates(xy2)$noneuc) # sigma(x,y) at mask points
euc <- edist(xy1, xy2)
sig <- outer (sig1, sig2, FUN = function(s1, s2) (s1 + s2)/2)
euc / sig
}
And used the following code for fitting my model:
fit.HN_e2_DHab_SigHab <- secr.fit(MRCH_e2, model = list(D~ Community, noneuc ~ Community), mask=clippedmask_e2_cov, detectfn = 'HN', trace = FALSE, details = list(autoini = 2, knownmarks = FALSE, userdist = fn2), fixed = list(sigma = 1))
predict(fit.HN_e2_DHab_SigHab, all.levels= TRUE)
After only a few seconds of process, I obtain NA results for all estimates. The same problem occurs if I only make sigma vary with community. I was therefore wondering if this meant that the model did not converge, or if my approach or code was incorrect.
Also, can I compare AIC values between models based on non-euclidean distances with the other standard models?
Thanks again!
Vanessa
Vanessa
Hi Murray,
I was able to obtain results by assigning a higher value to the fixed argument of sigma. I still need clarifications regarding this method, and will reformulate and simplify my last question, as followed:
- Is it possible to create a model in which sigma will vary from one habitat to another (i.e. from a two-class spatial covariate, not a continuous scale), while using an Euclidean distance calculation (i.e. without avoiding wholes that might exist in the habitat mask)? If so, how can this be done?
- Would this method assume that the animals have the same easiness to move in both types of environment? In our study area, the home range of animals might be larger or smaller due to food resource availability, only. It seems that the landscape connectivity theory from example 5 would assume that movement of animals would be altered (become easier or more difficult) depending on landscape structure, which is not the case in our situation.
- Can the AIC values from this model be compared with other standard Euclidean-distance based models?
Thank you very much in advance for your help!
Vanessa