###### # Main code for running RJMCMC algorithm for capture-recapture data # - dippers when only considering models C/C and C2/C. Created August 07 (RK) ###### # Define the function: with input parameters: # nt = number of iterations # nburn = burn-in dippercodeRJ <- function(nt,nburn){ # Define the parameter values: # ni = number of release years # nj = number of recapture years # nparam = maximum number of parameters ni = 6 nj = 6 nparam = 3 # Read in the data: data <- matrix(c( 11,2,0,0,0,0,9, 0,24,1,0,0,0,35, 0,0,34,2,0,0,42, 0,0,0,45,1,2,32, 0,0,0,0,51,0,37, 0,0,0,0,0,52,46),nrow=ni,byrow=T) # Read in the priors: # Beta priors on the parameters # Note: assume same prior for parameters independent of model alpha <- c(1,1,1) beta <- c(1,1,1) # Define prior model probabilities: # Set equal prior probability on each model: K <- 1/2 prob <- c(K,K) # Parameters for MH updates (Uniform random walk): delta <- c(0.1,0.1,0.1) # RJ proposal parameters for splitting survival rate # (mean and variance of Normal distribution): rjmu <- 0 rjsigma <- 0.01 # Set initial parameter values: param <- c(0.9,0.5,0.5) # param[1] = recapture rate # param[2] = survival rate - non-flood years # param[3] = survival rate - flood years # Set param[2] = param[3] for start in model C/C; # else chain is started in model C2/C # In order to calculate the initial model: model <- 1 if (param[2] != param[3]) { model <- 2 } # Model 1: C/C # Model 2: C2/C # sample is an array in which we put the sample from the posterior distribution. sample <- array(0, dim=c(nt, nparam+1)) sample <- data.frame(sample) # Label the columns of the array "sample" names(sample) <- c("p", "phin", "phif", "model") # Calculate log(likelihood) for initial state using a separate function "likhood": likhood <- calclikhood(ni, nj, data, param) # Set up a vector for parameter values and associated log(likelihood): output <- dim(nparam+2) # MCMC updates - MH algorithm: # Cycle through each iteration: for (t in 1:nt){ # Update the parameters in the model using function "updateparam": output <- updateparam(nparam, param, ni, nj, data, likhood, alpha, beta, delta, model) # Set parameter values and log(likelihood) value of current state to be the output from # the MH step: param <- output[1:nparam] likhood <- output[nparam+1] # RJ update: # So this updates the model within each iteration using function "updatemodel": output <- updatemodel(param, model, ni, nj, data, likhood, alpha, beta, prob, rjmu, rjsigma) # Set parameter values to be the output from the RJ step: # These are the parameter values, the model indicator and the corresponding # log(likelihood) value. param <- output[1:nparam] model <- output[nparam+1] likhood <- output[nparam+2] # Record the set of parameter values: for (i in 1:nparam) { sample[t,i] <- param[i] } # Record the model as the final column in the array: sample[t,nparam+1] <- model # Recall: # Model 1: C/C # Model 2: C2/C } # Calculate summary statistics using the function "summaryRJ": summaryRJ(nt,nburn,sample) # Output the sample from the posterior distribution: sample } ###### # This function calculates the log likelihood of capture-recapture data ###### calclikhood <- function(ni, nj, data, param){ # First we set the survival and recapture probs: # Set up the size of the array containing the survival and recapture probs: # Set up the set of cell probabilities (q), and initially set them to be all equal to zero: phi <- array(0,nj) p <- array(0,nj) q <- array(0,dim=c(ni,nj+1)) # Set the survival and recapture probs for each year of the study: # allowing for flood and non-flood years: for (i in 1:nj) { p[i] <- param[1] phi[i] <- param[2] } # For flood years: for (i in 2:3) { phi[i] <- param[3] } # Set initial value of the log(likelihood) to be zero likhood <- 0 # Calculate the Multinomial cell probabilities and corresponding contribution to the # log(likelihood): # Cycle through the number of years that there are releases: for (i in 1:ni){ # For diagonal elements: q[i,i] <- phi[i]*p[i] likhood <- likhood + data[i,i]*log(q[i,i]) # Calculate the elements above the diagonal: if (i <= (nj-1)) { for (j in (i+1):nj) { q[i,j] <- prod(phi[i:j])*prod(1-p[i:(j-1)])*p[j] likhood <- likhood + data[i,j]*log(q[i,j]) } } # Probability of an animal never being seen again q[i,nj+1] <- 1 - sum(q[i,i:nj]) likhood <- likhood + data[i,nj+1]*log(q[i,nj+1]) } # Output the log(likelihood) value: likhood } ###### # Function for updating the parameters values: ###### updateparam <- function(nparam, param, ni, nj, data, likhood, alpha, beta, delta, model){ # Cycle through each parameter in turn and propose to update using # random walk MH with Uniform proposal density: # nupdate is the number of parameters to be updated: nupdate <- 3 # If model = 1 then param[2] = param[3] (constant survival rate) if (model == 1) { nupdate <- 2 } for (i in 1:nupdate) { # Update the parameter # Keep a record of the current parameter value being updated oldparam <- param[i] # Propose a new value for the parameter using a random walk with # Uniform proposal density param[i] <- runif(1, param[i]-delta[i], param[i]+delta[i]) # Automatically reject any moves that are propose values outside [0,1] if (param[i] >= 0 & param[i] <= 1) { # For C/C model survival is constant in flood and non-flood years: if (model == 1 & i==2) { param[3] <- param[2] } # Calculate the log(acceptance probability): # Calculate the new likelihood value for the proposed move: # Calculate the numerator (num) and denominator (den) in turn: newlikhood <- calclikhood(ni, nj, data, param) # Include the likelihood and prior (Beta) terms in the acceptance probability num <- newlikhood + log(dbeta(param[i],alpha[i],beta[i])) den <- likhood + log(dbeta(oldparam,alpha[i],beta[i])) # All other prior terms (for other parameters) cancel in the acceptance probability. # Acceptance probability of MH step: A <- min(1,exp(num-den)) } else { A <- 0 } # To do the accept/reject step of the algorithm: # Simulate a random number in [0,1]: u <- runif(1) # Accept the move with probability A: if (u <= A) { # Accept the proposed move: # Update the log(likelihood) value: likhood <- newlikhood } else { # Reject proposed move so parameter stays at current value: param[i] <- oldparam if (model == 1 & i ==2) { param[3] <- oldparam } } } # Set the values to be outputted from the function to be the # set of parameter values and log(likelihood) value: output <- c(param, likhood) # Output the parameter values: output } ###### # This is the RJ step to update the model ###### updatemodel <- function(param, model, ni, nj, data, likhood, alpha, beta, prob, rjmu, rjsigma){ # We propose to update whether survival is different for # flood years compared to non-flood years # Record current value of the parameter and model indicator: oldparam <- array(0,2) oldparam[1] <- param[2] oldparam[2] <- param[3] oldmodel <- model reject <- 0 # If current model is C/C - propose model C2/C if (model == 1) { # Propose to add flood year dependence # Propose perturbation from current value using Normal proposal distribution # (recall that R uses the standard deviation in the Normal distribution). diff <- rnorm(1,rjmu,rjsigma) param[2] <- param[2] - diff param[3] <- param[3] + diff # Automatically reject if proposed parameter(s) are outside [0,1]: if (param[2] <= 0 | param[2] >= 1) { reject <- 1 } if (param[3] <= 0 | param[3] >= 1) { reject <- 1 } if (reject == 0) { # New proposed model indicator: model <- 2 # Calculate log of acceptance probability: # Consider numerator and denominator: # Calculate log(likelihood) for proposed state: newlikhood <- calclikhood(ni, nj, data, param) # Include the likelihood terms in the acceptance probability terms: num <- newlikhood den <- likhood # Prior in the numerator for new survival probabilities: for (i in 2:3) { num <- num + log(dbeta(param[i],alpha[i],beta[i])) } # Prior in the denominator for old survival probability: den <- den + log(dbeta(oldparam[1],alpha[i],beta[i])) # Proposal in the denominator for parameter "diff" den <- den + log(dnorm(diff,rjmu,sqrt(rjsigma))) # Priors on all other parameters cancel # P(m|m') = P(m'|m) = 1 # Jacobian = 2 num <- num + log(2) # Prior on the models: num <- num + log(prob[model]) den <- den + log(prob[oldmodel]) } } else { # Propose to move from model C2/C to C/C model <- 1 # Set new survval rate to be mean of flood and non-flood survival rate: for (i in 2:3) { param[i] <- (oldparam[1] + oldparam[2])/2 } # Calculate log of acceptance probability: # Consider the numerator and denominator. # Calculate log(likelihood) for proposed state: newlikhood <- calclikhood(ni, nj, data, param) # Include the likelihood terms in the acceptance probability terms: num <- newlikhood den <- likhood # Prior in the denominator for old survival probabilities: for (i in 1:2) { den <- den + log(dbeta(oldparam[i],alpha[i],beta[i])) } # Prior in the numerator for new survival probability: num <- num + log(dbeta(param[2],alpha[i],beta[i])) # Proposal in the numerator for parameter "diff" # Calculate the proposed Normal random variate for going from C2/C to C/C: diff <- 0.5*(oldparam[1] - oldparam[2]) num <- num + log(dnorm(diff,rjmu,sqrt(rjsigma))) # Priors on all other parameters cancel # P(m|m') = P(m'|m) = 1. # Jacobian = 1/2 den <- den + log(2) # Prior on the models: num <- num + log(prob[model]) den <- den + log(prob[oldmodel]) } # Calculate the acceptance probability: if (reject == 1) { A <- 0 } else { A <- min(1,exp(num-den)) } # Simulate u ~ U[0,1] u <- runif(1) # Accept the move with probability A, i.e. if u < A if (u <= A) { # Record new log(likelihood) value likhood <- newlikhood } else { # Reject proposed move so parameter and model stay at current values: param[2] <- oldparam[1] param[3] <- oldparam[2] model <- oldmodel } # Place set of parameter values, model indicator and log(likelihood) value in output output <- c(param, model, likhood) # Output the parameter values: output } ###### # This file uses the output from R-dippersRJ.R and calculates some # posterior summary statistics and plots. ###### summaryRJ <- function(nt,nburn,sample){ # Read in the number of iterations (nt) and burn-in (nburn). # To calculate the posterior probabilities for each model: # Model 1: C/C # Model 2: C2/C # Set up an array to calculate the number of times the chain # visits each model and for recording means/standard deviations: total <- array(0,2) mn <- array(0,3) std <- array(0,3) # Calculate the number of times the chain is in each model # after the burn-in period: subsample <- sample[((nburn+1):nt),] total[1] <- length(subsample[subsample$model=="1",4]) total[2] <- length(subsample[subsample$model=="2",4]) # Output the posterior model probabilities to the screen: cat("Posterior probabilities for each model: ", "\n") cat("\n") cat("Model 1: C/C: ", total[1]/(nt-nburn), "\n") cat("Model 2: C2/C ", total[2]/(nt-nburn), "\n") cat("\n") # To calculate the posterior density plots and estimates, conditional # on each model: for (i in 1:2) { mn[i] <- mean(subsample[subsample$model=="1",i]) std[i] <- sqrt(var(subsample[subsample$model=="1",i])) } cat("Posterior summary estimates for each model: ", "\n") cat("\n") cat("Model 1: C/C: ", "\n") cat("mean (SD)", "\n") cat("p: ", "\n") cat(mn[1], " (", std[1], ")", "\n") cat("\n") cat("phi: ", "\n") cat(mn[2], " (", std[2], ")", "\n") for (i in 1:3) { mn[i] <- mean(subsample[subsample$model=="2",i]) std[i] <- sqrt(var(subsample[subsample$model=="2",i])) } cat("\n") cat("\n") cat("Model 2: C2/C: ", "\n") cat("mean (SD)", "\n") cat("p: ", "\n") cat(mn[1], " (", std[1], ")", "\n") cat("\n") cat("phi_n: ", "\n") cat(mn[2], " (", std[2], ")", "\n") cat("phi_f: ", "\n") cat(mn[3], " (", std[3], ")", "\n") # To calculate the density plots and model-averaged summary estimates: # To make the graphic window show 6 plots: 3 rows and 1 column par(mfrow=c(3,1)) # Note to put this pack to a single plot use par(mfrow=c(1,1)) # Cycle through each parameter and calculate model-averaged mean and SD: for (i in 1:3) { mn[i] <- mean(sample[(nburn+1):nt,i]) std[i] <- sd(sample[(nburn+1):nt,i]) plot(density(sample[(nburn+1):nt,i])) } # Write out the summary statistics to the screen: cat("\n") cat("Posterior (model-averaged) summary estimates for parameters: ", "\n") cat("\n") cat("mean (SD)", "\n") cat("p: ", "\n") cat(mn[1], " (", std[1], ")", "\n") cat("\n") cat("phi_n: ", "\n") cat(mn[2], " (", std[2], ")", "\n") cat("phi_f: ", "\n") cat(mn[3], " (", std[3], ")", "\n") cat("\n") }