# Web appendix 1: Utility functions
#
###################################

# 1. Function to create a matrix of initial values for latent variable z
ch.init <- function(ch, f){
   for (i in 1:dim(ch)[1]){ch[i,1:f[i]] <- NA}
   return(ch)
   }


# 2. Functions to create matrices with information about known latent state z
# 2.1. For capture-recapture data
known.state.cjs <- function(ch){
   state <- ch
   for (i in 1:dim(ch)[1]){
      n1 <- min(which(ch[i,]==1))
      n2 <- max(which(ch[i,]==1))
      state[i,n1:n2] <- 1
      state[i,n1] <- NA
      }
   state[state==0] <- NA
   return(state)
   }

# 2.2. For mark-recovery data
known.state.mr <- function(mr){
   state <- matrix(NA, nrow = dim(mr)[1], ncol = dim(mr)[2])
   rec <- which(rowSums(mr)==2)
   for (i in 1:length(rec)){
      n1 <- min(which(mr[rec[i],]==1))
      n2 <- max(which(mr[rec[i],]==1))
      state[rec[i],n1:n2] <- 1
      state[rec[i],n1] <- NA
      state[rec[i],n2:dim(mr)[2]] <- 0
      }
   return(state)
   }

# 2.3. For multistate data
known.state.ms <- function(ms, notseen){
   # notseen: label for “not seen”
   state <- ms
   state[state==notseen] <- NA
   for (i in 1:dim(ms)[1]){
      m <- min(which(!is.na(state[i,])))
      state[i,m] <- NA
      }
   return(state)
   }


# 3. Functions to create matrices of initial values for latent variable z. To be used if the functions 2.1, 2.2 or 2.3 are applied
# 3.1. For capture-recapture data
cjs.init.z <- function(ch,f){
   for (i in 1:dim(ch)[1]){
      if (sum(ch[i,])==1) next
      n2 <- max(which(ch[i,]==1))
      ch[i,f[i]:n2] <- NA
      }
   for (i in 1:dim(ch)[1]){
   ch[i,1:f[i]] <- NA
   }
   return(ch)
   }

# 3.2. For mark-recovery data
mr.init.z <- function(mr){
   ch <- matrix(NA, nrow = dim(mr)[1], ncol = dim(mr)[2])
   rec <- which(rowSums(mr)==1)
   for (i in 1:length(rec)){
      n1 <- which(mr[rec[i],]==1)
      ch[rec[i],n1:dim(mr)[2]] <- 0
      ch[rec[i],n1] <- NA
      }
   return(ch)
   }

# 3.3. For multistate data
ms.init.z <- function(ch, f){
   for (i in 1:dim(ch)[1]){ch[i,1:f[i]] <- NA}
   states <- max(ch, na.rm = TRUE)
   known.states <- 1:(states-1)
   v <- which(ch==states)
   ch[-v] <- NA
   ch[v] <- sample(known.states, length(v), replace = TRUE)
   return(ch)
   }


# 4. Function to simulate capture-recapture data
simul.cjs <- function(PHI, P, marked){
   n.occasions <- dim(PHI)[2] + 1
   CH <- matrix(0, ncol = n.occasions, nrow = sum(marked))
   # Define a vector with the occasion of marking
   mark.occ <- rep(1:length(marked), marked[1:length(marked)])
   # Fill the CH matrix
   for (i in 1:sum(marked)){
      CH[i, mark.occ[i]] <- 1       # Write an 1 at the release occasion
      if (mark.occ[i]==n.occasions) next
      for (t in (mark.occ[i]+1):n.occasions){
         # Bernoulli trial: has individual survived occasion?
         sur <- rbinom(1, 1, PHI[i,t-1])
         if (sur==0) break		# If dead, move to next individual 
         # Bernoulli trial: has individual been recaptured? 
         rp <- rbinom(1, 1, P[i,t-1])
         if (rp==1) CH[i,t] <- 1
         } #t
      } #i
   return(CH)
   }


