Questions about sex-specific estimates

[Duane R. Diefenbach]

Thanks Murray for your recent fix to region.N and suggestions for
obtaining sex-specific abundance estimates.

I have a few questions about estimating abundance with a couple of the
sex-specific models.

1) I have coded sex-specific models both using male/female as a
'session' and using 'sex' as a covariate. The advantage of the 'session'
approach is that region.N gives you sex-specific abundance estimates. Is
there a way to get a total population estimate? I suppose I could use
the delta method to combine the male and female variance estimates, but
is that appropriate if some parameters are in common between sexes (so
the covariance is ^=0; e.g., g0 ~ 1, and sigma ~ sex)? Or at least close
enough?

2) With 'sex' as a covariate I used the derived() function to obtain
sex-specific density estimates. However, it appears the CIs are much
smaller than when I use 'session'. I assume that these derived density
estimates are missing some variance component? Any suggestions? In my
data it is pretty clear that the sex ratio (M:F) is 1:2 so I suspect a
model without a sex effect is going to have smaller variance so I'd
prefer to combine the M and F estimates into a total abundance estimate.

3)  The reason I'd like to use 'sex' as a covariate is that I have
surveys from different years and would like to share some parameters
across years (like having sigma different by sex but the same across
years) and I don't see how to do that using the 'session' approach -
unless I am missing some trick...  Suggestions?

Thanks,

Duane

--
Duane R.Diefenbach, Ph.D.
Unit Leader and Adjunct Professor of Wildlife Ecology
U.S. Geological Survey
Pennsylvania Cooperative Fish and Wildlife Research Unit
Penn State University
404 Forest Resources Bldg.
University Park, PA 16802

Email: DR...@psu.edu
Phone: 814-865-3992
Fax: 814-863-4710

http://ecosystems.psu.edu/deer
http://www.coopunits.org/Pennsylvania

[Murray Efford]

Hi Duane

Some responses inserted below. Hope this moves you forward!

Murray

On Thursday, September 6, 2018 at 8:10:37 AM UTC+12, drd11 wrote:

Thanks Murray for your recent fix to region.N and suggestions for
obtaining sex-specific abundance estimates.

I have a few questions about estimating abundance with a couple of the
sex-specific models.

1) I have coded sex-specific models both using male/female as a
'session' and using 'sex' as a covariate. The advantage of the 'session'
approach is that region.N gives you sex-specific abundance estimates. Is
there a way to get a total population estimate? I suppose I could use
the delta method to combine the male and female variance estimates, but
is that appropriate if some parameters are in common between sexes (so
the covariance is ^=0; e.g., g0 ~ 1, and sigma ~ sex)? Or at least close
enough?

I don't have good intuition on the limits of the delta method, but I expect it would be good enough. Also, there is now a direct approach using contrasts that I'll demonstrate with a subset of the ovenbird data.

library(secr)
ov1 <- subset(ovenCH, session = 1)
msk <- make.mask(traps(ov1), buffer = 300, nx = 25, type = 'trapbuffer')
ofit4 <- secr.fit (ov1MF, model = list(D~session, sigma~session), mask = msk,
                   details = list(contrasts = list(session = contr.sum)))
exp(coef(ofit4)['D',-2]) * 2 # double average for 2 sexes
#        beta       lcl      ucl
# D 0.9835145 0.4710242 2.053611

See the openCR vignette Section 7.8 for explanation. I haven't bothered to back-transform the SE from the log scale (more fiddly than 'exp').

2) With 'sex' as a covariate I used the derived() function to obtain
sex-specific density estimates. However, it appears the CIs are much
smaller than when I use 'session'. I assume that these derived density
estimates are missing some variance component? Any suggestions? In my
data it is pretty clear that the sex ratio (M:F) is 1:2 so I suspect a
model without a sex effect is going to have smaller variance so I'd
prefer to combine the M and F estimates into a total abundance estimate.

I don't understand - with the ovenbird data I get fairly similar CI. Maybe your data are exposing an issue I haven't seen...

