Simulations using groups for specifying sex-specific parameters

[Rebecca G]

Hi Murray,

I'm running some simulations to test some new proposed trap layouts with different parameters for males and females. Here is my code (using the example given in the secrdesign vignette, thanks for that).

library(secrdesign)
library(rgdal)
ciArea <- readOGR("ChurchillIs_1m_elev_area.shp")
ciScenarios <- make.scenarios(D = 0.944, g0 = 0.3, sigma = 30, noccasions = c(3,4), groups = c('F', 'M'), detectfn = "HR")
male <- ciScenarios$group == 'M'
ciScenarios$D[male] <- 1.847
ciScenarios$g0[male] <- 0.4
ciScenarios$sigma[male] <- 35
ciScenarios[,1:8]
ciGrid <- read.traps("ciTraps.txt")
ciMask <- make.mask(ciGrid, buffer = 160, type = "trapbuffer", poly = ciArea)
exfn <- function(x) {
  if (inherits(x, 'secr') & !is.null(x$fit)) {
    pred <- predict(x)
    pred[[1]][c('D','pmix'),]
  }
  else data.frame()
}
rawCiScen <- run.scenarios(100, ciScenarios, trapset = list(ciGrid))
summary(rawCiScen)
ciSims <- run.scenarios(100, scenarios = ciScenarios, trapset = list(ciGrid), fit = TRUE,
                        extractfn = exfn, fit.args = list(model = list(g0~h2, sigma~h2),
                        hcov = 'group'), maskset = list(ciMask))
summary(select.stats(ciSims, 'D'))$OUTPUT
summary(select.stats(ciSims, 'pmix'))$OUTPUT

And here is my OUTPUT

> summary(rawCiScen)
run.scenarios(nrepl = 100, scenarios = ciScenarios, trapset = list(ciGrid))

Replicates    100
Started       09:54:27 09 Nov 2017
Run time      1.034  minutes
Output class  selectedstatistics

$constant
            value
trapsindex      1
nrepeats        1
detectfn        1
recapfactor     1
popindex        1
detindex        1
fitindex        1

$varying
 scenario group noccasions     D  g0 sigma maskindex
        1     F          3 0.944 0.3    30         1
        1     M          3 1.847 0.4    35         2
        2     F          4 0.944 0.3    30         1
        2     M          4 1.847 0.4    35         2

$detectors
 trapsindex trapsname
          1    traps1

OUTPUT

$1
 group = F, noccasions = 3, D = 0.944, g0 = 0.3, sigma = 30, maskindex = 1 +
 group = M, noccasions = 3, D = 1.847, g0 = 0.4, sigma = 35, maskindex = 2
          n      mean      se
F.n    100 160.37000 1.14094
F.ndet 100 340.38000 2.63162
F.nmov 100 175.13000 1.59732
F.dpa  100   2.08314 0.00548
M.n    100 323.63000 1.92714
M.ndet 100 796.44000 4.84900
M.nmov 100 460.78000 2.99834
M.dpa  100   2.40987 0.00391

$2
 group = F, noccasions = 4, D = 0.944, g0 = 0.3, sigma = 30, maskindex = 1 +
 group = M, noccasions = 4, D = 1.847, g0 = 0.4, sigma = 35, maskindex = 2
          n       mean      se
F.n    100  166.07000 1.35790
F.ndet 100  457.97000 3.82341
F.nmov 100  283.97000 2.54104
F.dpa  100    2.68110 0.00787
M.n    100  323.13000 1.89532
M.ndet 100 1054.04000 6.11283
M.nmov 100  711.56000 4.22846
M.dpa  100    3.15587 0.00474

Looks good but then:

> ciSims <- run.scenarios(100, scenarios = ciScenarios, trapset = list(ciGrid), fit = TRUE,
+                         extractfn = exfn, fit.args = list(model = list(g0~h2, sigma~h2),
+                         hcov = 'group'), maskset = list(ciMask))
Completed scenario 1
Completed scenario 2
Completed in 14.951 minutes
> summary(select.stats(ciSims, 'D'))$OUTPUT
$`1`
            n mean se
estimate    0   NA NA
SE.estimate 0   NA NA
lcl         0   NA NA
ucl         0   NA NA
RB          0   NA NA
RSE         0   NA NA
COV         0   NA NA

$`2`
            n mean se
estimate    0   NA NA
SE.estimate 0   NA NA
lcl         0   NA NA
ucl         0   NA NA
RB          0   NA NA
RSE         0   NA NA
COV         0   NA NA

I have attached my detector layout. I have tried it without specifying 'poly' in make.mask too but no luck. 

[Rebecca G]

I updated my packages and now it runs but I get this warning and more NAs.

> ciSims <- run.scenarios(100, scenarios = ciScenarios, trapset = list(ciGrid), fit = TRUE,
+                         extractfn = exfn, fit.args = list(model = list(g0~h2, sigma~h2),
+                         hcov = 'group'), maskset = list(ciMask))
Completed scenario 1
Completed scenario 2

Completed in 7.717 minutes
There were 50 or more warnings (use warnings() to see the first 50)
> warnings()
Warning messages:
1: In start[parindx[[i]][1]] <- getdefault(names(parindx)[i]) :
  number of items to replace is not a multiple of replacement length
2: In start[parindx[[i]][1]] <- getdefault(names(parindx)[i]) :
  number of items to replace is not a multiple of replacement length

[Murray Efford]

I'm checking this out. It doesn't work for me, either.

Murray

--
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+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Murray Efford]

Hi again

Were you really needing to use detectfn = 'HR'? That is the source of the trouble.
1. When you simulate data with detectfn='HR' but do not specify the third parameter z, sim.capthist (called by run.scenarios) defaults to z = 1. With your other settings that results in a very fat tail, truncated at the mask buffer radius.
2. The fat tail results in very large numbers of detections - probably more than you really expect.
3. As you do not specify the detectfn to fit, secr.fit (again called by run.scenarios) defaults to HN. However the "true" HR value of sigma (30 m) is a terrible starting value for the HN sigma needed to model the flat, wide detections from the HR model with z=1, and the likelihood computation bombs.

I suggest keeping things simple (detectfn='HN' throughout). The alternative is to grapple with the various issues (specifying detectfn in fit.args etc.). Using HR (or HHR) is always risky.

The other warning message can probably be ignored - I'll clean it up in the next release.
Murray