Estimates, AIC when variance caclulation fails

[lourens]

Hi All,
I have seen on previous posts regarding AIC and estimates when variance calculation failed [url]
viewtopic.php?f=36&t=2382[/url] I have data on mark recapture for 3 different sites, where I modelled g0 and sigma affected by sex, b, and h; here is my code

Code: Select all

M0.welg<-secr.fit (welg.ch, model = g0~1, mask=mask.welg,CL=T, trace = FALSE)
Mb.welg<-secr.fit (welg.ch, model = g0~b,CL=T, mask=mask.welg, trace = TRUE, method = "Nelder-Mead")
Mh.welg<-secr.fit (welg.ch, model = g0~h2, mask=mask.welg,CL=T,trace = FALSE, details=list(hessian = "fdhess"))
Msex.welg<-secr.fit(welg.ch, model = g0~Sex, mask=mask.welg, CL=T, trace = FALSE)
Msex2.welg<-secr.fit(welg.ch, model = list(g0~Sex, sigma~Sex), mask=mask.welg, CL=T, trace = FALSE)


My question is as follows;
For one site, I get complete answers and estimates for all these models, while for the other 2 sites I get non-sense answers for both the heterogeneity (h) and behaviour (b) models (density estimates over milion etc). On these two sites it seems that the mark recapture data is not adequate for such an analysis. My question is how I report this or what do I do. I don’t think it is a maximisation problem, but rather a data problem.
Any suggestions and help will be welcome
Lourens

[murray.efford]

Hi Lourens
I can't solve your data problems, if that's what they are, but maybe these points will help
1. Make sure you have secr 2.6.0. This uses a more robust matrix inversion that helps when some parameters are not identifiable.
2. The 'approved' way of dealing with a mixture of adequate and inadequate data is to pool and declare the assumption that detection parameters were the same across sites. The mechanism for this in 'secr' is to treat sites as sessions in a multi-session capthist object. Then any secr.fit that does not specify 'session', or another session covariate, in the model argument implicitly pools data across sessions (and you'll see identical detection parameter estimates for each session = site in the output).
3. secr 2.6.0 has a new mechanism for dealing with sex differences that is not yet fully documented, but here is how it works:

Code: Select all

Msex.welg.h2 < -secr.fit(welg.ch, model = g0~h2, hcov = 'Sex', mask = mask.welg, CL = TRUE, trace = FALSE)
Msex2.welg.h2 < -secr.fit(welg.ch, model = list(g0~h2, sigma~h2), hcov = 'Sex', mask = mask.welg, CL = TRUE, trace = FALSE)

The fitted model estimates sex-specific detection parameters and the sex ratio. This won't solve your lack of data, but it's tidy and also allows some 'Sex' values to be NA. The catch: you need to include hcov = 'Sex' in all models you want to compare by AIC.
Murray

[lourens]

Hi Murry,
Thank you so much for the reply. I have downloaded version 2.6.0 and implemented the suggestions and pooled data across the different sites. Just to run code by you to make sure I have done it correctly:
I first made each reserve into a session, since they dont have the same number of occations;
then I made a session data frame where session 1 is reserve 1 etc;

Code: Select all

# rbind sessions, session1 is reserve 1
comb<-MS.capthist(session1.ch,session2.ch,session3.ch)
# session df
session.df <- read.table("session.txt",header=T,sep=" ",na.strings="NA", strip.white=TRUE)
# fit first model
model0 <- secr.fit(comb, model=g0~1, mask=list(mask.lap,mask.welg, mask.gf),CL=T,sessioncov = session.df,trace = TRUE)


While this works and variance calculations succeeded, I still had non-sense results for model=g0~b; e.g. 25000 animals, which does not make sense.
Looking at the data I think it is related to the low recaptures I had for 2 of the reserves. So, I quess I am out of luck.
Out of interest Tobler et al 2013 http://www.sciencedirect.com/science/article/pii/S0006320712005101 suggest to "borrow" correct values for sigma from study sites with adequate and more robust data and to fix these for other study sites with troubled data sets. I know by pooling the data as you suggested is doing something similar, but would fixing sigma be a better approach? In such a case how do you fix parameter sigma for each sex?
Also, the hcov='Sex' is still a bit unclear for me, e.g. if I had a model=list(g0~Sex, sigma~Sex), would it still be

Code: Select all

test<-secr.fit(ch, model=list(g0~sex, sigma~sex), hcov='Sex', mask=mask)

Lastly, which I should have done in first post is state my question:
I wanted to create 5 models (Model=g0~1, model=g0~b, model=g0~h2, model=g0~Sex, model=list(g0~Sex, sigma~Sex) to see if including sex in sigma and g0 produced better models. As I said before, for the 1 reserve with adequate data this worked well, but at other 2 model=g0~b, model=g0~h2 did not work.
Thanks again for the help.
Lourens

[murray.efford]

Hi Lourens

In your first block of code I don't understand what you are trying to do with session.df. At present it does nothing in the secr.fit (you name it as the sessioncov dataframe, but don't use any covariates in the model formulae).

The Tobler et al approach applied to MLE is more 'stealing' than 'borrowing'. It's something you might do to explore reasons for model failure, but should avoid for reporting results because it pretends that sigma is known without any error (this may not be the case in a Bayesian setting if you use the posterior from the first fit... I can't remember what they did). If you must fix the sexes separately then the procedure is rather technical and I'd rather dodge the question! (Hint: you can set details = list(fixedbeta = c(NA,NA,a,b-a,NA)) to fix the coefficients for sex-specific sigma where a = log(sigma1) and b = log(sigma2) in a CL g0~h2, sigma~h2 hcov model - aren't you sorry you asked?!)

For more on hcov and other ways to handle sex, see my post to Isoodon. Note the hcov bug in derived() when CL=T, to be fixed soon.

[lourens]

Hi Murray,
Once again, thanks for the help. R kept on crashing yesterday using derived (), CL and hcov, so I am relieved to see it’s a bug, I normally assume my code is incorrect!. I think I will rather stay with the pooled method.
Regards.
Lourens

[murray.efford]

I've submitted a corrected version (2.6.1) to CRAN - it will be available as Windows & Mac binaries in a day or two.
Murray