sim.capthist with sex-specific sigma

[PatrickCollins]

Hi there,

I am currently writing my thesis on camera trap designs for a region in the Western Cape, South Africa. I have hit a brick wall with regards to simulating a capture histories where sigma is dependent on sex. For context, the population is first simulated (with estimated parameters and a known pmix) before simulating the capture history and fitting the model with secr.fit and model args of D~elev and sigma~h2. I could not find any direct method of separating the sigma values by sex and so I tried ideas found on this message board which resulted in the following:

Simulating Population

Code: Select all

 my_pop <- sim.popn(D = D, core = masks, buffer = 0, Ndist = sim.popn.Ndist,poly =mask.poly$geometry,covariates = list(sex = c(M = 0.4375,F = 0.5625)))


Simulating Capture History

Code: Select all

    ch.M <- try(sim.capthist(traps = grids, pop = my_pop[covariates(my_pop)$sex=="M",],
                       noccasions = 1,
                       detectpar = list(lambda0 = 1, sigma =5100),
                       detectfn = "HHN"))
   
    ch.F <- try(sim.capthist(traps = grids, pop = my_pop[covariates(my_pop)$sex=="F",],
                       noccasions = 1,
                       detectpar = list(lambda0 = 1, sigma = 2360),
                       detectfn = "HHN"))
   
    covariates(ch.M)$sex <- try(factor(rep("M",dim(ch.M)[1])))
    covariates(ch.F)$sex <- try(factor(rep("F",dim(ch.F)[1])))
   
    ch.total<-try(MS.capthist(Females=ch.F,Males=ch.M,renumber=T))
    ch.total<-try(shareFactorLevels(ch.total))
    ch.total<-try(rbind(ch.total$Females,ch.total$Males))
   
    levels(covariates(ch.total)$sex)<-c("M","F")


Fitting Model

Code: Select all

 startvals <- list(D = mean(D), lambda0 = 1, sigma = sigma_pooled)
if (sum(covariates(ch.total)$sex=="F")==0){ #to avoid issues when no females caught
        model.args = list(D ~ elev, sigma~1)
      }else{
        model.args = list(D ~ elev, sigma~g)
      }
      mod <- try(secr.fit(capthist = ch.total,ncores=9,mask = masks_red, model = model.args, groups = "sex", binomN = secr.binomN, detectfn = "HHN", start = startvals, trace = FALSE, details = list(distribution = secr.dbn),method = "Nelder-Mead"))


This method resulted in decent estimations of the true simulated values of sigma but also resulted in systematically bad estimates of lambda0 as well as hilariously inconsistent results which make me inclined to believe that there is a huge issue with the way I've coded this (maybe theoretical, maybe code-based).

I would just like to ask if there's an easier or more straightforward method to do this. In the case that this method should be working, I'll then know to look for bugs in my code.

A sample of the results are shown below with estimates varying by the number of cameras used in the trap array where "B" is "Bias".



These may not seem abhorent but in we were getting some lambda0 estimates (using designs which are known to be good) with 200% bias or more. This, along with the fact that more cameras should mean better estimates which is not the case in my results, indicates to me that there is a structural issue at play. Any help will be appreciated beyond words.

edit: I should specify that the results shown are the aggregate of 50 iterations of the process described above.

[egc]

Murray will need to answer, but, thumbs up on using the code and tex markup tags. Made your post eminently readable.

[murray.efford]

Hi Patrick

Your code is more complicated than it needs to be, and there are lots of variables I can't see, so my mockup will differ. I agree simulating males and females separately, then combining is a bit clunky.

This looks OK to me:

Code: Select all

library(secr)
cameras <- make.grid(spacing = 2000, detector = 'proximity')
D <- 0.001
mask <- make.mask(cameras, buffer=10000)
   
my_pop <- sim.popn(D = D, core = mask, buffer = 0, Ndist = 'fixed',
    covariates = list(sex = c(M = 0.4375,F = 0.5625)))
# convert to factor; fix levels
covariates(my_pop)$sex <- factor(covariates(my_pop)$sex, levels = c('M','F'))

# use subset() to keep covariates
ch.M <- sim.capthist(traps = cameras, pop = subset(my_pop, covariates(my_pop)$sex=="M"),
    noccasions = 1, detectpar = list(lambda0 = 1, sigma = 5100), detectfn = "HHN")
ch.F <- sim.capthist(traps = cameras, pop = subset(my_pop, covariates(my_pop)$sex=="F"),
    noccasions = 1, detectpar = list(lambda0 = 1, sigma = 2360), detectfn = "HHN")

# rbind, not MS.capthist
ch.total <- rbind(ch.F, ch.M)
summary(ch.total)

secr.fit(ch.total, mask = mask, detectfn = 'HHN', groups = 'sex', model = sigma~g)


See if you can spot critical differences... I haven't used groups for > 10 years, but they still seem to work. Don't you expect bias in lambda0-hat if you fit a pooled sigma? How can you be sure the designs are good?

Murray

[murray.efford]

Actually, there may be a problem with reversing the factor levels. This works better.

Code: Select all

library(secr)
cameras <- make.grid(spacing = 2000, detector = 'proximity')
D <- 0.001
mask <- make.mask(cameras, buffer=10000)

set.seed(123)

my_pop <- sim.popn(D = D, core = mask, buffer = 0, Ndist = 'fixed',
    covariates = list(sex = c(F = 0.5625, M = 0.4375)))
covariates(my_pop)$sex <- factor(covariates(my_pop)$sex, levels = c('F','M'))

ch.M <- sim.capthist(traps = cameras, pop = subset(my_pop, covariates(my_pop)$sex=="M"),
    noccasions = 1, detectpar = list(lambda0 = 1, sigma =5100), detectfn = "HHN")
