# 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