sim.capthist with sex-specific sigma

questions concerning anlysis/theory using program DENSITY and R package secr. Focus on spatially-explicit analysis.

sim.capthist with sex-specific sigma

Postby PatrickCollins » Tue Oct 11, 2022 1:47 pm

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".

\begin{center}
\begin{tabular}{ l | c | c | c | c | c | c | c | } 
Cams & D & B(D) & lam0 & B(lam0) & sig_{pooled} & B(sig)  \\
 \hline
 \hline
36 & 1.82 & 8\% & 1.34 & 34\% & 3261 & 9\% \\ 
43 & 1.78 & 5\% & 1.28 & 28\% & 2910 & -3\%  \\
64 & 2.11 & 25\% & 0.94 & -6\% & 3237 & 8\% \\
81 & 1.56 & -8\% & 0.97 & -3\% & 3117 & 4\% \\
\hline
True Value & 1.69 &  & 1 & & 3000 &   \\
\hline
\end{tabular}
\end{center}

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.
Last edited by PatrickCollins on Tue Oct 11, 2022 3:44 pm, edited 2 times in total.
PatrickCollins
 
Posts: 3
Joined: Tue Oct 11, 2022 10:48 am

Re: sim.capthist with sex-specific sigma

Postby egc » Tue Oct 11, 2022 2:50 pm

Murray will need to answer, but, thumbs up on using the code and tex markup tags. Made your post eminently readable.
egc
Site Admin
 
Posts: 201
Joined: Thu May 15, 2003 3:25 pm

Re: sim.capthist with sex-specific sigma

Postby murray.efford » Tue Oct 11, 2022 7:56 pm

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
 
Posts: 686
Joined: Mon Sep 29, 2008 7:11 pm
Location: Dunedin, New Zealand

Re: sim.capthist with sex-specific sigma

Postby murray.efford » Tue Oct 11, 2022 8:03 pm

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)
murray.efford
 
Posts: 686
Joined: Mon Sep 29, 2008 7:11 pm
Location: Dunedin, New Zealand

Re: sim.capthist with sex-specific sigma

Postby PatrickCollins » Fri Oct 14, 2022 3:47 pm

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 sigma_{male} and sigma_{female} (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.

Image

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)

Image

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.
PatrickCollins
 
Posts: 3
Joined: Tue Oct 11, 2022 10:48 am

Re: sim.capthist with sex-specific sigma

Postby murray.efford » Sun Oct 16, 2022 1:42 pm

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
murray.efford
 
Posts: 686
Joined: Mon Sep 29, 2008 7:11 pm
Location: Dunedin, New Zealand

Re: sim.capthist with sex-specific sigma

Postby PatrickCollins » Sun Oct 23, 2022 9:24 am

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
PatrickCollins
 
Posts: 3
Joined: Tue Oct 11, 2022 10:48 am


Return to analysis help

Who is online

Users browsing this forum: No registered users and 15 guests