Trouble estimating density by group

[sixtystrat]

I am having trouble getting secr to output estimates of density by sex. I have set up my R script as such:
DeerCH <- read.capthist("CaptQualMF.txt","JoeTrapID.txt", detector = "proximity", covnames = "sex", fmt = "trapID")
and run the following model:
secr29 <- secr.fit(DeerCH, model = list(g0~1, sigma~t+h2), groups="sex", method = "BFGS", mask = AAFBhabmask, CL = FALSE, detectfn = 0, buffer=1500, trace = TRUE)

It seems to estimate everything fine but I am having trouble figuring out how to get it to output the density estimates by sex.

Thanks!
Joe

[howeer]

Hi Joe,

I think you need to include D~sex in the model description, e.g. model = list(g0~1, sigma~t+h2, D~sex). By omitting D from the model description, you're likely getting D~1 (constant across groups) as a default.

Eric

[sixtystrat]

Thanks Eric. I'll give it a try.

[sixtystrat]

One more question... I also am trying to estimate density by habitat type (n = 3) and ran the following model: D~habtype g0~time + h2 sigma~time + h2 pmix~h2.

The output shows only the density calculation for habitat #3 and my beta estimate for habitat type #2 has a high SE ( mean = -9.9878716, SE = 200.5290845). Should I suspect a problem? Will the "derived" command give me the overall density?

[murray.efford]

...I also am trying to estimate density by habitat type (n = 3) and ran the following model: D~habtype g0~time + h2 sigma~time + h2 pmix~h2.
The output shows only the density calculation for habitat #3

The default in predict.secr is to present 'real' parameter estimates for one level only. To get others you need to specify the covariate levels. In your case something like this should work

Code: Select all

predict(fittedmodel, newdata = data.frame(habtype=factor(1:3)))


I'm assuming that you have saved the result of secr.fit as fittedmodel, that your habtype covariate is a factor with levels 1,2,3, and that 'time' is a type for 't' rather than another covariate (another covariate would have to be included in newdata)

and my beta estimate for habitat type #2 has a high SE ( mean = -9.9878716, SE = 200.5290845). Should I suspect a problem?

Yes, the wild SE is a worry and probably means the model is not fully identifiable. I haven't much experience of fitting models that combine finite mixtures (h2) with explicit accounting for other sources of variation, and don't have all the facts here, so I am wary of providing specific advice. I suppose I would try a simpler model and see if it gave essentially the same density or N estimates (we're not chained to the model with smallest AIC).
Murray

[sixtystrat]

Thanks Murray. I applied a simpler model (D=~habtype, g0 ~1, sigma~1) and got these results using the derived command:
D = 0.05467124 SE = 0.01384239,
and with (D=~1, g0 ~1, sigma~1), got
D = 0.0563542 SE = 0.01437930, essentially the same but still with the high SE on the habitat 3 beta estimate.

With a more complicated model, (D = ~ habtype, g0~ t+h2, sigma~ t+h2), the density estimate with the derived command is
D = 0.1390182 SE = 0.04017682
and the estimate without the habitat covariate is (g0~ t+h2, sigma~ t+h2) is:
D = 0.1203805 SE = 0.03741416

Would that suggest that the estimates for the habtype models are okay?
BTW, the t represents 5 capture periods and I also have models where group is defined by sex. How would I code the predict command to get estimates by sex and time as well? Thanks.
Joe

[murray.efford]

Joe
I'm reluctant to offer an opinion on your results when I can't see the whole picture, and I don't have time just now to offer more, sorry. You can include sex and time in the newdata argument of predict - something like

Code: Select all

predict(..., newdata = data.frame (sex = c('F','M'), t = factor(1:5)))

We use 'factor' to ensure t is understood as a categorical predictor (one level per occasion). By definition, for a closed population density does not vary with t, but maybe you want to look at variation in g0 and sigma. The character vector c('F','M') is automatically converted to a factor with levels in alphabetical order.

For derived (used to get density from CL models that may include individual covariates like sex) you can specify groups post hoc with e.g. derived(..., groups = 'sex')
Murray