using secrdesign with groups

[Greg Distiller]

Hi Murry

I am doing some work on SCR design with Ian where we want to explore how different approaches perform when dealing with a population that has group-specific parameters (like different sigmas for males and females leopards). I could code everything up from scratch but I am hoping I can use the secrdesign functions. But I have read that you say "Simulation of structured populations....is still experimental" and so I'm trying to verify that I can in fact use your functions. BTW I am also planning on modelling sex using groups rather than hcov so I can get group-specific density estimates. I have a few different queries related to this:

1. I can see that the simulated CHs by default include a sex covariate (I think I read it is simulated as a 2:1 ratio). If I want to have control over the sex-specific parameters I assume I should rather use groups and then essentially ignore the default sex covariate? 

2. I have written an extract function that binds together the estimates for each group. But when I try to run the summary function I can see that I get mean and se values for the 1st group but for the 2nd group there are no values for RB and COV. Is this because the function does not know what the true parameter values are for the 2nd group? Is there any way to provide this info to the function? If not I can obviously just calculate things manually. 

3. I'm also not sure things are correct for the first group because if I manually calculate say the RB from the output I am not getting a matching value. For example here is the summary output for D from a fit to 5 reps:

             n     mean      se

estimate    5  0.00042 0.00004

SE.estimate 5  0.00013 0.00001

lcl         5  0.00024 0.00002

ucl         5  0.00075 0.00007

RB          5 -0.39827 0.05299

RSE         5  0.30381 0.02029

COV         5  0.60000 0.24495

The value for RB that is shown is -0.39 but if I try to calculate it manually with the true D of 0.0005 I get (0.00042-0.0005) / 0.0005 = -0.16. Any idea what may be going on?

4. I have read about the alternative approach to improve speed by specifying method = "none" with start = TRUE but I am a bit confused by the statement that says this only works when there is a 1:1 relationship between real and beta parameters. Surely the general situation would be that there is not such a relationship? I also thought about supplying the true values to start but suspect there will be an issue with trying to supply true values for both groups. Do you think it is better to not use this alternative approach with groups?

5. Is there anything else you can think of that I should be aware of, wrt the limitations of using secrdesign for this setup with two groups?

Many thanks, and my apologies for such a long message.

Greg

[Murray Efford]

Hi Greg

This is just a partial answer, but we'll see how far it gets us...

1. The 'secr' function sim.popn() generates a population covariate called 'sex'. The default is 50:50, but that can be set (I think you got 2:1 from an example in the secrdesign vignette here). However, that is no use to you as sim.capthist() (and hence secrdesign) doesn't use population covariates.

2. Yes, 'secrdesign' struggles to keep track of the true values - a bit of a [headache]. However, select.stats() has argument 'true' that might help - I'll see if I can get that working with a custom extract function.

More to come.

I would say: I doubt you gain anything from group-specific analyses if density and detection parameters are all group-specific. D usually needs to be group-specific because ~1 implies 50:50, which is arbitrary. Group-specific sigma to me implies group-specific g0/lambda0. But maybe you have a specific need.

Murray

[Murray Efford]

Picking this up again -

3. When the scenario has a group structure, select.stats (used within summary) automatically sums the group densities for a single, overall, true density. This is right for the use case I had in mind, but not so for separate group-specific estimates. That probably explains your result.

4. I think you are citing the secrdesign vignette. In general, secr.fit() accepts a start vector specified as beta parameters, and method = "none" works just fine. I think this also works in secrdesign ::run.scenarios() (i.e. providing a full vector of beta values as the start fit arg). The 1:1 statement applies when start is specified via the shortcut 'true': probably the coder cut his losses rather than deal with the general case of constructing the vector for more complex scenarios and models: it's up to you (then need unique fitindex for each scenario).

5. I have added new functions 'estimateArray' and 'estimateSummary' to the pre-release version of secrdesign (2.8.1; on GitHub and as source secrdesign_2.8.1.tar.gz and Windows binary secrdesign_2.8.1.zip). estimateSummary() is an alternative mechanism for summarising output from run.scenarios() that deals with group structure (when coupled with extractfn = predict). Example below. 2.8.1 may be tweaked before it goes to CRAN, but estimateSummary() should stay. Let me know if you find bugs or there are gaps in the help (?estimateSummary)!

Murray

