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)

Now, let’s do it again, but this time, consider that the population might be 950 individuals strong and we sample 85 in our first outing and 75 in the second.

lincoln_petersen(n=85,K=75,950)

That’s pretty much the idea….