# 5. Function to simulate mark-recovery data
simul.mr <- function(S, R, marked){
   n.occasions <- dim(S)[2]
   MR <- matrix(NA, ncol = n.occasions+1, nrow = sum(marked))
   # Define a vector with the occasion of marking
   mark.occ <- rep(1:n.occasions, marked)
   # Fill the CH matrix
   for (i in 1:sum(marked)){
      MR[i, mark.occ[i]] <- 1    # Write an 1 at the release occasion
      for (t in mark.occ[i]:n.occasions){
         # Bernoulli trial: has individual survived occasion? 
         sur <- rbinom(1, 1, S[i,t])
         if (sur==1) next    # If still alive, move to next occasion 
         # Bernoulli trial: has dead individual been recovered? 
         rp <- rbinom(1, 1, R[i,t])
         if (rp==0){
            MR[i,t+1] <- 0
            break
            }
         if (rp==1){
            MR[i,t+1] <- 1
            break
            }
      } #t
   } #i
      # Replace the NA in the file by 0
      MR[which(is.na(MR))] <- 0
      return(MR)
   }


# 6. Function to simulate multistate capture-recapture data
simul.ms <- function(PSI.STATE, PSI.OBS, marked, unobservable = NA){
   # Unobservable: number of state that is unobservable
   n.occasions <- dim(PSI.STATE)[4] + 1
   CH <- CH.TRUE <- matrix(NA, ncol = n.occasions, nrow = sum(marked))
   # Define a vector with the occasion of marking
   mark.occ <- matrix(0, ncol = dim(PSI.STATE)[1], nrow = sum(marked))
   g <- colSums(marked)
   for (s in 1:dim(PSI.STATE)[1]){
      if (g[s]==0) next  # To avoid error message if nothing to replace
      mark.occ[(cumsum(g[1:s])-g[s]+1)[s]:cumsum(g[1:s])[s],s] <-
      rep(1:n.occasions, marked[1:n.occasions,s])
      } #s
   for (i in 1:sum(marked)){
      for (s in 1:dim(PSI.STATE)[1]){
         if (mark.occ[i,s]==0) next
         first <- mark.occ[i,s]
         CH[i,first] <- s
         CH.TRUE[i,first] <- s
         } #s
      for (t in (first+1):n.occasions){
         # Multinomial trials for state transitions
         if (first==n.occasions) next
         state <- which(rmultinom(1, 1, PSI.STATE[CH.TRUE[i,t-1],,i,t-1])==1)
         CH.TRUE[i,t] <- state
         # Multinomial trials for observation process
         event <- which(rmultinom(1, 1, PSI.OBS[CH.TRUE[i,t],,i,t-1])==1)
         CH[i,t] <- event
         } #t
      } #i
   # Replace the NA and the highest state number (dead) in the file by 0
   CH[is.na(CH)] <- 0
   CH[CH==dim(PSI.STATE)[1]] <- 0
   CH[CH==unobservable] <- 0
   id <- numeric(0)
   for (i in 1:dim(CH)[1]){
      z <- min(which(CH[i,]!=0))
      ifelse(z==dim(CH)[2], id <- c(id,i), id <- c(id))
      }
   return(list(CH=CH[-id,], CH.TRUE=CH.TRUE[-id,]))
   # CH: capture histories to be used
   # CH.TRUE: capture histories with perfect observation
   }


# 7. Function to simulate capture-recapture data for JS analysis
simul.js <- function(PHI, P, b, N){
   B <- rmultinom(1, N, b) # Generate no. of entering ind. per occasion
   n.occasions <- dim(PHI)[2] + 1
   CH.sur <- CH.p <- matrix(0, ncol = n.occasions, nrow = N)
   # Define a vector with the occasion of entering the population
   ent.occ <- numeric()
   for (t in 1:n.occasions){
      ent.occ <- c(ent.occ, rep(t, B[t]))
      }
   # Modeling survival
   for (i in 1:N){
      CH.sur[i, ent.occ[i]] <- 1   # Write 1 when ind. enters the pop.
      if (ent.occ[i] == n.occasions) next
      for (t in (ent.occ[i]+1):n.occasions){
         # Bernoulli trial: has individual survived occasion?
         sur <- rbinom(1, 1, PHI[i,t-1])
         ifelse (sur==1, CH.sur[i,t] <- 1, break)
         } #t
      } #i
   # Modeling capture
   for (i in 1:N){
      CH.p[i,] <- rbinom(n.occasions, 1, P[i,])
      } #i
   # Full capture-recapture matrix
   CH <- CH.sur * CH.p
   
   # Remove individuals never captured
   cap.sum <- rowSums(CH)
   never <- which(cap.sum == 0)
   CH <- CH[-never,]
   Nt <- colSums(CH.sur)    # Actual population size
   return(list(CH=CH, B=B, N=Nt))
   }