library(secrdesign)

# 2-scenario, 2-group simulation
scen8 <- make.scenarios (D = 8, g0 = 0.3, sigma = 30, noccasions = c(4,8),
    groups = c('F','M'))
   
# replace density and sigma values of males to make it interesting
male <- scen8$group == 'M'
scen8$D[male] <- 4
scen8$sigma[male] <- 40

# run a few simulations with group model
setNumThreads(7)
grid <- make.grid(8, 8, spacing = 30)
mask <- make.mask(grid, buffer = 160, type = 'trapbuffer')
sims <- run.scenarios(10, scen8, trapset = grid, fit = TRUE,
    fit.args = list(model = list(D~g, g0~1, sigma~g), groups = 'group'),
    extractfn = predict, maskset = mask)

# format by scenario and group, and pre-pend details --
options(digits = 3)
estimateSummary(sims, 'D',  format = 'data.frame', cols = c(1,2,4,6:8))

  scenario group noccasions D  g0 sigma true.D nvalid  EST seEST       RB   seRB   RSE  RMSE  rRMSE COV
1        1     F          4 8 0.3    30      8     10 8.27 0.163  0.03371 0.0204 0.120 0.559 0.0699 1.0
2        1     M          4 4 0.3    40      4     10 4.26 0.170  0.06530 0.0425 0.155 0.573 0.1431 1.0
3        2     F          8 8 0.3    30      8     10 7.95 0.435 -0.00617 0.0544 0.111 1.305 0.1632 0.8
4        2     M          8 4 0.3    40      4     10 4.23 0.155  0.05846 0.0387 0.138 0.520 0.1301 1.0

On Friday, February 24, 2023 at 12:43:52 AM UTC+13 greg.distiller wrote:

[Greg Distiller]

Hi Murray

As always,thanks for the rapid response and for the various clarifications. I will play around with the new functions when I get a chance and let you know if I have any further queries. Just one follow up query regarding the use of start values, is the order they must be provided in per group i.e. if there was group-specific D and sigma then one would provide starting values in this order: D1, lambda0, sigma1, D2, Sigma2?

I understand what you are saying about group-specific analyses when all parameters are group-specific but we want to explore how different SCR designs perform when there is a group structure i.e.  is there a sensible way to modify the OF used in GAoptim()  so that the resulting design performs well for both sexes.

Thanks again

Greg

--
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/3970ffee-49e6-4593-9161-1c1387146b31n%40googlegroups.com.

[Murray Efford]

When you provide 'start' as a vector of coefficients on the link scale ('beta parameters') the order follows the order of coefficients in a fitted model. For a group model, each group effect follows the intercept. For example: 

msk <- make.mask(traps(ovenCH[[1]]), buffer = 300, nx = 20)
fit <- secr.fit(ovenCH[[1]],  mask=msk, trace = FALSE,
    model = list(g0~g, sigma~g), groups = 'Sex')
coef(fit)

beta SE.beta lcl ucl D -0.7342255 0.3494230 -1.4190820 -0.04936908 g0 -3.4077967 0.6622468 -4.7057767 -2.10981676 g0.gM -0.4071990 0.7889274 -1.9534683 1.13907036 sigma 4.0227402 0.3303179 3.3753291 4.67015132 sigma.gM 0.5767424 0.3681590 -0.1448361 1.29832081

[Murray Efford]

Good to hear of work on designs that work well for disparate groups. For the record, algorithmic optimization has deeper problems (foreshadowed in Dupont et al. and Durbach et al. papers), but the disparate-group question is a good one and arises also in clustered and random designs.

On Wednesday, March 1, 2023 at 6:10:11 PM UTC+13 greg.distiller wrote:

[Greg Distiller]

Hi again Murray

I've been playing with the secrdesign group sims and am getting some strange results that I wanted to show you. Essentially, when I simulate two populations independently under a scenario that generates lots of data, I get unbiased performance for all parameters. But when I do the exact same thing using groups I am getting biased estimates for the detection function parameters (large positive bias in lambda0 and slight negative bias in sigma). See below:

###############################################

#simulation using 2 independent pops with same trap array

#all parameters unbiased (1% or less)

################################################

scen.testF <- make.scenarios (D = 8, lambda0 = 1.5, sigma = 30, detectfn = 'HHN', noccasions = 1)
scen.testM <- make.scenarios (D = 6, lambda0 = 1.2, sigma = 40, detectfn = 'HHN', noccasions = 1)

