AIC for models with and without group effects in secr

[murray.efford]

An off-list conversation prompts me to post this warning. Care is needed when comparing models with and without an effect of group (e.g. sex).

In 'secr.fit', group effects are specified in two stages: grouping is defined with the 'groups' argument, and the term 'g' may be used in models for D, g0, sigma etc.

Defining groups fundamentally changes the structure of the likelihood, so AIC should not be compared between models fitted with different 'groups' arguments, or with and without groups defined. To evaluate the effect of adding a group effect to a model, compare models with and without 'g', making sure to define 'groups' in both cases.

For example:

Code: Select all

library(secr)
tempgrid <- make.grid()
tempCH <- sim.capthist(tempgrid)
dummysex <- sample(size=nrow(tempCH), c('F','M'), replace = T, prob=c(0.5,0.5))
covariates(tempCH) <- data.frame(sex=dummysex)

fit.0 <- secr.fit (tempCH, model = list(g0 ~ 1, sigma ~ 1), groups = 'sex', trace = F)
fit.sex <- secr.fit (tempCH, model = list(g0 ~ g, sigma ~ g), groups = 'sex', trace = F)
AIC (fit.0, fit.sex)


This should give a result something like

Code: Select all

                   model   detectfn npar    logLik     AIC    AICc dAICc  AICwt
fit.0   D~1 g0~1 sigma~1 halfnormal    3 -229.0988 464.198 465.698 0.000 0.9654
fit.sex D~1 g0~g sigma~g halfnormal    5 -229.0356 468.071 472.357 6.659 0.0346


Do NOT do this:

Code: Select all

fit.1 <- secr.fit (tempCH, model = list(g0 ~ 1, sigma ~ 1), trace = F)
AIC(fit.1, fit.sex)


As an aside - the issue does not arise when 'sex' is used as a categorical individual covariate in a model fitted by maximising the conditional likelihood (CL = T).

Murray

[mobear]

Sorry, somehow submitted a post before it was complete, please see below.

[mobear]

I am using secr 2.5.0 in Rstudio on Windows 7.

I am trying to estimate density of a black bear population using DNA-based encounter history data from 5 hair snare arrays with Maximum Likelihood. Arrays were 210 km2 and divided into 81, 2.6 km2 cells, each with one hair snare (ie, proximity detector; n = 403). The spatial arrangement of the array is surprisingly similar to Efford et al. 2005 with brushtail possums except I am using solid grids rather than hollow.

This is my first ‘secr’ analysis so I’m trying to grasp an understanding of ‘secr’ by generating models for just one of my arrays.

I have successfully ran a number of basic models using combinations of the provided effects on detection parameters (b, t, h2, k), as well as models with sex as a group-level covariate (g).

Only 4 of my models include this group effect, so following this warning I included the groups=’Sex’ argument in all of my models regardless of whether I used ‘g’ in my model.

Curious if including versus excluding the ‘groups’ argument altered my real parameter estimates, I ran my full suite of models both with and without the argument (obviously sparing the models using ‘g’ as an effect).

As I expected, summing the density estimates of each group for a model with the ‘groups’ argument equaled the density estimate for the same model without the ‘groups’ argument…

…But with two exceptions:

1)

>A.grp <- secr.fit(bsex, model=list(D ~ 1, g0 ~ h2, sigma ~ h2), groups = "Sex", method = "BFGS", mask = maskb, trace = F)

Density result:

_________D____SE__lcl___ucl
>A.grp.g1 0.13 0.04 0.06 0.24
>A.grp.g2 0.13 0.04 0.06 0.24
>A.sum12 0.25 0.09 0.13 0.48

>A.nogrp <- secr.fit(bsex, model=list(D ~ 1, g0 ~ h2, sigma ~ h2), method = "BFGS", mask = maskb, trace=F)

Density result:

_________D____SE__lcl___ucl
>A.nogrp 1.45 0.89 0.48 4.40

2)

>B.grp <- secr.fit(bsex, model=list(D ~ 1, g0 ~ b*h2, sigma ~ h2), groups = "Sex", method = "BFGS", mask = maskb, trace = F)

Density result:

_________D____SE__lcl___ucl
>B.grp.g1 0.39 0.20 0.15 0.99
>B.grp.g2 0.39 0.20 0.15 0.99
>B.sum12 0.78 0.39 0.31 1.97

>B.nogrp <- secr.fit(bsex, model=list(D ~ 1, g0 ~ b*h2, sigma ~ h2), method = "BFGS", mask = maskb, trace=F)

Density result:

_________D____SE__lcl___ucl
>B.nogrp 0.13 0.03 0.08 0.20

My question is, why don’t the summed group densities match with the no group densities for only these two models, all else being equal? Does it have something to do with the influence of h2 on sigma and g0? The large difference in densities using the same models with/without groups is concerning and the 1.45 and 0.78 bears/km2 densities are illogically high.

Apologies if this is an oversight in reading the manual or a really novice ML question, but thank you for any help. Also forgive my terrible tables, obviously haven't taken the time yet to learn the insert code and image functions.

