variance of Lincoln-Petersen

questions concerning analysis/theory using program MARK

variance of Lincoln-Petersen

Postby markmiller » Sat Dec 30, 2017 7:23 pm

I cannot locate the formula for the variance of the Lincoln-Petersen Method. Every source I check provides the variance of the Chapman estimator. I think, but am not sure, that the variance of the Lincoln-Petersen Method can be expressed in R syntax as:

Code: Select all
var <- (n1)*(n2)*(n1-m)*(n2-m)/(m)^2/(m)


where n1, n2, and m are as conventionally defined:

n1 = number captured on first visit
n2 = number captured on second visit
m = number captured on both visits

If I try using the Huggins model in MARK with one primary period and two secondary periods the SE based on the proposed variance formula above matched the SE returned by MARK. So, I think the above formula is correct. However, I would feel much more comfortable if someone verified that formula above or if I could find that formula in a paper or book.

I suspect MARK is using maximum likelihood to obtain the variance. I have not yet tried programming a multinomial version of the Lincoln-Petersen into optim in R to see whether I can obtain the same SE that MARK returns.

Could someone email me a copy of Chapman 1951, Lincoln 1930 (and maybe Petersen (1896)? I cannot locate those publications either. I suspect one or more of them includes the above proposed formula for the variance.

Chapman, D.G. (1951). Some properties of the hypergeometric distribution with applications to zoological sample censuses.

Lincoln, F. C. (1930). Calculating Waterfowl Abundance on the Basis of Banding Returns. United States Department of Agriculture Circular. 118: 1–4.

Petersen, C. G. J. (1896). The Yearly Immigration of Young Plaice Into the Limfjord From the German Sea, Report of the Danish Biological Station (1895), 6, 5–84.

If I get the R code in optim to work and return a SE that matches MARK I will post it here.

Thank you sincerely,

Mark

mark.wayne.miller@gmail.com
markmiller
 
Posts: 49
Joined: Fri Nov 08, 2013 6:23 pm

Re: variance of Lincoln-Petersen

Postby markmiller » Wed Jan 03, 2018 8:07 pm

So far I have boiled the problem of the Lincoln-Petersen down to a 4-sided die problem. This code estimates N = 79.49601 when the true value of N = 80. I have also written code that uses a multinomial logit link (not shown) which also estimates N = 79.496. I have not yet tried Lagrange multipliers to see whether they might help.

I have put in a request to Inter-library loan to try to obtain the three papers mentioned above.

Code: Select all
N  <- 80     # number of trials
pA <- 0.30   # probability of obtaining State A
pB <- 0.45   # probability of obtaining State B
pC <- 0.10   # probability of obtaining State C
pD <- 0.15   # probability of obtaining State D
nA <- N * pA # number of trials resulting in State A
nB <- N * pB # number of trials resulting in State B
nC <- N * pC # number of trials resulting in State C
nD <- N * pD # number of trials resulting in State D

# check that log-likelihood is correct
llh <- (  log(factorial(N)) -
         (log(factorial(nA)) + log(factorial(nB)) + log(factorial(nC)) + log(factorial(N - (nA+nB+nC)))) +
         (nA               * log(pA) +
          nB               * log(pB) +
          nC               * log(pC) +
          (N - (nA+nB+nC)) * log(1 - (pA+pB+pC)))  )
llh <- -1 * llh

# check that likelihood is correct
llh.2 <- (  (factorial(N) / (factorial(nA) * factorial(nB) * factorial(nC) * factorial(N - (nA+nB+nC)))) *
             pA^nA * pB^nB * pC^nC * (1 - (pA+pB+pC)) ^(N - (nA+nB+nC)))

# identify number of trials, N, that corresponds to minimum values of -llh
my.data <- data.frame(N = N, llh = llh, llh.2 = llh.2)
my.data[my.data$llh   == min(my.data$llh),   ]
#   N      llh       llh.2
#1 80 6.250862 0.001928791
my.data[my.data$llh.2 == min(my.data$llh.2), ]
#   N      llh       llh.2
#1 80 6.250862 0.001928791

# estimate with log-likelihood
my.function <- function(betas, nA, nB, nC){
     N  = betas[1]
     llh <- (  log(factorial(N)) -
              (log(factorial(nA)) + log(factorial(nB)) + log(factorial(nC)) + log(factorial(N - (nA+nB+nC)))) +
              (nA               * log(pA) +
               nB               * log(pB) +
               nC               * log(pC) +
               (N - (nA+nB+nC)) * log(1 - (pA+pB+pC)))  )
     -1 * llh
}

my.output <- optim(c((nA+nB+nC)), my.function, nA=nA, nB=nB, nC=nC, method="Brent", lower = 0, upper = 200, hessian=TRUE)
my.output$par
#[1] 79.49601

# estimate with likelihood
my.function <- function(betas, nA, nB, nC){
     N  = betas[1]
     llh <- (  (factorial(N) / (factorial(nA) * factorial(nB) * factorial(nC) * factorial(N - (nA+nB+nC)))) *
               pA^nA * pB^nB * pC^nC * (1 - (pA+pB+pC)) ^(N - (nA+nB+nC)))
     -1 * llh
}

my.output3 <- optim(c((nA+nB+nC)), my.function, nA=nA, nB=nB, nC=nC, method="Brent", lower = 0, upper = 200, hessian=TRUE)
my.output3$par
#[1] 79.49601
markmiller
 
Posts: 49
Joined: Fri Nov 08, 2013 6:23 pm

Re: variance of Lincoln-Petersen

Postby markmiller » Fri Jan 05, 2018 5:47 am

I now realize that the estimated N = 79.496 in my previous post matches the estimate MARK returns when using the unconditional Robust Design with one primary period and two secondary periods after I constrain p1 = 0.75 and p2 = 0.40 (their true values).

From here it is straight-forward to obtain the MLE of the variance of N using the R code I posted above and to generalize that R model to also estimate p1 and p2.

I imagine I also will be able to write similar code that does not include N in the likelihood to match the Huggins model, but I have not tried that yet.

Once I locate the variance formula from my first post in the literature I will provide that reference.
markmiller
 
Posts: 49
Joined: Fri Nov 08, 2013 6:23 pm


Return to analysis help

Who is online

Users browsing this forum: No registered users and 10 guests