# 8. Function to create a m-array based on capture-recapture data (CH)
marray <- function(CH){
   nind <- dim(CH)[1]
   n.occasions <- dim(CH)[2]
   m.array <- matrix(data = 0, ncol = n.occasions+1, nrow = n.occasions)
   # Calculate the number of released individuals at each time period
   for (t in 1:n.occasions){
      m.array[t,1] <- sum(CH[,t])
      }
   for (i in 1:nind){
      pos <- which(CH[i,]!=0)
      g <- length(pos)
      for (z in 1:(g-1)){
         m.array[pos[z],pos[z+1]] <- m.array[pos[z],pos[z+1]] + 1
         } #z
      } #i
   # Calculate the number of individuals that is never recaptured
   for (t in 1:n.occasions){
      m.array[t,n.occasions+1] <- m.array[t,1] - sum(m.array[t,2:n.occasions])
      }
   out <- m.array[1:(n.occasions-1),2:(n.occasions+1)]
   return(out)
   }


# 9. Function to create a m-array based on mark-recovery data (MR)
marray.dead <- function(MR){
   nind <- dim(MR)[1]
   n.occasions <- dim(MR)[2]
   m.array <- matrix(data = 0, ncol = n.occasions+1, nrow = n.occasions)
   # Create vector with occasion of marking 
   get.first <- function(x) min(which(x!=0))
   f <- apply(MR, 1, get.first)
   # Calculate the number of released individuals at each time period
   first <- as.numeric(table(f))
   for (t in 1:n.occasions){
      m.array[t,1] <- first[t]
      }
   # Fill m-array with recovered individuals
   rec.ind <- which(apply(MR, 1, sum)==2)
   rec <- numeric()
   for (i in 1:length(rec.ind)){
      d <- which(MR[rec.ind[i],(f[rec.ind[i]]+1):n.occasions]==1)
      rec[i] <- d + f[rec.ind[i]]
      m.array[f[rec.ind[i]],rec[i]] <- m.array[f[rec.ind[i]],rec[i]] + 1
      }
   # Calculate the number of individuals that is never recovered
   for (t in 1:n.occasions){
      m.array[t,n.occasions+1] <- m.array[t,1]-sum(m.array[t,2:n.occasions])
      }
   out <- m.array[1:(n.occasions-1),2:(n.occasions+1)]
   return(out)
   }


# 10. Function to monitor BUGS running time (BRT)
# Function written by Mike Meredith
bugs.run.time <- function(look.in=c(tempdir(), getwd())){
   for(dir in look.in){
      logfile <- file.path(dir, "log.txt")
      logtime <- file.info(logfile)$mtime  # Returns NA if file not found.
      if(!is.na(logtime)) break
   }
   if(!is.na(logtime)){
      # Get run time:
      runtime <- round(difftime(logtime,  
            file.info(file.path(dir,"script.txt"))$mtime, units="min"), 2)
      # Get name of model file:
      scriptfile <- file.path(dir, "script.txt")
      line2 <- readLines(scriptfile, 2)[2]
      model <- strsplit(basename(line2), "'")[[1]][1]
      # Display results:
      if(runtime < 0) {
         cat("WinBUGS run for model \"", model,
            "\" terminated abnormally\n", sep="")
         runtime <- NA
      } else {
         cat("WinBUGS run for model \"", model,
             "\" completed at ", strftime(logtime), "\n", sep="")
         cat('Posterior computed in', runtime, 'mins\n')
      }
   } else {
      cat("'log.txt' and 'script.txt' not found.\n")
      runtime <- NA
   }
   return(invisible(runtime))
}


# 11. Code fragments for some figures and calculations in chapter 11
# 11.1. Code for Fig. 11-7
# Needed output: imp.hoopoe
par(mfrow = c(2, 2), cex.axis = 1.2, cex.lab = 1.2, mar = c(5, 6, 1.5, 2), las = 1)
lower <- upper <- numeric()
year <- 2002:2010
for (i in 1:nyears){
   lower[i] <- quantile(ipm.hoopoe$sims.list$Ntot[,i], 0.025)
   upper[i] <- quantile(ipm.hoopoe$sims.list$Ntot[,i], 0.975)}
