# The Barbary lion model. The data is from pooling expert opinions across different factors, and different experts. See paper below to # see how the distribution was derived. # Lee, T. E., et al. (2015) "Assessing uncertainty in sighting records: an example of the Barbary lion" # Input: random numbers drawn from the individual distributions for the quality of each sighting # Output: The probability that the lion is extinct in each year after the last 'certain' sighting. #============================================================================================================================================== library(rjags) library(R2jags) Years = seq(1926,2016) Tall <- 122 t0<-1895 N <- 91 # The model runs for every year after the last certain sighting (N years) n <- Tall-N qDist = read.csv("Individual_q_Dist_SumW.csv", header=FALSE); # Random numbers drawn from the distribution for the quality of individual sightings RR = dim(qDist)[1]; Randi = sample(1:RR, Tall,replace=T); yyall <- c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) # Certain sightings yall <- matrix(0,1,Tall); yall[1:length(yyall)] = yyall; zzall <- c(0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0) # Uncertain sightings zall <- matrix(0,1,Tall); zall[1:length(zzall)] = zzall; # When an uncertain sighting occurs, the quality is the mean of the quality distribution. # Years following this sighting have a quality drawn randomly from the quality distribution of the most recent uncertain sighting. Thisqq <- c(1,1,1,1, 1,2,3,3, 3,3,4,4, 4,4,4,5, 6,7,7,7, 7,7,8,8, 8,9,9,10, 10,10,10,10, 10,10,11,12, 12,12,12,13, 14,14,14,14, 15,15,15,16, 17,17,17,17, 17,17,18,18, 18,18,18,18, 18,19,19,19, 19,19,19,19, 19,19,19); Thisq <- matrix(19,1,Tall); Thisq[1:length(Thisqq)] = Thisqq; qzall = matrix(nrow=Tall, ncol=1) for (k in 1:Tall) { qzall[k]<-qDist[Randi[k],Thisq[k]]; if (zall[k]==1){qzall[k] = mean(qDist[,Thisq[k]])} } #============================================================================================================================================== Extantmean = matrix(nrow=N, ncol=1) Extantsd = matrix(nrow=N, ncol=1) temean = matrix(nrow=N, ncol=1) tesd = matrix(nrow=N, ncol=1) for (k in 1:N) { y = yall[0:(n+k)] z = zall[0:(n+k)] qz = qzall[0:(n+k)] T = n+k model1.string <-" model { # set hyper priors as in Congdon BSM p170. Exponential and Jeffreys non-informative nu11~dexp(1) nu12~dexp(1) nu21~dexp(1) nu22~dexp(1) mu11~dbeta(0.5,0.5) mu12~dbeta(0.5,0.5) mu21~dbeta(0.5,0.5) mu22~dbeta(0.5,0.5) # parameterise parameters for prior of beta parameters as in Stroud(1994) # hay11<-nu11*mu11 # hby11<-nu11*(1-mu11) # hay12<-nu12*mu12 # hby12<-nu12*(1-mu12) haz21<-nu21*mu21 hbz21<-nu21*(1-mu21) haz22<-nu22*mu22 hbz22<-nu22*(1-mu22) # then program from Lee, T. E.,et al. (2013). Inferring extinctions from sighting records of variable reliability. Journal of Applied Ecology, 51 251-258 te ~ dunif(0,T) m1 ~ dunif(0,100) f1 <-0 # no false sightings in certain sighting record m2 ~ dunif(0,100) f2 ~ dunif(0,100) #Jeff~dbeta(0.5,0.5) #Extant ~ dbern(Jeff) # Jeff trick Extant ~ dbern(0.9) for (i in 1:T) { mm1[i] <- Extant * m1+(1-Extant)*step(te-i) * m1 + f1 p1[i] <- 1- exp(-mm1[i]) y[i] ~ dbern(p1[i]) mm2[i] <- Extant * m2 + (1-Extant)*step(te-i) * m2 + f2 p2[i] <- 1 - exp(-mm2[i]) z[i] ~ dbern(p2[i]) # for first couplet truth table all [i] get same hay or hby. if declared currently extant you are extant and get 1st chunk params; if declared # currently non-extant then depending upon current Te this is overruled and get other chunk params . Otherwise f2 trick. # for second couplet truth table all [i] get same haz or hbz as logic above az[i]<-(Extant*haz21)+(1-Extant)*(step(te-i)*haz21)+ haz22 bz[i]<-(Extant*hbz21)+(1-Extant)*(step(te-i)*hbz21)+ hbz22 # state what stochastic model for quality observation [i=1..n] is - note beta is conjugate for bernouilli so later Bayes Factor # multiplication or joint MCMC easy. It needs two parameters scale and shape # quality between 0 and 1, beta dist. 2 parameters qz[i]~dbeta(az[i],bz[i]) # track mean or each beta ... diagnostic... fitted quality zalive[i] <- (haz21+haz22)/((haz21+hbz21+haz22+hbz22)) znotional[i] <- haz22/(haz22+hbz22) } } " model1.spec<-textConnection(model1.string) inits <- function() {list(nu11=1,nu12=1,nu21=1,nu22=1,mu11=0.5,mu12=0.5,mu21=0.5,mu22=0.5,te=10,m1=10,m2=10,f2=10,Extant=1)} #list(te = 10, # m1 = runif(1, 0, 100), # m2 = 10, # f2 = 10, # Extant = 1)} jags <- jags.model(model1.spec, data = list('y' = y,'z' = z, 'qz'=qz,'T' = T), inits=inits, n.chains=1) message(k, sprintf(" out of %s\n", N)) res <- coda.samples(jags, c('Extant', 'te'), n.iter=60000, n.burnin= 30000) #class(res) #summary(res) rmat <- as.matrix(res) Extantmean[k] = mean(rmat[,1]) Extantsd[k] = sd(rmat[,1]) temean[k] = mean(rmat[,2]) tesd[k] = sd(rmat[,2]) } Adf = data.frame(Extantmean,Extantsd, temean,tesd) col_headings<-c('Extant_mean','Extant_sd', 'te_mean','te_sd') names(Adf)<-col_headings write.csv(Adf, "Lion_QualityChangePoint_.csv") #========== plotting================================ png('Lion_QualityChangePoint_.png') First05 = min(which(Extantmean <0.5)); Year1 = t0+round(temean[First05]); Year2 = t0+round(temean[N]); Extinctmean = 1- Extantmean; plot(Years, Extinctmean,type="n", ylim=range(c(0,1)), pch=19, xlab="Years", ylab="Prob. extinct",cex.lab=1.4, xaxs="i",yaxs="i", # xaxs=i is to flush frame to axis limits # main="Scatter plot with std.dev error bars" ) lines(Years, Extinctmean, type="l") # line for mean # hack: draw arrows but with very special "arrowheads" arrows(Years, Extinctmean-Extantsd, Years, Extinctmean+Extantsd, length=0.05, angle=90, code=3) # errorbars for (k in 1:N){ if (zall[n+k]==1) {points(Years[k],Extinctmean[k],pch = 19, bg = "black", cex = 1, lwd = 2)} } polygon(c(Years[1],Year1,Year1,Years[1]), c(1,1,0,0),col=rgb(0.5, 0.5, 0.5,0.25), border=NA) polygon(c(Years[1],Year2,Year2,Years[1]), c(1,1,0,0),col=rgb(0.5, 0.5, 0.5,0.25), border=NA) dev.off()