[murray.efford]

Hi
I'm on holiday just now so can't reply fully... Although you seem to have been careful & systematic about this, I would caution that h2 models are not to be used lightly, and combining them with groups and learned responses may be part of the problem (and I'm not ruling out a bug). I expect you have only modest numbers of animals & recaptures, and any benefits of tuning the detection model are likely to get lost in the wide confidence intervals. There may be a specific problem with the numerical maximisation that I could troubleshoot if you send that data subset offline ('save' capthist and mask, as RData file) for me to check next week. Sorry not to be more helpful at this point.
Murray

[mobear]

Thank you very much for the quick reply and hope you enjoyed your holiday.

Yes, I have very modest captures and recaptures. The array I am using here to learn is actually the one with the most captures/recaptures, although still lots of 0’s.

I didn’t actually combine both ‘g’ and h2 as an effect on g0 or sigma. But maybe I am misunderstanding how specifying groups=’Sex’ influences the detection model, even when ‘g’ is not included as an effect.

These are the four models where I actually used sex as an effect:

bsex1 <- secr.fit(bsex, model=list(D ~ 1, g0 ~ g, sigma ~ 1), groups = "Sex", method = "BFGS", mask = maskb, trace = F)
bsex2 <- secr.fit(bsex, model=list(D ~ 1, g0 ~ 1, sigma ~ g), groups = "Sex", ……)
bsex3 <- secr.fit(bsex, model=list(D ~ 1, g0 ~ g, sigma ~ g), groups = "Sex", ……)
bsex4 <- secr.fit(bsex, model=list(D ~ g, g0 ~ g, sigma ~ g), groups = "Sex", ……)

Would it be more efficacious to only choose sex or h2 to include in my suite of models? I included both because I read a bear paper (Drewry et al. 2012) that found greater support for h2 over sex so I wanted to test if my data supported this as well.

I came across a question you wrote in the FAQ section of your ‘secr’ manual – “Do you really need to fit that complex model?” Given the sparseness of my data, perhaps this applies to some or all of my h2 models, which I included to try to account for realistic bear behavior. My primary goal is to obtain density estimates of the greatest precision possible, so I don’t have any particular attachment to using both h2 and sex regarding a specific hypothesis.

Perhaps this is too philosophical for the analysis board, but thought I would throw it out there.

I sent my capthist and mask files to murray.efford@otago.ac.nz so you can troubleshoot.

Many thanks for your help!

[murray.efford]

After some investigation I've concluded that Clay's problem arises from difficulties in maximizing the likelihood of SECR finite mixture models. It's well known (e.g. Brooks et al. 1997, Pledger 2000) that in finite mixture models the likelihood can have multiple modes, with an ever-present risk that the numerical maximisation algorithm will get stuck on a local maximum. In Clay's case flipping between grouped and ungrouped data was enough to make a difference (perhaps because of the way starting values are computed). I haven't come across the problem before (I don't use the models much), and it is certainly data-dependent (in Clay's data there are a few very prominent individuals).

When in doubt, try specifying starting values. You can use the alternate model if you like. For instance, from Clays first example:

Code: Select all

secr.fit (bsex, start = A.grp, etc.

A more extensive investigation is needed to provide general advice on how to recognise and resolve the problem. Any takers? One approach is to compute the profile likelihood across values of the mixing proportion (Brooks et al. 1997 Biometrics), but this is slow (an overnight job, probably; see code below). My first trial suggests fitting more groups is not a solution (contra Pledger 2000). It is better to avoid h2 and h3 models. The next release of secr will have a new mechanism for 2-class data where at least some individuals are of known class (e.g. sex).
Murray

Some R code to compute & plot profile likelihood for varying mixture proportion in h2 model (data in bsex capthist object, mask in maskb).

Code: Select all

library(secr)
pmvals <- seq(0.01, 0.99, 0.01)
npm <- length(pmvals)
out <- vector('list', npm)
for (pm in 1:length(pmvals)) {
    ## fit model with fixed beta value for mixing proportion
    out[[pm]] <- secr.fit(bsex, model = list(g0~h2, sigma~h2),
                          mask = maskb, trace = F,
                          details=list(fixedbeta=c(NA,NA,NA,NA,NA,logit(pmvals[pm]))))
    save(out, file='pmix.profile.RData')
    cat('Completed ', pmvals[pm], '\n')
    flush.console()
}

## later, or in a different R session
load(file = 'pmix.profile.RData')
out <- out[!sapply(out, is.null)]  ## in case not all finished

par(mfrow=c(1,2))

plot(pmvals[1:length(out)], sapply(out, logLik),
     xlim = c(0,1),
     xlab = 'Fixed pmix', ylab = 'Profile log-likelihood')

plot(pmvals[1:length(out)], 100*sapply(out, function(x) predict(x)[[1]]['D','estimate']),
     xlim = c(0,1),
     xlab = 'Fixed pmix', ylab = 'Density estimate / km^2')