ch.F <- sim.capthist(traps = cameras, pop = subset(my_pop, covariates(my_pop)$sex=="F"),
    noccasions = 1, detectpar = list(lambda0 = 1, sigma =2360), detectfn = "HHN")

ch.total <- rbind(ch.F, ch.M)
summary(ch.total)

secr.fit(ch.total, mask = mask, detectfn = 'HHN', groups = 'sex', model = sigma~g)


[PatrickCollins]

Hi Dr Efford,

Thank you for your response, it has been extremely helpful in improving the model. Your advice with regards to altering the capture history, paired with swapping out groups for hcov, has resulted in the model excellently predicting lambda as well as and (which are now separated out). The issue now, which I can't seem to fix for the life of me, is that the density estimates are being systematically over-estimated. I have been working on this issue periodically for a few days now but to no avail. I honestly can't think of a reason for this, especially as all other predictions are good. If you have the time to give some advice, I'd greatly, greatly appreciate it. I've attached the code and results below.

For the sake of simplicity, I've only attached code used to fit the model of D~1, sigma~h2. Both constant density as well as spatially varying density will be considered in the thesis.

Simulating the Population

Code: Select all

        my_pop <- sim.popn(D = 0.000169, core = masks, buffer = 0, Ndist = "poisson",
poly =mask.poly$geometry,covariates = list(sex = c(M = 0.4375,F = 0.5625)))}

covariates(my_pop)$sex <- factor(covariates(my_pop)$sex, levels = c('F','M'))


Simulating the Capture History

Code: Select all

ch.M <- sim.capthist(traps = grids, popn = subset(my_pop, covariates(my_pop)$sex=="M"),
                       noccasions = 1,
                       detectpar = list(lambda0 = 1, sigma = 5100),
                       detectfn = "HHN")
   
    ch.F <- sim.capthist(traps = grids, popn = subset(my_pop, covariates(my_pop)$sex=="F"),
                       noccasions = 1,
                       detectpar = list(lambda0 = 1, sigma = 2360),
                       detectfn = "HHN")
   
    ch.total<-rbind(ch.F,ch.M)

Fitting the model

Code: Select all

   startvals <- list(D = 0.000169, lambda0 = 1, sigma = 3000, pmix=0.5)
masks_red<-masks

      mod <- try(secr.fit(capthist = ch.total , ncores=9 , mask = masks_red , model = list(D~1, sigma~h2) , hcov="sex" , binomN = 0,detectfn = "HHN" ,verify=T, start = startvals , trace = F , details = list(distribution = "poisson")))
     


Results
Please note that the results shown below are based off of the aggregation of 25 iterations of the entire process described above.



As can be seen, density seems to be consistently over-estimated while the other variables have reasonable bias. Additionally, an indication that there is issue is that the estimates don't necessarily improve with the number of cameras in the trap array which is rather counter-intuitive.

For the sake of completeness, I've attached an image of the 4 designs which the above results are based upon. Note that these plots impose the camera placements over the entire buffer used (which is why no cameras are near the edge)



Thank you for your time reading this and any further advice,
All the best,
Patrick

edit: The design method used is the one developed by Durbach et. al. (DOI:10.1111/2041-210X.13517) and is known to perform far better than is depicted above.

[murray.efford]

Hi Patrick

I can't diagnose your problem - I suspect it's in the inputs that we can't see (masks, grids).
I adapted your code for my own example and it looks fine (below).

An algorithmic approach to optimizing detector placement will be included in a future release of secrdesign and a development version is already on GitHub. The function GAminnr() implements Ian Durbach's approach.

Murray

Code: Select all

library(secr)
masks <- make.mask(type = 'polygon', poly = possumarea)
mod <- list()
for (r in 1:10) {
    grids <- make.systematic(n = 9, region = attr(masks,'polygon'),
        cluster = make.grid(nx=4, ny=4, spacing=60, detector='count'))
   
    my_pop <- sim.popn(D = 2, core = masks, buffer = 0, Ndist = "poisson",
        poly = attr(masks,'polygon'),
        covariates = list(sex = c(F = 0.5625, M = 0.4375)))
   
    covariates(my_pop)$sex <- factor(covariates(my_pop)$sex, levels = c('F','M'))
   
    ch.M <- sim.capthist(traps = grids, popn = subset(my_pop, covariates(my_pop)$sex=="M"),
        noccasions = 1, detectpar = list(lambda0 = 0.5, sigma = 60),
        detectfn = "HHN")
   
    ch.F <- sim.capthist(traps = grids, popn = subset(my_pop, covariates(my_pop)$sex=="F"),
        noccasions = 1, detectpar = list(lambda0 = 0.5, sigma = 40),
        detectfn = "HHN")
   
    ch.total<-rbind(ch.F,ch.M)
   
    mod[[r]] <- predict(secr.fit(capthist = ch.total , ncores=7 , mask = masks,
        model = list(D~1, sigma~h2) , hcov="sex" , binomN = 0, detectfn = "HHN" ,
        verify = FALSE, trace = FALSE ,
        details = list(distribution = "poisson")))
   
    cat("Completed ", r, "\n")
}

mod[[1]]
extract <- function(x, parm='D') x[[1]][parm,'estimate']
mean((sapply(mod, extract, 'D') - 2.0) / 2.0)
# [1]0.06192489


[PatrickCollins]

Hi Murray,

Thank you for all of your assistance. Our main issue was with the stability of our estimates as we weren't correctly reading in the starting values. The median improved our results (over the mean) as well as the proper implementation of the starting values to the point where we would get biases for all parameters below 5%.

This isn't a post asking for help, I just thought it would be rude not to give the update. Additionally, would you consent to us thanking you in our acknowledgements?

All the best,
patrick