traps.testF <- make.grid(8,8, detector = 'count', spacing = 30)
mask.testF <- make.mask(traps.testF, buffer = 120, type = 'trapbuffer')

sims.testF.30 <- run.scenarios(100, scen.testF, trapset = traps.testF, fit = TRUE, fit.args = list(detectfn = 'HHN'), extractfn = predict, maskset = mask.testF)

sims.testM.30 <- run.scenarios(100, scen.testM, trapset = traps.testF, fit = TRUE, fit.args = list(detectfn = 'HHN'), extractfn = predict, maskset = mask.testF)

> D.ests <- select.stats(sims.testF.30, parameter = "D", statistics = c("estimate","SE.estimate","RB", "RSE", "COV","lcl","ucl"))
> L0.ests <- select.stats(sims.testF.30, parameter = "lambda0", statistics = c("estimate","SE.estimate","RB", "RSE", "COV","lcl","ucl"))
> Sig.ests <- select.stats(sims.testF.30, parameter = "sigma", statistics = c("estimate","SE.estimate","RB", "RSE", "COV","lcl","ucl"))

> summary(D.ests, dec = 6, fields = c("n", "mean", "se"), alpha = 0.05, type = "list")

OUTPUT
               n      mean       se
estimate    100  8.056524 0.082877
SE.estimate 100  0.916811 0.005037
lcl         100  6.450952 0.073358
ucl         100 10.062951 0.092633
RB          100  0.007066 0.010360
RSE         100  0.114375 0.000584
COV         100  0.970000 0.017145

summary(L0.ests, dec = 4, fields = c("n", "mean", "se"), alpha = 0.05, type = "list")

OUTPUT
               n   mean     se
estimate    100 1.5156 0.0111
SE.estimate 100 0.1114 0.0010
lcl         100 1.3126 0.0099
ucl         100 1.7501 0.0126
RB          100 0.0104 0.0074
RSE         100 0.0736 0.0005
COV         100 0.9500 0.0219

summary(Sig.ests, dec = 4, fields = c("n", "mean", "se"), alpha = 0.05, type = "list")

OUTPUT
               n    mean     se
estimate    100 29.9356 0.0950
SE.estimate 100  0.9496 0.0078
lcl         100 28.1319 0.0903
ucl         100 31.8556 0.1022
RB          100 -0.0021 0.0032
RSE         100  0.0317 0.0002
COV         100  0.9500 0.0219

D.ests2 <- select.stats(sims.testM.30, parameter = "D", statistics = c("estimate","SE.estimate","RB", "RSE", "COV","lcl","ucl"))
L0.ests2 <- select.stats(sims.testM.30, parameter = "lambda0", statistics = c("estimate","SE.estimate","RB", "RSE", "COV","lcl","ucl"))
Sig.ests2 <- select.stats(sims.testM.30, parameter = "sigma", statistics = c("estimate","SE.estimate","RB", "RSE", "COV","lcl","ucl"))

summary(D.ests2, dec = 6, fields = c("n", "mean", "se"), alpha = 0.05, type = "list")

OUTPUT
               n     mean       se
estimate    100 6.059114 0.071818
SE.estimate 100 0.719551 0.004657
lcl         100 4.805290 0.063148
ucl         100 7.641427 0.080808
RB          100 0.009852 0.011970
RSE         100 0.119529 0.000689
COV         100 0.970000 0.017145

summary(L0.ests2, dec = 4, fields = c("n", "mean", "se"), alpha = 0.05, type = "list")

OUTPUT
               n    mean     se
estimate    100  1.1998 0.0080
SE.estimate 100  0.0869 0.0008
lcl         100  1.0412 0.0073
ucl         100  1.3827 0.0090
RB          100 -0.0002 0.0067
RSE         100  0.0725 0.0006
COV         100  0.9700 0.0171

summary(Sig.ests2, dec = 3, fields = c("n", "mean", "se"), alpha = 0.05, type = "list")

OUTPUT

               n   mean    se
estimate    100 39.937 0.154
SE.estimate 100  1.355 0.012
lcl         100 37.369 0.143
ucl         100 42.683 0.169
RB          100 -0.002 0.004
RSE         100  0.034 0.000
COV         100  0.900 0.030

