#  CountRemoval.r - A REMOVAL MODEL FOR ESTIMATING DETECTION PROBABILITIES
#                   FROM POINT COUNT SURVEYS
rm(list=ls())

CntRem <- function(x,npoints,point_radius=50) {
  # x            = counts in 1st, 2nd and 3rd intervals
  # npoints      = number of count locations
  # point_radius = radius of observable area around each point
  #
  # returns estimates of
  #   p-hat = prob of detection over 10-minute interval
  #   N     = number of individuals present around all points
  #   D     = number of individuals per unit area

  #  llcomb = log-likelihood combinatorial term
  llcomb=ll=lgamma(sum(x)+1)-lgamma(x[1]+1)-lgamma(x[2]+1)-lgamma(x[3]+1)


  cellprobs <- function(c,q) {  # computes cell probs for each time interval
    q10=q^(10);
    p1=(1.-c*q*q*q)/(1.-c*q10);
    p2=c*q*q*q*(1.-q*q)/(1.-c*q10);
    p3=c*q*q*q*q*q*(1.-q*q*q*q*q)/(1.-c*q10);
    return(c(p1,p2,p3))
  }

  est_qc <- function(S) {           #  estimation routine #1 - estimates
    S=plogis(S); q=S[1];  c=1;      #     q = 1 - Pr(detection in a minute)
    if (length(S)>1) c=S[2]         #     c = heterogeneity parameter (prob in grp2)
    prbs=cellprobs(c,q)
    ll=llcomb + sum(x*log(prbs))
    return(-ll)
  }
  est_phat <- function(S) {               #  estimation routine #2 - same as #1
    S=plogis(S); q10=S[1]; q=q10^.1; c=1; #  except that it is re-parameterized to
    if (length(S)>1)  c=(1-S[2])/q10;     #  estimate p-hat instead of q
    if (c>1) c=1
    prbs=cellprobs(c,q)
    ll=llcomb + sum(x*log(prbs))
    return(-ll)
  }

  cat('\nCountRemoval - input: X=',x,'Count locations:',npoints,'radius:',point_radius,'\n')

  #  use R nlm function to find minimum neg. log-likelihood estimates
  #    for two-param model (Mc)
  suppressWarnings(z1<-nlm(est_qc,rep(0,2),hessian=T));

  S=plogis(z1$estimate);                 #  get ML estimates
  aic_mc=2*z1$minimum+2*length(S)        # compute aic
  ec=try(vc<-solve(z1$hessian),silent=T) #  compute var-cov matrix from hessian
  if (is.character(ec))
    vc2=matrix(NA,2,2)
  else {
    g=diag(exp(-z1$estimate)/(1+exp(-z1$estimate))^2); vc2=g %*% vc %*% t(g)
  }
  #  compute std. errors and print
  suppressWarnings(se2 <- sqrt(diag(vc2)))
  out=cbind(S,se2); rownames(out)=c('q','c');
  colnames(out)=c('estimate','std.err');
  cat('\nmodel M(c): ll=',-z1$minimum,'aic=',aic_mc,'\n')
  print(out)

  #  use R nlm function to find minimum neg. log-likelihood estimates
  #    for one-param model (M)
  suppressWarnings(z1a<-nlm(est_qc,1,hessian=T));

  S=plogis(z1a$estimate);             # get ML estimates
  aic_m0=2*z1a$minimum+2*length(S)    # compute aic
  #     compute var-cov matrix from hessian
  vc=1/z1a$hessian; g=exp(-z1a$estimate)/(1+exp(-z1a$estimate))^2
  vc2=g * vc * g; suppressWarnings(se2<-sqrt(vc2))
  #  print estimates and std. errors
  out=matrix(c(S,1,se2,NA),2,2); rownames(out)=c('q','c');
  colnames(out)=c('estimate','std.err');
  cat('\nModel M():ll=',-z1a$minimum,'aic=',aic_m0,'\n'); print(out)

  #   do lokelihood ratio test to compare models
  cat('\nLikelihood Ratio Test - M(c) vs M()\n')
  chisq=2*(z1a$minimum-z1$minimum); df=1; pr=1-pchisq(chisq,1); ec2=1
  cat('  Chi-square=',chisq,' df=',1,' Pr(Larger chi-sq)=',pr,'\n')

  #  Instead of using likelihood-ratio test results,
  #   we'll compare aic values from the two models
  #    and use the model with the lowest aic
  if (aic_mc<aic_m0) {
    #  model Mc had lower aic, so re-run model using
    #  p-hat parameterization to get estimate of p-hat
    suppressWarnings(z2<-nlm(est_phat,c(0,0),hessian=T))
    mname='M(c)-reparameterized'
    S=plogis(z2$estimate);

    #  Compute var-cov matrix using SURVIV 2nd partial derivitave method.
    #  This method is sometimes more accurate than the inverse
    #   of the Hessian matrix.
    ppar=S; p=cellprobs(c=(1-ppar[2])/ppar[1], q=ppar[1]^0.1)
    parpi=cbind(p,p)   #  parpi = derivatives of cell probs w.r.t params

    for (i in 1:2) {  # compute partial derivatives...
      parval=ppar[i]
      ppar[i]=ppar[i]*.9999
      p=cellprobs((1-ppar[2])/ppar[1], ppar[1]^.1)
      parpi[,i]=(p-parpi[,i])/(parval-ppar[i])
      ppar[i]=parval
    }
      # construct info matrix in array H
    h=matrix(0,2,2)
    for (i in 1:2) {
      for (j in 1:2) {
        parval=sum(parpi[,i]*parpi[,j]/p)
        h[i,j]=sum(sum(x)*parval)
      }
    }
    #  var-cov matrix is inverse of info matrix, H
    ec2=try(vc2<-solve(h))

    #  save p-hat and std.error to print later
    suppressWarnings(se2<-sqrt(diag(vc2)))
    out=cbind(S,se2); rnames=c('IGNORE','p-hat')

  } else {

    #  model M had lower aic, so re-run model using
    #  p-hat parameterization to get estimate of p-hat
    suppressWarnings(z2<-nlm(est_phat,qlogis(.7),hessian=T));
    mname='M()-reparameterized';
    S=1-plogis(z2$estimate)

    #  Compute var-cov matrix using SURVIV 2nd partial derivitave method.
    ppar=S[1]
    parpi=cellprobs(1, (1-ppar)^.1)
    parval=ppar; ppar=ppar+.00000000001
    p=cellprobs(1, (1-ppar)^.1)
    parpi=(p-parpi)/(parval-ppar)
    # construct info matrix in array H
    parval=sum(parpi*parpi/p)
    h=sum(sum(x)*parval); vc2=matrix(1/h,1,1)

    #  save p-hat and std.error to print later
    suppressWarnings(se2<-sqrt(diag(vc2)))
    out=matrix(c(S,se2),1,2); rnames='p-hat'

  }

  #    Compute N, and D estimates with std. errors and print
  cat('\nLow AIC model:',mname,'\n')
  if (is.character(ec2)) {
     cat('\n=======================\n',ec2,'======================\n');
  } else {
     l=length(S); phat=S[l]
     n=sum(x); N=n/phat; A=npoints*pi*point_radius^2/10000; D=N/A
     varN=n^2/phat^4*vc2[l,l]
     varD=n^2*vc2[l,l]/A^2/phat^4+n*(1-phat)/A^2/phat^2
     suppressWarnings(seN<-sqrt(varN)); suppressWarnings(seD<-sqrt(varD))
     out=rbind(out,matrix(c(N,D,seN,seD),2,2))
     rownames(out)=c(rnames,'N','D'); colnames(out)=c('estimate','std.err')
     print(round(out,4))
  }
}

#  enter observation data:
#          x1=detected within the first 3-minutes,
#          x2=no. first detected within the following 2-minute interval,
#          x3=no. first detected within the last interval of 5-minutes.
#          npoints = no. count locations
#          point_radius = max distance to an observation

#  Check estimates using data from Farnsworth et. al., Table 4.

# CntRem(x=c(141,29,39),npoints=155,point_radius=50)
# CntRem(x=c(152,23,29),npoints=155,point_radius=50)
# CntRem(x=c(96, 18,26),npoints=155,point_radius=50)
# CntRem(x=c(118,26,32),npoints=155,point_radius=50)
# CntRem(x=c(64,19,25),npoints=155,point_radius=50)
# CntRem(x=c(94,15,25),npoints=155,point_radius=50)
# CntRem(x=c(43,11,9),npoints=155,point_radius=50)
# CntRem(x=c(35,9,12),npoints=155,point_radius=50)