huggins, how to calculate N

[j_p]

Hi, I am using a huggins model in RMark to analyse abundance of some point count data. The data are 10 minute point counts, split into 2 minute intervals. My goal is to predict relative abundance a species at each point count, based on the observer and two individual covariates, noise and vegetation density. I can create model averaged estimates of p, my question is how do I either a) use these to predict N for each point count or b) is there another way to predict N for each point count?
Thanks for your help!
JP

[bacollier]

Hi, I am using a huggins model in RMark to analyse abundance of some point count data. The data are 10 minute point counts, split into 2 minute intervals. My goal is to predict relative abundance a species at each point count, based on the observer and two individual covariates, noise and vegetation density. I can create model averaged estimates of p, my question is how do I either a) use these to predict N for each point count or b) is there another way to predict N for each point count?
Thanks for your help!
JP

JP,
Huggins style closed capture models are typically not focused on point-specific density estimates but rather on overall abundance of the sampled area, I don't think you will be able to estimate point density(cleanly) like you are thinking. I suppose you could brute force it and estimate a point-specific p based on your best model and the covariate data from each point location and then apply that to the count at that point, but that is fairly ad hoc and likely not trustworthy.

Bret

[jlaake]

I've been in the field with no email contact or would have replied sooner. Bret, makes a good point but if you still want to do it you have 2 options. First, you can specify each point as a group and the N estimates should be available in the derived list (model$results$derived). You can model average those values with model.average.list in RMark or use the aicmodavg(sp?) package in R. Alternatively, you can use your model.averaged values of p and construct your own estimate from the n's. How you do that depends on the model for the data. You'll also need to derive the variance estimate this way. Was this removal or independent point counts? As a side note, the "N" you get in RMark from any of the closed models is really f0 (number not caught) and must be added to number caught to get N. The values in $derived are from MARK and include f0 and number caught(seen).

--jeff

[marshbirder]

I'm trying to use Huggins as well, based on Thompson and La Sorte 2008. I have looked through the forum here for advice before posting and have found some previous responses helpful.

Right now I have one year of data, but I will eventually have three. So, for now, I just want to be able to estimate abundance without looking at changes over time. I see that you are saying here I shouldn't expect to get abundance at each point count station with Huggins, which is fine. My covariates should tease out the differences I'm looking for (which would be based on point count locations). My research is looking at songbird abundance relative to silviculture, with focal species I can use, rather than the entire community so that I can use the more abundant birds for analysis. I have two sites that are already harvested and two that will have one year of pre- and at least two years of post-harvest sampling.

*123 point count stations
*2013 data are two sites harvested and two sites pre-harvest
*3 visits per year
*10-min point counts with five 2-min intervals
*We will use only detections within 50m
*I want to relate abundance to slope, aspect, and basal area gradient (measured in 4 prism plots at each point count station) with more variables to be added later

It was suggested that I run Huggins to get p and c and then use a GLMM to model abundance relative to my linear variables. The MARK portion of this seems much simpler than what I was initially thinking I would have to do, although I'm not sure how to code for more than 3 species (0 1, 1 0, 0 0, ?). Thompson and La Sorte used 6 species in their analyses. Does this seem like a good approach or is there a model I can use that would incorporate my continuous variables?

I have only ~6 weeks of a 1-credit graduate level MARK course under my belt.

Thanks!

[ctlamb]

As a side note, the "N" you get in RMark from any of the closed models is really f0 (number not caught) and must be added to number caught to get N. The values in $derived are from MARK and include f0 and number caught(seen).

Jeff, just to follow up on this,

when I do the following:

Code: Select all

###SR MODEL
SR1<-mark(data=SR.proc,ddl=SR.ddl,model.parameters = list(p=list(formula=~Sex+TrapNights+TrapType+time,share=TRUE),
                                                         Phi=list(formula=~1),
                                                         f=list(formula=~1)),output=FALSE,model="RDPdfHuggins")

##View Derived Pop size
a<-SR1$results$derived$`N Population Size`

###Take mean of males and females across all years
PopSR<-mean(a[1:8,1])+mean(a[9:16,1])

##View MARK output
SR1  ##same as a




I receive identical abundance estimates, from RMark or from the MARK output. I interpreted these as the true population abundance, as opposed to f0. Am I correct?

Also, how am I able to get the standard deviation for the population sizes? trying to calculate CV's, but I have been taking the mean of male and female abundances for each year, and I only have SE's. I guess I would just need the sample size for each year? and can go SD=sqrt(n)*SE?

[jlaake]

The derived parameter for N is the total abundance and N=f0+M_t+1 where M_t+1 is the total number of unique individuals caught and f0 is the estimate of the number never caught.

With the new release of RMark all derived parameters are now extracted back to R in the result. They should have the se and conf limit. If you have two separate groups and you want to combine them you want to take the square root of the sum of all of the elements in the derived.vcv. Try the example below. It is a common mistake to think you need to do something with n but that relationship between std dev of the data and the std dev (std error) of the mean is only for means. For a parameter the std error is the measure of precision of the parameter and the std error^2 is the variance of the parameter.

Code: Select all

data(edwards.eberhardt)
edwards.eberhardt$sex=factor(c(rep(c("M","F"),38)))
                             
ee.huggins.m0=mark(edwards.eberhardt,groups="sex",model="Huggins",model.parameters=list(p=pdotshared))

ee.huggins.m0$results$derived
ee.huggins.m0$results$derived.vcv

#std errors for each are

sqrt(diag(ee.huggins.m0$results$derived.vcv$`N Population Size`))


# total abundance is
sum(ee.huggins.m0$results$derived$`N Population Size`$estimate)
# se for total abundance is square root of sum of variances and covariances
sqrt(sum(ee.huggins.m0$results$derived.vcv$`N Population Size`))



[ctlamb]

Thanks for the explanation, Jeff. I certainly had been making this misconception regarding n and parameter estimates, luckily I hadn't been successful in deriving any such misconceived numbers.

So, to be clear, I am looking to derive a measure of spread around my Huggins estimate. I have seen coefficients of variations used, which I understood to be the standard deviation/mean, but am I understanding correctly that the SE and SD for a parameter estimate are equivalent? Wherein the CV for the example provided above would be CV=7% based on this:

Code: Select all

> sum(ee.huggins.m0$results$derived$`N Population Size`$estimate)
[1] 97.19319
> sqrt(sum(ee.huggins.m0$results$derived.vcv$`N Population Size`))
[1] 6.97389
> sqrt(sum(ee.huggins.m0$results$derived.vcv$`N Population Size`))/sum(ee.huggins.m0$results$derived$`N Population Size`$estimate)
[1] 0.07175287



FYI I thin your pdothsared call was missing so I amended the mark call a such:

Code: Select all

ee.huggins.m0=mark(edwards.eberhardt,groups="sex",model="Huggins",model.parameters=list(p=list(formula=~1,share=TRUE)))


[cooch]

You need to be somewhat careful in constructing a valid CI for abundance, especially f you're model averaging (which you almost inevitably will be). In particular, evidence suggests that f0 follows a log-normal distribution, and this approaches to constructing the CI based on any normal asymptotics are probably wrong.

You are referred to section 14.1`0.1 in chapter 14 of the MARK book. While MARK (and thus, presumably, RMark) correctly estimates the unconditional variance, and SE, the reported UCI and LCI (which are based on normal asymptotic) are incorrect, and the correct CI for model averaged estimate of N require you do them by hand. I'm guessing that RMark doesn't do this correction, but Jeff will need to advise.

Note that this issue applies regardless of whether N is included in the likelihood, or conditioned out (i.e., Huggins models).