##############################

#simulation using groups

#D unbiased but large bias in lambda0 (over 100%)

#and negative bias in sigma

###############################
scen <- make.scenarios (D = 8, lambda0 = 1.5, sigma = 30, detectfn = 'HHN', noccasions = 1, groups = c('F','M'))

# replace density, lambda0 and sigma values of males
male <- scen$group == 'M'
scen$D[male] <- 6
scen$lambda0[male] <- 1.2
scen$sigma[male] <- 40

grid <- make.grid(8, 8, spacing = 30, detector = "count")

mask <- make.mask(grid, buffer = 160, type = 'trapbuffer')

sims <- run.scenarios(100, scen, trapset = grid, fit = TRUE, fit.args = list(detectfn = 'HHN', model = list(D~g, lambda0~g, sigma~g), groups = 'group'),

                      extractfn = predict, maskset = mask)

> estimateSummary(sims, 'D',  format = 'data.frame', cols = c(1,2))
  scenario group true.D nvalid      EST      seEST           RB        seRB       RSE      RMSE      rRMSE  COV
1        1     F      8    100 7.926621 0.07648465 -0.009172418 0.009560582 0.1154263 0.7645423 0.09556778 0.97
2        1     M      6    100 6.114034 0.07738254  0.019005690 0.012897090 0.1193914 0.7783453 0.12972422 0.94

> estimateSummary(sims, 'lambda0',  format = 'data.frame',cols = c(1,2))
  scenario group true.lambda0 nvalid      EST      seEST       RB       seRB       RSE     RMSE    rRMSE  COV
1        1     F          1.5    100 3.861306 0.07620475 1.574204 0.05080316 0.2016264 2.480055 1.653370 0.00
2        1     M          1.2    100 2.419413 0.03470464 1.016177 0.02892053 0.1447882 1.267361 1.056134 0.01

> estimateSummary(sims, 'sigma',  format = 'data.frame',cols = c(1,2))
  scenario group true.sigma nvalid      EST     seEST         RB        seRB        RSE     RMSE     rRMSE  COV
1        1     F         30    100 25.98717 0.1222425 -0.1337609 0.004074750 0.04228685 4.193110 0.1397703 0.07
2        1     M         40    100 35.66837 0.1440439 -0.1082908 0.003601097 0.03870661 4.562582 0.1140646 0.16

I expected these two approaches to essentially be equivalent and give very similar results. Any ideas as to what may be going on?

Many thanks

Greg

To view this discussion on the web, visit https://groups.google.com/d/msgid/secrgroup/55eea0fe-467d-4c76-8466-7f7e5dba0742n%40googlegroups.com.

[Murray Efford]

The problem seems to arise only with 'count' detectors, not (binary) 'proximity' detectors. After some probing I suspect  the issue goes beyond secrdesign, perhaps to the way counts are modelled (simulated as Poisson for each group separately, but fitted in secr.fit as one Poisson?). I'll have to search some more.

Murray

[Murray Efford]

A possible workaround: treat the groups as sessions. Simply specify multisession = TRUE in run.scenarios(). Each 'group' row in a scenario becomes a session.

mask32 <- make.mask(grid, buffer = 160, nx = 32, type = 'trapbuffer')  # for speed

sims4 <- run.scenarios(10, scen, trapset = grid, fit = TRUE, multisession = TRUE,
    fit.args = list(detectfn = 'HHN', model = list(D~session, lambda0~session, sigma~session)),
    extractfn = predict, maskset = mask32)

lapply(c(D='D', lambda0='lambda0', sigma='sigma'),    # everything all at once
    estimateSummary, object = sims4, form = 'data')

# $D
# scenario group true.D nvalid      EST     seEST          RB       seRB       RSE      RMSE      rRMSE COV
# 1        1     F      8     10 7.826335 0.2397127 -0.02170815 0.02996408 0.1162635 0.7398102 0.09247627   1
# 2        1     M      6     10 6.282612 0.2163749  0.04710204 0.03606248 0.1174163 0.7079777 0.11799628   1
#
# $lambda0
# scenario group true.lambda0 nvalid      EST      seEST            RB       seRB        RSE       RMSE      rRMSE COV
# 1        1     F          1.5     10 1.499969 0.03055152 -0.0000209231 0.02036768 0.07469113 0.09165457 0.06110305 1.0
# 2        1     M          1.2     10 1.244710 0.02463430  0.0372581521 0.02052858 0.07048493 0.08637477 0.07197898 0.9
#
# $sigma
# scenario group true.sigma nvalid      EST     seEST          RB        seRB        RSE      RMSE      rRMSE COV
# 1        1     F         30     10 29.85066 0.2482750 -0.00497794 0.008275832 0.03253589 0.7596486 0.02532162 1.0
# 2        1     M         40     10 39.52914 0.4907521 -0.01177146 0.012268802 0.03297868 1.5457186 0.03864296 0.9