m1 <- min(c(ipm.hoopoe$mean$Ntot, popcount, lower), na.rm = T)
m2 <- max(c(ipm.hoopoe$mean$Ntot, popcount, upper), na.rm = T)
plot(0, 0, ylim = c(0, m2), xlim = c(1, nyears), ylab = "Population size", xlab = " ", col = "black", type = "l",  axes = F, frame = F)
axis(2)
axis(1, at = 1:nyears, labels = year)
polygon(x = c(1:nyears, nyears:1), y = c(lower, upper[nyears:1]), col = "grey90", border = "grey90")
points(popcount, type = "l", col = "grey30", lwd = 2)
points(ipm.hoopoe$mean$Ntot, type = "l", col = "blue", lwd = 2)
legend(x = 2, y = 25, legend = c("Counts", "Estimates"), lty = c(1, 1),lwd = c(2, 2), col = c("grey30", "blue"), bty = "n", cex = 1)

lower <- upper <- numeric()
T <- nyears-1
for (t in 1:T){
   lower[t] <- quantile(ipm.hoopoe$sims.list$phij[,t], 0.025)
   upper[t] <- quantile(ipm.hoopoe$sims.list$phij[,t], 0.975)}
par(mgp=c(3.8,1,0))
plot(y = ipm.hoopoe$mean$phij, x = (1:T)+0.5, xlim= c(1, 9), type = "b", pch = 16, ylim = c(0, 0.6), ylab = "Annual survival probability", xlab = "", axes = F, cex = 1.5, frame = F, lwd = 2)
axis(2)
axis(1, at = 1:(T+1), labels = 2002:2010)
segments((1:T)+0.5, lower, (1:T)+0.5, upper, lwd = 2)
segments(1, ipm.hoopoe$mean$mphij, T+1, ipm.hoopoe$mean$mphij, lty = 2, lwd = 2, col = "red")
segments(1, quantile(ipm.hoopoe$sims.list$mphij, 0.025), T+1, quantile(ipm.hoopoe$sims.list$mphij, 0.025), lty = 2, col = "red") 
segments(1, quantile(ipm.hoopoe$sims.list$mphij, 0.975), T+1, quantile(ipm.hoopoe$sims.list$mphij, 0.975), lty = 2, col = "red")
for (t in 1:T){
   lower[t] <- quantile(ipm.hoopoe$sims.list$phia[,t], 0.025)
   upper[t] <- quantile(ipm.hoopoe$sims.list$phia[,t], 0.975)}
points(y=ipm.hoopoe$mean$phia, x = (1:T)+0.5, type = "b", pch = 1, cex = 1.5, lwd = 2)
segments((1:T)+0.5, lower, (1:T)+0.5, upper, lwd = 2)
segments(1, ipm.hoopoe$mean$mphia, T+1, ipm.hoopoe$mean$mphia, lty = 2, lwd = 2, col = "red")
segments(1, quantile(ipm.hoopoe$sims.list$mphia, 0.025), T+1, quantile(ipm.hoopoe$sims.list$mphia, 0.025), lty = 2, col = "red") 
segments(1, quantile(ipm.hoopoe$sims.list$mphia, 0.975), T+1, quantile(ipm.hoopoe$sims.list$mphia, 0.975), lty = 2, col = "red")
legend(x = 4.5, y = 0.66, legend = c("Adults", "Juveniles"), pch = c(1, 16), bty = "n")

lower <- upper <- numeric()
T <- nyears-1
for (t in 1:T){
   lower[t] <- quantile(ipm.hoopoe$sims.list$f[,t], 0.025)
   upper[t] <- quantile(ipm.hoopoe$sims.list$f[,t], 0.975)}
plot(y=ipm.hoopoe$mean$f, x = (1:T), type = "b", pch = 16, ylim = c(5, 9), ylab = "Fecundity (fledgling / female)", xlab = "", axes = F, cex = 1.5, frame = F, lwd = 2)
axis(2)
axis(1, at = 1:T, labels = 2003:2010)
segments((1:T), lower, (1:T), upper)
segments(1, ipm.hoopoe$mean$mfec, T, ipm.hoopoe$mean$mfec, lty = 2, lwd = 2, col = "red")
segments(1, quantile(ipm.hoopoe$sims.list$mfec, 0.025), T, quantile(ipm.hoopoe$sims.list$mfec, 0.025), lty = 2, col = "red") 
segments(1, quantile(ipm.hoopoe$sims.list$mfec, 0.975), T, quantile(ipm.hoopoe$sims.list$mfec, 0.975), lty = 2, col = "red")

