Monday, October 28, 2024

Visual illustration of Lincoln-Petersen capture-mark-recapture method for estimating population size

Today in an intermediate-level ecology/evolution class I was explaining the Lincoln-Petersen method for estimating population size from a capture-mark-recapture study.

As I worked through the concepts, I realized a visual illustration could be useful – so I made one.

Here it is. If you’re scrutinizing the code, you may notice that I used a loop to sleep the simulation on certain frames (rather than a long sleep). This is so that each plot would be rendered so I could export them & convert the set to a GIF.

lincoln_petersen<-function(n,K,N=200,...){
  if(hasArg(sleep.long)) sleep.long<-list(...)$sleep.long
  else sleep.long<-1.2
  if(hasArg(sleep.short)) sleep.short<-list(...)$sleep.short
  else sleep.short<-0.04
  n.sleep<-ceiling(sleep.long/sleep.short)
  X<-cbind(runif(n=N),runif(n=N))
  par(mar=c(1.1,1.1,2.1,1.1))
  for(i in 1:n.sleep){
    dev.hold()
    plot(NA,asp=1,xlab="",ylab="",bty="n",axes=FALSE,
      xlim=c(-0.25,1.25),ylim=c(0,1))
    mtext("(a) unmarked population",line=0)
    polygon(c(0,1,1,0),c(0,0,1,1),col="lightgreen")
    points(X,pch=21,bg=phytools::make.transparent("blue",0.75),
      cex=1.5)
    legend(x=1.05,y=1,c("unmarked","captured or\nre-captured",
      "marked"),cex=0.7,pt.cex=c(1.5,2.5,1.5),pch=c(21,1,21),
      pt.bg=c("blue",NA,"red"),lwd=c(1,2,1),bty="n",lty=NA)
    legend(x=-0.25,y=1,c("n =","K = ","k = ",
      expression(paste(hat(N),' ='))),bty="n",cex=0.8)
    Sys.sleep(sleep.short)
    dev.flush()
  }
  ii<-sample(1:nrow(X),n)
  for(i in 1:n.sleep){
    dev.hold()
    plot(NA,asp=1,xlab="",ylab="",bty="n",axes=FALSE,
      xlim=c(-0.25,1.25),ylim=c(0,1))
    mtext("(b) capture sample for marking",line=0)
    polygon(c(0,1,1,0),c(0,0,1,1),col="lightgreen")
    points(X,pch=21,bg=phytools::make.transparent("blue",0.75),
      cex=1.5)
    legend(x=1.05,y=1,c("unmarked","captured or\nre-captured",
      "marked"),cex=0.7,pt.cex=c(1.5,2.5,1.5),pch=c(21,1,21),
      pt.bg=c("blue",NA,"red"),lwd=c(1,2,1),bty="n",lty=NA)
    legend(x=-0.25,y=1,c("n =","K = ","k = ",
      expression(paste(hat(N),' ='))),bty="n",cex=0.8)
    points(X[ii,],cex=2.5,lwd=2)
    dev.flush()
    Sys.sleep(sleep.short)
  }
  for(i in 1:n.sleep){
    plot(NA,asp=1,xlab="",ylab="",bty="n",axes=FALSE,
      xlim=c(-0.25,1.25),ylim=c(0,1))
    mtext("(c) mark all captured individuals",line=0)
    polygon(c(0,1,1,0),c(0,0,1,1),col="lightgreen")
    points(X[-ii,],pch=21,bg=phytools::make.transparent("blue",0.75),
      cex=1.5)
    points(X[ii,],pch=21,bg=phytools::make.transparent("red",0.75),
      cex=1.5)
    legend(x=1.05,y=1,c("unmarked","captured or\nre-captured",
      "marked"),cex=0.7,pt.cex=c(1.5,2.5,1.5),pch=c(21,1,21),
      pt.bg=c("blue",NA,"red"),lwd=c(1,2,1),bty="n",lty=NA)
    legend(x=-0.25,y=1,c(paste("n =",n),"K = ","k = ",
      expression(paste(hat(N),' ='))),bty="n",cex=0.8)
    dev.flush()
    Sys.sleep(sleep.short)
  }
  for(i in 1:100){
    for(j in 1:N){
      X[j,]<-X[j,]+rnorm(n=2,sd=0.02)
      while(any(X[j,]>1)||any(X[j,]<0)){
        if(X[j,1]>1) X[j,1]<-2-X[j,1]
        if(X[j,2]>1) X[j,2]<-2-X[j,2]
        if(X[j,1]<0) X[j,1]<--X[j,1]
        if(X[j,2]<0) X[j,2]<--X[j,2]
      }
    }
    dev.hold()
    plot(NA,asp=1,xlab="",ylab="",bty="n",axes=FALSE,
      xlim=c(-0.25,1.25),ylim=c(0,1))
    mtext("(d) allow population to mix freely",line=0)
    polygon(c(0,1,1,0),c(0,0,1,1),col="lightgreen")
    points(X[-ii,],pch=21,bg=phytools::make.transparent("blue",0.75),
      cex=1.5)
    points(X[ii,],pch=21,bg=phytools::make.transparent("red",0.75),
      cex=1.5)
    legend(x=1.05,y=1,c("unmarked","captured or\nre-captured",
      "marked"),cex=0.7,pt.cex=c(1.5,2.5,1.5),pch=c(21,1,21),
      pt.bg=c("blue",NA,"red"),lwd=c(1,2,1),bty="n",lty=NA)
    legend(x=-0.25,y=1,c(paste("n =",n),"K = ","k = ",
      expression(paste(hat(N),' ='))),bty="n",cex=0.8)
    dev.flush()
    Sys.sleep(sleep.short)
  }
  jj<-sample(1:nrow(X),K)
  k<-length(intersect(ii,jj))
  for(i in 1:n.sleep){
    dev.hold()
    plot(NA,asp=1,xlab="",ylab="",bty="n",axes=FALSE,
      xlim=c(-0.25,1.25),ylim=c(0,1))
    mtext("(e) re-capture sample",line=0)
    polygon(c(0,1,1,0),c(0,0,1,1),col="lightgreen")
    points(X[-ii,],pch=21,bg=phytools::make.transparent("blue",0.75),
      cex=1.5)
    points(X[ii,],pch=21,bg=phytools::make.transparent("red",0.75),
      cex=1.5)
    points(X[jj,],cex=2.5,lwd=2)
    legend(x=1.05,y=1,c("unmarked","captured or\nre-captured",
      "marked"),cex=0.7,pt.cex=c(1.5,2.5,1.5),pch=c(21,1,21),
      pt.bg=c("blue",NA,"red"),lwd=c(1,2,1),bty="n",lty=NA)
    legend(x=-0.25,y=1,c(paste("n =",n),paste("K = ",K),
      "k = ",bquote(paste(hat(N),' = '))),
      bty="n",cex=0.8)
    Sys.sleep(sleep.short)
    dev.flush()
  }
  for(i in 1:n.sleep){
    dev.hold()
    plot(NA,asp=1,xlab="",ylab="",bty="n",axes=FALSE,
      xlim=c(-0.25,1.25),ylim=c(0,1))
    mtext("(f) count number marked in first sample & estimate N",
      line=0)
    polygon(c(0,1,1,0),c(0,0,1,1),col="lightgreen")
    points(X[-ii,],pch=21,bg=phytools::make.transparent("blue",0.75),
      cex=1.5)
    points(X[ii,],pch=21,bg=phytools::make.transparent("red",0.75),
      cex=1.5)
    points(X[jj,],cex=2.5,lwd=2)
    points(X[intersect(ii,jj),],cex=2.5,lwd=2,
      pch=13,col="red")
    legend(x=1.05,y=1,c("unmarked","captured or\nre-captured",
      "marked"),cex=0.7,pt.cex=c(1.5,2.5,1.5),pch=c(21,1,21),
      pt.bg=c("blue",NA,"red"),lwd=c(1,2,1),bty="n",lty=NA)
    legend(x=-0.25,y=1,c(paste("n =",n),paste("K = ",K),
      paste("k = ",length(intersect(ii,jj))),
      bquote(paste(hat(N),' = ',.(round(n*K/k),1)))),
      bty="n",cex=0.8)
    Sys.sleep(sleep.short)
    dev.flush()
  }
  Sys.sleep(1)
  invisible(setNames(c(n,K,k,n*K/k),
    c("n","K","k","Nhat")))
}

OK, let’s test it out. First, I’ll imagine a true population size of 250, which I’m going to sample by marking 37 in my first sample and 44 in the second. (These numbers don’t mean anything. I’m just imagining two different random samples of similar effort.)

lincoln_petersen(n=37,K=44,250)