# sex as covariate
ofit <- secr.fit (ov1, CL = TRUE, sigma~Sex, mask = msk)
derived(ofit, groups = 'Sex')
# $`F`
# estimate SE.estimate       lcl      ucl       CVn       CVa       CVD
# esa 13.7693154          NA        NA       NA        NA        NA        NA
# D    0.5083768   0.3144909 0.1665909 1.551387 0.3779645 0.4897253 0.6186178
#
# $M
# estimate SE.estimate       lcl     ucl       CVn       CVa      CVD
# esa 27.3297900          NA        NA      NA        NA        NA       NA
# D    0.4756714   0.2035574 0.2129357 1.06259 0.2773501 0.3258942 0.427937

# sex as session
ov1MF <- split(ov1, covariates(ov1)$Sex)
ofit2 <- secr.fit (ov1MF, model = list(D~session, sigma~session), mask = msk)
predict(ofit2)
# $`session = F`
# link    estimate  SE.estimate         lcl          ucl
# D       log  0.50838366  0.347195208  0.15119452   1.70941349
# g0    logit  0.02399149  0.008476867  0.01194642   0.04759607
# sigma   log 60.68164425 21.615116367 30.82217364 119.46795161
#
# $`session = M`
# link    estimate  SE.estimate         lcl          ucl
# D       log  0.47567130  0.213242719  0.20561063   1.10044499
# g0    logit  0.02399149  0.008476867  0.01194642   0.04759607
# sigma   log 96.19470928 24.996740460 58.28609125 158.75866601

 

3)  The reason I'd like to use 'sex' as a covariate is that I have
surveys from different years and would like to share some parameters
across years (like having sigma different by sex but the same across
years) and I don't see how to do that using the 'session' approach -
unless I am missing some trick...  Suggestions?

If you do want to go down the path of having a 'session' for each sex and year combination then it is straightforward to code session covariates to get what you want. For example (not tested)

sesscovdf <- data.frame(year = factor(rep(2005:2009, each = 2)), sex = rep(c('F','M'), 5))
sesscovdf
#    year sex
# 1  2005   F
# 2  2005   M
# 3  2006   F
# 4  2006   M
# 5  2007   F
# 6  2007   M
# 7  2008   F
# 8  2008   M
# 9  2009   F
# 10 2009   M
ofit5 <- secr.fit (ovMF, model = list(D~year, sigma~sex), mask = msk, sessioncov = sesscovdf)

[Eric H]

Good day,

We're trying to estimate sex-specific and total black bear densities from hair snag data from multiple independent trap arrays. We'd also like to be able to estimate average density across several arrays, with standard errors and asymmetrical (back-transformed) confidence limits. However, I'm having difficulty extracting or calculating the standard error of total bear density on the link scale from the individual trap arrays.

We have already fit sets of candidate models to each data set using the conditional likelihood with sex as a covariate, and estimated density from AICc-minimizing models. We can easily calculate total and sex-specific densities, with SEs and CIs on the "real" (response) scale, using derived(). I can also use the delta method to calculate average variances and SEs of total density across across multiple traplines, but since these are on the real scale associated confidence limits are symmetrical around the average density.  I don't know how to calculate or extract an SE for density on the link scale from the output from derived(). 

I found this thread and modeled sexes as sessions in a full likelihood model to see if I could use the contrasts method Murray provided an example of.  I used the detection model from the AICc-minimizing CL model and tried both D~sex (as sessions) and D~1. The model with sex-specific density yielded the same density estimates as the CL model; the model with constant density yielded a slightly lower estimate.  When I tried to apply the contrast method to the model with sex-specific density, it reports density as half the total density from the model with constant density.  The model with sex-specific density was supported by AICc, and we expect higher female than male density.  Also there is still no SE on the link scale, so I still wouldn't be able to use the delta method to calculate average D across traplines. The contrasts method therefore doesn't seem useful in this particular case(??).

I also calculated total density (both sexes) and used the delta method to obtain the SE and CIs for total density from the fitted FL model with sex-specific density, as follows:

fit.LR=secr.fit(Dat.LR, model= list(D~session,g0~1,sigma~T+session), mask=mask.LR, details=list(contrasts=list(session=contr.sum)))
reals =predict(fit.LR,realnames='D',type="response")
total.D=log(reals[[1]][1,2]+reals[[2]][1,2]) # yields total D on link scale; backtransforming gives same estimate as derived(CL model).
#
delta.mo = function(se) {
  vr=se^2
  mean.vr=sum(vr) / length(vr)
  mean.se=sqrt(mean.vr)
  return(mean.se)
}
#
betas=predict(fit.LR, realnames='D', type="link")
mean.se = delta.mo(se=c(betas[[1]][1,3],betas[[2]][1,3]))
CI = c(mean.d - (1.96*mean.se), mean.d + (1.96*mean.se))
#
real.D.CI = c(exp(mean.d), exp(CI))
real.D.CI
# 0.0015536801  0.0006751119  0.0035755877

The resulting density estimate matches the one obtained using derived(CL model) (below), but my delta method confidence intervals are quite a bit wider.

D = 0.001553678, SE = 0.0005069814, CI = 0.0008329 - 0.00289816.

So, my approach would provide the SE around total density on the link scale needed to calculate SEs and asymmetrical confidence limits around average D across traplines, but seems to sacrifice precision.

We know that we could fit models to multisession data with sex as a covariate and trap arrays as sessions to obtain average density across the arrays. I've haven't done so (yet) because detection models with slightly different forms were supported for the different data sets. This is probably simply because the true model is not (ever) in the candidate set, and similar but not identical models were supported due to sampling variation between the traplines. However, perhaps bear behaviour actually differed among arrays, in which case we'd need a detection model with a lot of interactions with session to model detection appropriately. This seems less than ideal. We will likely try fitting multisession models with with the most commonly-supported forms of variation in detection parameters, and additive session effects (if supported), to pooled data, but I'd still like to find a solution to the problem described above if possible.

I think what would be best for us is if we could calculate or extract the SE around total density on the link scale from the CL models. Is this possible?  Then we could estimate density from the AICc-supported detection model for each trapline, and calculate average density across traplines (all used the same # of traps and covered a similar area) and SEs and asymmetrical CIs using the delta method (right?). 

[Murray Efford]

Hi Eric

I have yet to digest your whole post, but starting at the bottom...

- variances for derived D are inherently slightly ad hoc, as they use my application of the Huggins non-spatial variance that no statistician has spent time on (see the 'supplements' to our 2008 paper that I wrote, not David).

- have you noticed that the default (loginterval = TRUE) is for derived() to compute CL on the log (i.e. default link) scale? The code for this is buried in an internal function (called by derived.secr, see getAnywhere(derived.secr)) that you can access with secr:::add.cl. This uses a (delta method) estimate of the log-scale SE. You could easily adapt this.

- in case it's relevant and you forgot, there is the opportunity to define groups post-hoc in derived.secr (groups argument).

Murray

[Murray Efford]

I see secr:::derived.secr also works - you learn something every day.

[Eric H]

Thanks very much Murray,

I looked at the code for derived() and add.cl(),

It seems that obtaining the SE, on the (log) link scale, around the density estimated by applying derived() to an secr model object created by fitting the conditional likelihood is as simple as:

# where "der" is the data frame created when one applies derived() to a fitted model,
# such that der[2,1] and der[2,2] are the estimated density and it's SE on the real scale, respectively:
#
SE.derived.D= function(der) {
  se.D.link = sqrt(log(1+(der[2,2] / der[2,1])^2))
  return(se.D.link)
}

I used the SE I obtained to manually calculate confidence limits on the link scale and then back-transformed them to the real scale. The results matched those from derived(). 

I also compared the result obtained using the above to the beta coefficients obtained by fitting a full likelihood model (same detection model, sexes modeled as sessions) to the same data.

The CL model, derived, and the above function yielded an SE on the link scale of 0.1625689 (this is around the total density of both sexes, but density differs between sexes because the detection model included sex effects on g0 and sigma).

The FL model yielded SE.beta of 0.16144109 around the first density parameter (intercept), and 0.16135984 around the second density parameter (effect of session, i.e. of sex).

This should serve our purposes well.

Thanks again,

Eric Howe