lower <- upper <- numeric()
T <- nyears-1
for (t in 1:T){
   lower[t] <- quantile(ipm.hoopoe$sims.list$omega[,t], 0.025)
   upper[t] <- quantile(ipm.hoopoe$sims.list$omega[,t], 0.975)}
plot(y = ipm.hoopoe$mean$omega, x = (1:T)+0.5, xlim = c(1, 9), type = "b", pch = 16, ylim = c(0, 1.1), ylab = "Immigration rate", xlab = "", axes = F, cex = 1.5, frame = F, lwd = 2)
axis(2)
axis(1, at = 1:(T+1), labels = 2002:2010)
segments((1:T)+0.5, lower, (1:T)+0.5, upper)
segments(1, ipm.hoopoe$mean$mim, T+1, ipm.hoopoe$mean$mim, lty = 2, lwd = 2, col = "red")
segments(1, quantile(ipm.hoopoe$sims.list$mim, 0.025), T+1, quantile(ipm.hoopoe$sims.list$mim, 0.025), lty = 2, col = "red") 
segments(1, quantile(ipm.hoopoe$sims.list$mim, 0.975), T+1, quantile(ipm.hoopoe$sims.list$mim, 0.975), lty = 2, col = "red")


# 11.2. Code for Fig. 11-8 and some descriptive statistics
# Needed output: ipm.hoopoe
nyears <- 9
lambda.h <- lam.lower.h <- lam.upper.h <- numeric()
Fitted.h <- lower.h <- upper.h <- matrix(NA, nrow = nyears-1, ncol = 4)

for (i in 1:(nyears-1)){
   lambda.h[i] <- mean(ipm.hoopoe$sims.list$lambda[,i])
   lam.lower.h[i] <- quantile(ipm.hoopoe$sims.list$lambda[,i], 0.025)
   lam.upper.h[i] <- quantile(ipm.hoopoe$sims.list$lambda[,i], 0.975)
   }

for (i in 1:(nyears-1)){
   Fitted.h[i,1] <- mean(ipm.hoopoe$sims.list$phij[,i])
   lower.h[i,1] <- quantile(ipm.hoopoe$sims.list$phij[,i], 0.025)
   upper.h[i,1] <- quantile(ipm.hoopoe$sims.list$phij[,i], 0.975)
   }

for (i in 1:(nyears-1)){
   Fitted.h[i,2] <- mean(ipm.hoopoe$sims.list$phia[,i])
   lower.h[i,2] <- quantile(ipm.hoopoe$sims.list$phia[,i], 0.025)
   upper.h[i,2] <- quantile(ipm.hoopoe$sims.list$phia[,i], 0.975)
   }

for (i in 1:(nyears-1)){
   Fitted.h[i,3] <- mean(ipm.hoopoe$sims.list$f[,i])
   lower.h[i,3] <- quantile(ipm.hoopoe$sims.list$f[,i], 0.025)
   upper.h[i,3] <- quantile(ipm.hoopoe$sims.list$f[,i], 0.975)
   }

for (i in 1:(nyears-1)){
   Fitted.h[i,4] <- mean(ipm.hoopoe$sims.list$omega[,i])
   lower.h[i,4] <- quantile(ipm.hoopoe$sims.list$omega[,i], 0.025)
   upper.h[i,4] <- quantile(ipm.hoopoe$sims.list$omega[,i], 0.975)
   }


# Calculate some correlation coefficients
correl.h <- matrix(NA, ncol = 4, nrow = 5001)
for (i in 1:5001){
   correl.h[i,1] <- cor(ipm.hoopoe$sims.list$lambda[i,], ipm.hoopoe$sims.list$phij[i,])
   correl.h[i,2] <- cor(ipm.hoopoe$sims.list$lambda[i,], ipm.hoopoe$sims.list$phia[i,])
   correl.h[i,3] <- cor(ipm.hoopoe$sims.list$lambda[i,], ipm.hoopoe$sims.list$f[i,])
   correl.h[i,4] <- cor(ipm.hoopoe$sims.list$lambda[i,], ipm.hoopoe$sims.list$omega[i,])
   }

# Credible intervals of correlation coefficients
quantile(correl.h[,1], c(0.05, 0.5, 0.95), na.rm = TRUE)
quantile(correl.h[,2], c(0.05, 0.5, 0.95), na.rm = TRUE)
quantile(correl.h[,3], c(0.05, 0.5, 0.95), na.rm = TRUE)
quantile(correl.h[,4], c(0.05, 0.5, 0.95), na.rm = TRUE)