[Greg Distiller]

OK thanks a lot Murray, that should work easily. One point of clarification, when you say:

"I suspect  the issue goes beyond secrdesign, perhaps to the way counts are modelled (simulated as Poisson for each group separately, but fitted in secr.fit as one Poisson?)"

Does that mean there is a bug within secr.fit when using groups? If group-specific detection functions are specified surely they should not be fitted as one PP?

Many thanks

Greg

To view this discussion on the web, visit https://groups.google.com/d/msgid/secrgroup/b9d96461-29f6-4b37-8658-ef0b2aeb9679n%40googlegroups.com.

[Murray Efford]

I'm investigating. Yes, I think there probably is a bug in secr.fit for count detectors with groups. PP = Poisson process? It now looks more like a coding error. More later.

Murray

[Murray Efford]

The bad news is that there was  a longstanding bug in secr.fit() for Poisson-distributed counts ('count' detectors with default binomN=0). The good news is that this arose only when fastproximity = FALSE (always the case when groups defined, but usually defaults to TRUE) and even then had little effect on density estimates. The bug was this: when computing the likelihood component for each capture history, probabilities g(d) instead of hazards h(d) were passed to an internal C++ function (pski). The effect is most noticeable for large hazards. Fitted lambda0 were overestimated and sigma underestimated, approximately balancing the effect on effective sampling area and hence D-hat.

My earlier advice about using multisession=TRUE was only correct by accident: in that case (session rather than group model) the fit defaulted to fastproximity = TRUE and gave the right answer. (Using multisession = TRUE while forcing fastproximity = FALSE gives the wrong answer).

I have fixed the bug in the working version  secr 4.5.9 on GitHub and in the source and Windows binary on the Density website.

Many thanks to Greg for spotting and reporting this. Always check the results make sense!

Murray

[Murray Efford]

A few simulations following the previous example -

sims1fit <- run.scenarios(10, scen, trapset = grid, fit = TRUE, extractfn = predict, fit.args = list(model = list(D~g, lambda0~g, sigma~g), detectfn = 'HHN', groups = 'group'), seed = 123) Completed scenario 1 Completed in 9.63 minutes lapply(c(D='D', lambda0='lambda0', sigma='sigma'), # everything all at once estimateSummary, object = sims1fit, form = 'data') $D scenario group true.D nvalid EST seEST RB seRB RSE RMSE rRMSE COV 1 1 F 8 10 8.119147 0.2793553 0.01489338 0.03491941 0.1147488 0.8464930 0.10581162 1 2 1 M 6 10 6.170429 0.1713804 0.02840489 0.02856340 0.1185942 0.5416525 0.09027541 1 $lambda0 scenario group true.lambda0 nvalid EST seEST RB seRB RSE RMSE rRMSE COV 1 1 F 1.5 10 1.477465 0.02952726 -0.015023482 0.01968484 0.07484968 0.09140332 0.06093555 1.0 2 1 M 1.2 10 1.208778 0.02850831 0.007314748 0.02375693 0.07123948 0.08597420 0.07164517 0.9 $sigma scenario group true.sigma nvalid EST seEST RB seRB RSE RMSE rRMSE COV 1 1 F 30 10 29.52809 0.1880911 -0.01573039 0.006269705 0.03161213 0.7355984 0.02451995 1.0 2 1 M 40 10 39.45292 0.5189599 -0.01367694 0.012973997 0.03289875 1.6502025 0.04125506 0.9

>

[Greg Distiller]

Thanks Murray! 

To view this discussion on the web, visit https://groups.google.com/d/msgid/secrgroup/55fcc7e9-8fc9-457a-9511-8cc6912b0e21n%40googlegroups.com.