# Compute the posterior modes of correlation coefficients
m <- density(correl.h[,1], na.rm = TRUE)
m$x[which(m$y==max(m$y))]

m <- density(correl.h[,2], na.rm = TRUE)
m$x[which(m$y==max(m$y))]

m <- density(correl.h[,3], na.rm = TRUE)
m$x[which(m$y==max(m$y))]

m <- density(correl.h[,4], na.rm = TRUE)
m$x[which(m$y==max(m$y))]


# Probability that correlation coefficients (r) > 0
sum(correl.h[!is.na(correl.h[,1]),1]>0)/5001
sum(correl.h[!is.na(correl.h[,2]),2]>0)/5001
sum(correl.h[!is.na(correl.h[,3]),3]>0)/5001
sum(correl.h[!is.na(correl.h[,4]),4]>0)/5001


# Plot Fig. 11-8
par(mfrow = c(2, 2), mar = c(5.5, 4, 1.5, 1), mgp=c(3, 1, 0), las = 1, cex = 1)
linecol <- c("grey70")
plot(y = lambda.h, Fitted.h[,1], type = "n", xlim = c(0.05, 0.25), ylim = c(0.6, 1.8), ylab = "Population growth rate", xlab = "Juvenile survival", frame = FALSE, pch = 19)
segments(Fitted.h[,1], lam.lower.h, Fitted.h[,1], lam.upper.h, col = linecol)
segments(lower.h[,1], lambda.h, upper.h[,1], lambda.h, col = linecol)
points(y = lambda.h, Fitted.h[,1], pch = 19, col = "black")
text(x = 0.13, y = 0.75, "r = 0.75 (0.11, 0.90)", pos = 4, font = 3, cex = 0.8)
text(x = 0.13, y = 0.65, "P(r>0) = 0.97", pos = 4, font = 3, cex = 0.8)

par(mar = c( 5.5, 4, 1.5, 1))
plot(y = lambda.h, Fitted.h[,2], type = "n", xlim = c(0.3, 0.51), ylim = c(0.6, 1.8),  ylab = "", xlab = "Adult survival", frame.plot = FALSE, pch = 19)
segments(Fitted.h[,2], lam.lower.h, Fitted.h[,2], lam.upper.h, col = linecol)
segments(lower.h[,2], lambda.h, upper.h[,2], lambda.h, col = linecol)
points(y = lambda.h, Fitted.h[,2], pch = 19, col = "black")
text(x = 0.295, y = 1.75, "r = 0.26 (-0.44, 0.73)", pos = 4, font = 3, cex = 0.8)
text(x = 0.295, y = 1.65, "P(r>0) = 0.71", pos = 4, font = 3, cex = 0.8)

par(mar = c(5, 4, 2, 1))
plot(y = lambda.h, Fitted.h[,3], type = "n", xlim = c(5, 8.5), ylim = c(0.6, 1.8), ylab = "Population growth rate", xlab = "Fecundity", frame.plot = FALSE, pch = 19)
segments(Fitted.h[,3], lam.lower.h, Fitted.h[,3], lam.upper.h, col = linecol)
segments(lower.h[,3], lambda.h, upper.h[,3], lambda.h, col = linecol)
points(y=lambda.h, Fitted.h[,3], pch = 19, col = "black")
text(x = 5, y = 1.75, "r = 0.70 (0.21, 0.87)", pos = 4,  font = 3, cex = 0.8)
text(x = 5, y = 1.65, "P(r>0) = 0.99", pos = 4, font = 3, cex = 0.8)

par(mar = c( 5, 4, 2, 1))
plot(y = lambda.h, Fitted.h[,4], type = "n", xlim = c(0, 0.8), ylim = c(0.6, 1.8),  ylab = "", xlab = "Immigration rate", frame.plot = FALSE, pch = 19)
segments(Fitted.h[,4], lam.lower.h, Fitted.h[,4], lam.upper.h, col = linecol)
segments(lower.h[,4], lambda.h, upper.h[,4], lambda.h, col = linecol)
points(y=lambda.h, Fitted.h[,4], pch = 19, col = "black")
text(x = 0.35, y = 0.8, "r = 0.83 (-0.29, 0.94)", pos = 4, font = 3, cex = 0.8)
text(x = 0.35, y = 0.7, "P(r>0) = 0.86", pos = 4, font = 3, cex = 0.8)