###############################################################
#
# Dynamics of a goshawk population across half a century is driven by the variation of first-year survival
#
# Michael Schaub, Volkher Looft, Floriane Plard, Jan A. C. von Rönn
#
###############################################################



###############################################################
#
# Code for fitting the integrated population models 
#
###############################################################
#
# Created: 14.12.2021
# Revised: 06.03.2024
# 
# Written by: Michael Schaub & Floriane Plard
#
###############################################################
#
# R-version: 4.2.2
#
###############################################################



#########################################
#
# Load libraries
#
#########################################

library(IPMbook)
library(nimble)
library(MCMCvis)



#########################################
#
# Set path
#
#########################################

path <- '...'
setwd(path)



#########################################
#
# Load functions
#
#########################################

source('Functions.txt')



#########################################
#
# Load data
#
#########################################

# 1. Feather capture-recapture data  (1968-2014)
dat1 <- read.csv('Capture-recapture_data.csv', header=TRUE, sep=';')

# 1.1. Females
CH  <- as.matrix(rbind(dat1[which(dat1$sex=='F'),2:48]))
age <- dat1[which(dat1$sex=='F'),49]

# Remove first captures in all birds. Now we just have survival of at least 2-years old individuals which have adult survival (originally 3y survival). 
ch.new <- rmFirst(CH)
z <- which(rowSums(ch.new)==0)

# Create capture-recapture m-array
marrF <- marray(as.matrix(ch.new[-z,]))

# 1.2. Males
CH  <- as.matrix(rbind(dat1[which(dat1$sex=='M'),2:48]))
ageM <- dat1[which(dat1$sex=='M'),49]

# Remove first captures in all birds. Now we just have survival of at least 2-years old individuals which have adult survival (originally 3y survival). 
ch.new <- rmFirst(CH)
z <- which(rowSums(ch.new)==0)

# Create capture-recapture m-array
marrM <- marray(as.matrix(ch.new[-z,]))

# Combine m-arrays
marr <- array(NA, dim=c(nrow(marrM), ncol(marrM), 2))
marr[,,1] <- marrF
marr[,,2] <- marrM



# 2. Dead-recovery data of birds from the study area (1969-2014)
dat2 <- read.csv('Dead-recovery_data.csv', header=TRUE, sep=';')

# 2.1. Females
recov <- subset(dat2, dat2$sex==2)
for (lh in 1:nrow(recov)){
  y <- which(recov[lh,1:46]==2)
  if (length(y>0)){
    recov[lh, c(y, y+1)] <- 1
  }
}
recov <- cbind(rep(0, nrow(recov)), recov)

# Create dead-recovery m-array
marrF <- marrayDead(recov[,1:47], freq=recov$n)

# 2.2. Males
recov <- subset(dat2, dat2$sex==1)
for (lh in 1:nrow(recov)){
  y <- which(recov[lh,1:46]==2)
  if (length(y>0)){
    recov[lh, c(y, y+1)] <- 1
  }
}
recov <- cbind(rep(0, nrow(recov)), recov)

# Create dead-recovery m-array
marrM <- marrayDead(recov[,1:47], freq=recov$n)

# Combine m-arrays
marrD <- array(NA, dim=c(nrow(marrM), ncol(marrM), 2))
marrD[,,1] <- marrF
marrD[,,2] <- marrM



# 3. Probability to successfully raise a brood (to have at least on fledgling)
dat3 <- read.csv('Productivity_1_data.csv', header=TRUE, sep=';')

s1 <- matrix(as.numeric(table(dat3$success, dat3$year, dat3$f_age)[2,,]), nrow=3, byrow=TRUE)
s2 <- matrix(as.numeric(table(dat3$success, dat3$year, dat3$f_age)[2,,]) + as.numeric(table(dat3$success, dat3$year, dat3$f_age)[1,,]), nrow=3, byrow=TRUE)
psucc <- array(NA, dim=c(3, 2, 47))
dimnames(psucc) <- list(c('1y', '2y', 'older'), c('successful', 'total'), 1968:2014)  # dimension names
psucc[,1,] <- s1
psucc[,2,] <- s2



# 4. Sex ratio of the chicks and productivity, given success
dat4 <- read.csv('Productivity_2_data.csv', header=TRUE, sep=';')

dat4 <- subset(dat4, dat4$chick_age>9)         # Only use broods when the age of the chicks is at least 10 days
for (id in unique(dat4$ind)){
  afr <- unique(dat4$age_e[dat4$ind==id])
  dat4$agef[dat4$ind==id] <- dat4$year[dat4$ind==id] - min(dat4$year[dat4$ind==id]) + afr   # Actual age of female
}
dat4$age3 <- dat4$agef
dat4$age3[dat4$agef>3] <- 3            # allocate actual age of female to 3 age classes

# 4.1. Sex ratio
csr <- subset(dat4, (!is.na(dat4$fchicks)|!is.na(dat4$mchicks)))    
csr$mchicks[is.na(csr$mchicks)] <- 0
csr$fchicks[is.na(csr$fchicks)] <- 0
csr$total <- csr$mchicks + csr$fchicks

# 4.2. Reproductive output, given success
dch <- subset(dat4, dat4$chicks!=0)  # Remove data with missing information
y.pro <- cbind(dch$chicks, dch$age3, dch$year)
colnames(y.pro) <- c('chicks', 'agef', 'year')



# 5. Population count data (all eyries in every year) 1968-2014
dat5 <- read.csv('Population-count_data.csv', header=TRUE, sep=';')



# 6. Age ratio from observed feathers of females
dat6 <- read.csv('Age-distribution_data.csv', header=TRUE, sep=';')
pr1y <- rbind(dat6[,1]/dat6[,5], dat6[,2]/dat6[,6])
pr2y <- rbind(dat6[,3]/dat6[,5], dat6[,4]/dat6[,6])
pr1y[pr1y==0] <- 0.001
pr2y[pr2y==0] <- 0.001

w1 <- rbind(dat6[,1], dat6[,2])
w2 <- rbind(dat6[,3], dat6[,4])
wT <- rbind(dat6[,1]+dat6[,3]+dat6[,5], dat6[,2]+dat6[,4]+dat6[,6])

w <- array(NA, dim=c(2,nrow(dat6),3))
w[1,,] <- cbind(dat6$X1yF, dat6$X2yF, dat6$X3yplusF)
w[2,,] <- cbind(dat6$X1yM, dat6$X2yM, dat6$X3yplusM)



#########################################
#
# Bundle data
#
#########################################

# Priors for the initial stage-structured population sizes
p1 <- c(dUnif(1,28), rep(0, 50-28))
p2 <- c(dUnif(1,8), rep(0, 50-8))
p3 <- c(dUnif(1,12), rep(0, 50-12))
p4 <- c(dUnif(1,14), rep(0, 50-14))
p5 <- c(dUnif(1,44), rep(0, 50-44))
pNinit <- rbind(p1, p2, p3, p4, p5)

# Bundle data and constants
gh.data <- list(marr=marr, marrD=marrD, lcountf=log(dat5[,2]), lcountm=log(dat5[,2]), w=w, psy=psucc[, 1, ], fchicks=csr[, 'fchicks'], yc=y.pro[, 'chicks'])

gh.constants <- list(rel=cbind(rowSums(marr[,,1]), rowSums(marr[,,2])), relD=cbind(rowSums(marrD[,,1]), rowSums(marrD[,,2])), pst=psucc[, 2, ], tchicks=csr[, 'total'], nyears=dim(marr)[2], agef.chicks=y.pro[, 'agef'], year.chicks=y.pro[,'year']-1967, nchicks=nrow(y.pro), agef.sr=csr[, 'age3'], year.sr=csr[,'year']-1967, nsr=nrow(csr), pNinit=pNinit, wT=apply(w, c(1,2), sum))



#########################################
#
# Wite the IPM file in NIMBLE
#
#########################################

# Model 1 
# s(s*a3*t), p(a2*t), r(t), eta(a3*t), rho(a3*t), xi(a3*t), alpha(s*a2*t), omega(f:-; m:t)
# (This is the main model)


# Write NIMBLE model file
code_ipm1 <- nimbleCode({

# 1. Priors and linear models
# 1.1. CMR data
for (t in 1:(nyears-1)){
  for (k in 1:2){
    phi[k,t] <- s[3,k,t]
    for (a in 1:2){
      p[a,k,t] <- ilogit(lp[a,t])
    } #a

# 1.2. Dead recovery data
    for (a in 1:3){
      ls[a,k,t] ~ dnorm(mean.ls[a,k], sd=sigma.s[a,k])
      s[a,k,t] <- ilogit(ls[a,k,t])
    } #a
    r[k,t] <- ilogit(lr[t])
  } #k
  for (a in 1:2){
    lp[a,t] ~ dnorm(mean.lp[a], sd=sigma.p[a])
  } #a
  lr[t] ~ dnorm(mean.lr, sd=sigma.r)
} #t
for (k in 1:2){
  for (a in 1:3){
    mean.s[a,k] ~ dunif(0, 1)
    mean.ls[a,k] <- logit(mean.s[a,k])
    sigma.s[a,k] ~ dunif(0, 2)
  } #a
} #k
for (a in 1:2){
  mean.p[a] ~ dunif(0, 1)
  mean.lp[a] <- logit(mean.p[a])
  sigma.p[a] ~ dunif(0, 3)
}
mean.r ~ dunif(0, 1)
mean.lr <- logit(mean.r)
sigma.r ~ dunif(0, 3)

# 1.3. Probability of breeding success
for (a in 1:3){
  for (t in 1:nyears){
    leta[a,t] ~ dnorm(mean.leta[a], sd=sigma.eta[a])
  } #t
  eta[a,1:nyears] <- ilogit(leta[a,1:nyears])
  mean.eta[a] ~ dunif(0, 1)
  mean.leta[a] <- logit(mean.eta[a])
  sigma.eta[a] ~ dunif(0, 2)
} #a

# 1.4. Number of chicks given success
for (a in 1:3){
  for (t in 1:nyears){
    rho[a,t] ~ dnorm(mean.rho[a], sd=sigma.rho[a])
  } #t
  mean.rho[a] ~ dnorm(2, 0.01)
  sigma.rho[a] ~ dunif(0, 2)
} #a
sigma.chicks ~ dunif(0, 2)

# 1.5. Sex ratio of chicks
for (a in 1:3){
  for (t in 1:nyears){
    lxi[a,t] ~ dnorm(mean.lxi[a], sd=sigma.xi[a])
  } #t
  xi[a,1:nyears] <- ilogit(lxi[a,1:nyears])
  mean.xi[a] ~ dunif(0, 1)
  mean.lxi[a] <- logit(mean.xi[a])
  sigma.xi[a] ~ dunif(0, 2)
} #a

# 1.6. Recruitment probability (hidden parameter)
for (k in 1:2){
  for (a in 1:2){
    for (t in 1:(nyears-1)){
      lalpha[a,k,t] ~ dnorm(mean.lalpha[a,k], sd=sigma.alpha[a,k])
      alpha[a,k,t] <- ilogit(lalpha[a,k,t])
    } #t
    mean.alpha[a,k] ~ dunif(0, 1)
    mean.lalpha[a,k] <- logit(mean.alpha[a,k])
    sigma.alpha[a,k] ~ dunif(0, 2)
  } #a
} #k

# 1.7. Male immigration
for (t in 1:nyears){
  log.omega[t] ~ dnorm(mean.lomega, sd=sigma.omega)
  omega[t] <- exp(log.omega[t])
}
sigma.omega ~ dunif(0, 2)
mean.omega ~ dunif(0, 20)
mean.lomega <- log(mean.omega)

# 1.8. Residual / observation error
sigma.obs ~ dunif(0.02, 0.3)

# 1.9. Priors for the initial population size: discrete uniform distributions
for (a in 1:5){
  N[a,1,1] ~ dcat(pNinit[a,])
  N[a,2,1] ~ dcat(pNinit[a,])
}
I[1] ~ dpois(omega[1])


# 2. Likelihoods
# 2.1 State-space model for population counts
# Process model of the state-space model: our model of population dynamics
for (t in 1:(nyears-1)){
  # Females
  F[1,t] ~ dpois(N[2,1,t] * eta[1,t] * rho[1,t] * xi[1,t])              # Total number of female fledglings produced by 1y
  F[2,t] ~ dpois(N[4,1,t] * eta[2,t] * rho[2,t] * xi[2,t])              # Total number of female fledglings produced by 2y
  F[3,t] ~ dpois(N[5,1,t] * eta[3,t] * rho[3,t] * xi[3,t])              # Total number of female fledglings produced by adults
  
  N[1,1,t+1] ~ dbin(s[1,1,t] * (1-alpha[1,1,t]), (F[1,t] + F[2,t] + F[3,t]))   # 1y NB
  N[2,1,t+1] ~ dbin(s[1,1,t] * alpha[1,1,t], (F[1,t] + F[2,t] + F[3,t]))       # 1y B 
  N[3,1,t+1] ~ dbin(s[2,1,t] * (1-alpha[2,1,t]), N[1,1,t])
  n[1,1,t] ~ dbin(s[2,1,t] * alpha[2,1,t], N[1,1,t])
  n[2,1,t] ~ dbin(s[2,1,t], N[2,1,t])  
  N[4,1,t+1] <- n[1,1,t] + n[2,1,t]  
  N[5,1,t+1] ~ dbin(s[3,1,t], (N[3,1,t] + N[4,1,t] + N[5,1,t]))
  
  # First-time breeders of different ages
  FB[1,1,t] <- N[2,1,t+1]
  FB[2,1,t] <- n[1,1,t]
  FB[3,1,t] ~ dbin(s[3,1,t], N[3,1,t])
  
  # Males
  M[1,t] ~ dpois(N[2,1,t] * eta[1,t] * rho[1,t] * (1-xi[1,t]))            # Total number of male fledglings produced by 1y
  M[2,t] ~ dpois(N[4,1,t] * eta[2,t] * rho[2,t] * (1-xi[2,t]))            # Total number of male fledglings produced by 2y
  M[3,t] ~ dpois(N[5,1,t] * eta[3,t] * rho[3,t] * (1-xi[3,t]))            # Total number of male fledglings produced by adults

  N[1,2,t+1] ~ dbin(s[1,2,t] * (1-alpha[1,2,t]), (M[1,t] + M[2,t] + M[3,t]))   # 1y NB
  N[2,2,t+1] ~ dbin(s[1,2,t] * alpha[1,2,t], (M[1,t] + M[2,t] + M[3,t]))       # 1y B  
  N[3,2,t+1] ~ dbin(s[2,2,t] * (1-alpha[2,2,t]), N[1,2,t])
  n[1,2,t] ~ dbin(s[2,2,t] * alpha[2,2,t], N[1,2,t])
  n[2,2,t] ~ dbin(s[2,2,t], N[2,2,t])
  N[4,2,t+1]  <- n[1,2,t] + n[2,2,t]
  N[5,2,t+1] ~ dbin(s[3,2,t], (N[3,2,t] + N[4,2,t] + N[5,2,t] + I[t]))
  I[t+1] ~ dpois(omega[t+1])

  # First-time breeders of different ages
  FB[1,2,t] <- N[2,2,t+1]
  FB[2,2,t] <- n[1,2,t]
  FB[3,2,t] ~ dbin(s[3,2,t], N[3,2,t])
}

# Number of fledglings produced in the last study year
F[1,nyears] ~ dpois(N[2,1,nyears] * eta[1,nyears] * rho[1,nyears] * xi[1,nyears])              
F[2,nyears] ~ dpois(N[4,1,nyears] * eta[2,nyears] * rho[2,nyears] * xi[2,nyears])              
F[3,nyears] ~ dpois(N[5,1,nyears] * eta[3,nyears] * rho[3,nyears] * xi[3,nyears])              
M[1,nyears] ~ dpois(N[2,1,nyears] * eta[1,nyears] * rho[1,nyears] * (1-xi[1,nyears]))        
M[2,nyears] ~ dpois(N[4,1,nyears] * eta[2,nyears] * rho[2,nyears] * (1-xi[2,nyears]))
M[3,nyears] ~ dpois(N[5,1,nyears] * eta[3,nyears] * rho[3,nyears] * (1-xi[3,nyears]))   

# Observation model of the state-space model
for (t in 1:nyears){
  lcountf[t] ~ dnorm(log(N[2,1,t] + N[4,1,t] + N[5,1,t]), sd=sigma.obs)
  lcountm[t] ~ dnorm(log(N[2,2,t] + N[4,2,t] + N[5,2,t] + I[t]), sd=sigma.obs)
}


# 2.2. Observed age distribution of breeding individuals
for (t in 1:nyears){
  for (i in 1:2){       # sex
    w[i,t,1:3] ~ dmulti(prAge[i,t,1:3], wT[i,t])
  } #i
  # females
  prAge[1,t,1] <- N[2,1,t] / (N[2,1,t] + N[4,1,t] + N[5,1,t])
  prAge[1,t,2] <- N[4,1,t] / (N[2,1,t] + N[4,1,t] + N[5,1,t])
  prAge[1,t,3] <- N[5,1,t] / (N[2,1,t] + N[4,1,t] + N[5,1,t])	
  # males
  prAge[2,t,1] <- N[2,2,t] / (N[2,2,t] + N[4,2,t] + N[5,2,t] + I[t])
  prAge[2,t,2] <- N[4,2,t] / (N[2,2,t] + N[4,2,t] + N[5,2,t] + I[t])
  prAge[2,t,3] <- (N[5,2,t]  + I[t])/ (N[2,2,t] + N[4,2,t] + N[5,2,t] + I[t])	
} #t


# 2.3. Capture-recapture model (CJS model with multinomial likelihood)
for (k in 1:2){
  for (t in 1: (nyears-1)){
  marr[t,1:nyears,k] ~ dmulti(pr[t,1:nyears,k], rel[t,k])
  } #t
  # Define the cell probabilities of the m-arrays
  for (t in 1:(nyears-1)){
    # Main diagonal
    q[1,k,t] <- 1-p[1,k,t]
    q[2,k,t] <- 1-p[2,k,t]
    pr[t,t,k] <- phi[k,t] * p[1,k,t]
    # Further above main diagonal
    for (j in (t+2):(nyears-1)){
      pr[t,j,k] <- prod(phi[k,(t):j]) * q[1,k,t] * prod(q[2,k,(t+1):(j-1)]) * p[2,k,j] 
    } #j
    # Below main diagonal
    for (j in 1:(t-1)){
      pr[t,j,k] <- 0
    } #j
  } #t
  # One above main diagonal
  for (t in 1:(nyears-2)){
    pr[t,t+1,k] <- phi[k,t] * phi[k,t+1] * q[1,k,t] * p[2,k,t+1]
  } #t
  # Last column: probability of non-recapture
  for (t in 1:(nyears-1)){
    pr[t,nyears,k] <- 1-sum(pr[t,1:(nyears-1),k])
  } #t
} #k


# 2.4. Dead-recovery model
for (k in 1:2){
  for (t in 1:(nyears-1)){
    marrD[t,1:nyears,k] ~ dmulti(prD[t,1:nyears,k], relD[t,k])
  } #t
  # Define the cell probabilities of the m-array
  for (t in 1:(nyears-1)){
    # Main diagonal
    prD[t,t,k] <- (1-s[1,k,t]) * r[k,t]
    # Further than three above main diagonal
    for (j in (t+3):(nyears-1)){
      prD[t,j,k] <- s[1,k,t] * s[2,k,t+1] * prod(s[3,k,(t+2):(j-1)]) * (1-s[3,k,j]) * r[k,j]
    } #j
    # Below main diagonal
    for (j in 1:(t-1)){
      prD[t,j,k] <- 0
    } #j
  } #t
  # One above main diagonal
  for (t in 1:(nyears-2)){
    prD[t,t+1,k] <- s[1,k,t] * (1-s[2,k,t+1]) * r[k,t+1]
  } #t
  # Two above main diagonal
  for (t in 1:(nyears-3)){
    prD[t,t+2,k] <- s[1,k,t] * s[2,k,t+1] * (1-s[3,k,t+2]) * r[k,t+2]
  } #t
  # Last column: probability of non-recovery
  for (t in 1:(nyears-1)){
    prD[t,nyears,k] <- 1-sum(prD[t,1:(nyears-1),k])
  } #t
} #k


# 2.5. Probability of breeding success
# Define the multinomial likelihood
for (a in 1:3){
  for (t in 1:nyears){
    psy[a,t] ~ dbinom(eta[a,t], pst[a,t])
  } #t
} #a


# 2.6. Sex ratio data
for (i in 1:nsr){
  fchicks[i] ~ dbin(xi[agef.sr[i], year.sr[i]], tchicks[i])
}

# 2.7. Number of fledglings, given success
# Define the normal likelihood
for (i in 1:nchicks){
  yc[i] ~ dnorm(rho[agef.chicks[i], year.chicks[i]], sd=sigma.chicks)
}
})



# Parameters to monitor
parameters <- c('N', 'F', 'M', 'I', 'FB', 's', 'alpha', 'phi', 'xi', 'eta', 'rho', 'p', 'r', 'mean.s', 'sigma.s', 'mean.alpha', 'sigma.alpha', 'mean.xi', 'sigma.xi', 'mean.eta', 'sigma.eta', 'mean.rho', 'sigma.rho', 'sigma.chicks', 'mean.omega', 'sigma.omega', 'mean.p', 'sigma.p', 'mean.r', 'sigma.r', 'sigma.obs')

# Build model
model_ipm1 <- nimbleModel(code_ipm1, data=gh.data, constants=gh.constants, inits=inits1(), calculate=FALSE)

# Compile the model
Cmodel_ipm1 <- compileNimble(model_ipm1)

# Build MCMC
conf_ipm1 <- configureMCMC(Cmodel_ipm1, enableWAIC=TRUE, useConjugacy=FALSE, monitors=parameters)
mcmc_ipm1 <- buildMCMC(conf_ipm1, useConjugacy=FALSE)

# Compile MCMC
Cmcmc_ipm1 <- compileNimble(mcmc_ipm1, project=model_ipm1)

# Run the MCMC
res_ipm1 <- runMCMC(Cmcmc_ipm1, niter=110000, nburnin=10000, thin=50, nchains=3, WAIC=TRUE, samplesAsCodaMCMC=TRUE)

# Inspect results
MCMCsummary(res_ipm1$samples, round=3)

# Save results
save(res_ipm1, code_ipm1, gh.data, gh.constants, file='Model1.Rdata')



#############################################

# Model 2
# s(s*a3*t), p(a2*t), r(t), eta(a3*t), rho(a3*t), xi(a3*t), alpha(s*a2*t), omega(f:-; m:t)
# - This model is the same as model 1, but it includes density-dependence in all demographic parameters (with the exception of immigration)


# Write NIMBLE model file
code_ipm2 <- nimbleCode({

# 1. Priors and linear models
# 1.1. CMR data
for (t in 1:(nyears-1)){
  for (k in 1:2){
    phi[k,t] <- s[3,k,t]
    for (a in 1:2){
      p[a,k,t] <- ilogit(lp[a,t])
    } #a

  # 1.2. Dead recovery data
    for (a in 1:3){
      ls[a,k,t] <- mean.ls[a,k] + beta.s[a,k] * (sum(N[1:5,1,t]) - 44)/10 + eps.s[a,k,t]
      eps.s[a,k,t] ~ dnorm(0, sd=sigma.s[a,k])
      s[a,k,t] <- ilogit(ls[a,k,t])
    } #a
    r[k,t] <- ilogit(lr[t])
  } #k
  for (a in 1:2){
    lp[a,t] ~ dnorm(mean.lp[a], sd=sigma.p[a])
  } #a
  lr[t] ~ dnorm(mean.lr, sd=sigma.r)
} #t
for (k in 1:2){
  for (a in 1:3){
    mean.s[a,k] ~ dunif(0, 1)
    mean.ls[a,k] <- logit(mean.s[a,k])
    sigma.s[a,k] ~ dunif(0, 2)
	beta.s[a,k] ~ dnorm(0, sd=1)
  } #a
} #k
for (a in 1:2){
  mean.p[a] ~ dunif(0, 1)
  mean.lp[a] <- logit(mean.p[a])
  sigma.p[a] ~ dunif(0, 3)
}
mean.r ~ dunif(0, 1)
mean.lr <- logit(mean.r)
sigma.r ~ dunif(0, 3)

# 1.3. Probability of breeding success
for (a in 1:3){
  for (t in 1:nyears){
	leta[a,t] <- mean.leta[a] + beta.eta[a] * (sum(N[1:5,1,t]) - 44)/10 + eps.eta[a,t]
	eps.eta[a,t] ~ dnorm(0, sd=sigma.eta[a])    
    eta[a,t] <- ilogit(leta[a,t])
  } #t
  mean.eta[a] ~ dunif(0, 1)
  mean.leta[a] <- logit(mean.eta[a])
  sigma.eta[a] ~ dunif(0, 2)
  beta.eta[a] ~ dnorm(0, sd=1)
} #a

# 1.4. Number of chicks given success
for (a in 1:3){
  for (t in 1:nyears){
	rho[a,t] <- mean.rho[a] + beta.rho[a] * (sum(N[1:5,1,t]) - 44)/10 + eps.rho[a,t]
	eps.rho[a,t] ~ dnorm(0, sd=sigma.rho[a])    
  } #t
  mean.rho[a] ~ dnorm(2, 0.01)
  sigma.rho[a] ~ dunif(0, 2)
  beta.rho[a] ~ dnorm(0, sd=1)
} #a
sigma.chicks ~ dunif(0, 2)

# 1.5. Sex ratio of chicks
for (a in 1:3){
  for (t in 1:nyears){
	lxi[a,t] <- mean.lxi[a] + beta.xi[a] * (sum(N[1:5,1,t]) - 44)/10 + eps.xi[a,t]
	eps.xi[a,t] ~ dnorm(0, sd=sigma.xi[a])    
    xi[a,t] <- ilogit(lxi[a,t])
  } #t
  mean.xi[a] ~ dunif(0, 1)
  mean.lxi[a] <- logit(mean.xi[a])
  sigma.xi[a] ~ dunif(0, 2)
  beta.xi[a] ~ dnorm(0, sd=1)
} #a

# 1.6. Recruitment probability (hidden parameter)
for (k in 1:2){
  for (a in 1:2){
    for (t in 1:(nyears-1)){
  	  lalpha[a,k,t] <- mean.lalpha[a,k] + beta.alpha[a,k] * (sum(N[1:5,1,t]) - 44)/10 + eps.alpha[a,k,t]
      eps.alpha[a,k,t] ~ dnorm(0, sd=sigma.alpha[a,k])    
      alpha[a,k,t] <- ilogit(lalpha[a,k,t])
    } #t
    mean.alpha[a,k] ~ dunif(0, 1)
    mean.lalpha[a,k] <- logit(mean.alpha[a,k])
    sigma.alpha[a,k] ~ dunif(0, 2)
    beta.alpha[a,k] ~ dnorm(0, sd=1)
  } #a
} #k

# 1.7. Male immigration
for (t in 1:nyears){
  log.omega[t] ~ dnorm(mean.lomega, sd=sigma.omega)
  omega[t] <- exp(log.omega[t])
}
sigma.omega ~ dunif(0, 2)
mean.omega ~ dunif(0, 20)
mean.lomega <- log(mean.omega)

# 1.8. Residual / observation error
sigma.obs ~ dunif(0.02, 0.3)

# 1.9. Priors for the initial population size: discrete uniform distributions
for (a in 1:5){
  N[a,1,1] ~ dcat(pNinit[a,])
  N[a,2,1] ~ dcat(pNinit[a,])
}
I[1] ~ dpois(omega[1])


# 2. Likelihoods
# 2.1 State-space model for population counts
# Process model of the state-space model: our model of population dynamics
for (t in 1:(nyears-1)){
  # Females
  F[1,t] ~ dpois(N[2,1,t] * eta[1,t] * rho[1,t] * xi[1,t])              # Total number of female fledglings produced by 1y
  F[2,t] ~ dpois(N[4,1,t] * eta[2,t] * rho[2,t] * xi[2,t])              # Total number of female fledglings produced by 2y
  F[3,t] ~ dpois(N[5,1,t] * eta[3,t] * rho[3,t] * xi[3,t])              # Total number of female fledglings produced by adults
  
  N[1,1,t+1] ~ dbin(s[1,1,t] * (1-alpha[1,1,t]), (F[1,t] + F[2,t] + F[3,t]))   # 1y NB
  N[2,1,t+1] ~ dbin(s[1,1,t] * alpha[1,1,t], (F[1,t] + F[2,t] + F[3,t]))       # 1y B  
  N[3,1,t+1] ~ dbin(s[2,1,t] * (1-alpha[2,1,t]), N[1,1,t])
  n[1,1,t] ~ dbin(s[2,1,t] * alpha[2,1,t], N[1,1,t])
  n[2,1,t] ~ dbin(s[2,1,t], N[2,1,t])  
  N[4,1,t+1] <- n[1,1,t] + n[2,1,t]  
  N[5,1,t+1] ~ dbin(s[3,1,t], (N[3,1,t] + N[4,1,t] + N[5,1,t]))
  
  # First-time breeders of different ages
  FB[1,1,t] <- N[2,1,t+1]
  FB[2,1,t] <- n[1,1,t]
  FB[3,1,t] ~ dbin(s[3,1,t], N[3,1,t])
  
  # Males
  M[1,t] ~ dpois(N[2,1,t] * eta[1,t] * rho[1,t] * (1-xi[1,t]))            # Total number of male fledglings produced by 1y
  M[2,t] ~ dpois(N[4,1,t] * eta[2,t] * rho[2,t] * (1-xi[2,t]))            # Total number of male fledglings produced by 2y
  M[3,t] ~ dpois(N[5,1,t] * eta[3,t] * rho[3,t] * (1-xi[3,t]))            # Total number of male fledglings produced by adults

  N[1,2,t+1] ~ dbin(s[1,2,t] * (1-alpha[1,2,t]), (M[1,t] + M[2,t] + M[3,t]))   # 1y NB
  N[2,2,t+1] ~ dbin(s[1,2,t] * alpha[1,2,t], (M[1,t] + M[2,t] + M[3,t]))       # 1y B 
  N[3,2,t+1] ~ dbin(s[2,2,t] * (1-alpha[2,2,t]), N[1,2,t])
  n[1,2,t] ~ dbin(s[2,2,t] * alpha[2,2,t], N[1,2,t])
  n[2,2,t] ~ dbin(s[2,2,t], N[2,2,t])
  N[4,2,t+1]  <- n[1,2,t] + n[2,2,t]
  N[5,2,t+1] ~ dbin(s[3,2,t], (N[3,2,t] + N[4,2,t] + N[5,2,t] + I[t]))
  I[t+1] ~ dpois(omega[t+1])

  # First-time breeders of different ages
  FB[1,2,t] <- N[2,2,t+1]
  FB[2,2,t] <- n[1,2,t]
  FB[3,2,t] ~ dbin(s[3,2,t], N[3,2,t])
}

# Number of fledglings produced in the last study year
F[1,nyears] ~ dpois(N[2,1,nyears] * eta[1,nyears] * rho[1,nyears] * xi[1,nyears])              
F[2,nyears] ~ dpois(N[4,1,nyears] * eta[2,nyears] * rho[2,nyears] * xi[2,nyears])              
F[3,nyears] ~ dpois(N[5,1,nyears] * eta[3,nyears] * rho[3,nyears] * xi[3,nyears])              
M[1,nyears] ~ dpois(N[2,1,nyears] * eta[1,nyears] * rho[1,nyears] * (1-xi[1,nyears]))        
M[2,nyears] ~ dpois(N[4,1,nyears] * eta[2,nyears] * rho[2,nyears] * (1-xi[2,nyears]))
M[3,nyears] ~ dpois(N[5,1,nyears] * eta[3,nyears] * rho[3,nyears] * (1-xi[3,nyears]))   

# Observation model of the state-space model
for (t in 1:nyears){
  lcountf[t] ~ dnorm(log(N[2,1,t] + N[4,1,t] + N[5,1,t]), sd=sigma.obs)
  lcountm[t] ~ dnorm(log(N[2,2,t] + N[4,2,t] + N[5,2,t] + I[t]), sd=sigma.obs)
}


# 2.2. Observed age distribution of breeding individuals
for (t in 1:nyears){
  for (i in 1:2){       # sex
    w[i,t,1:3] ~ dmulti(prAge[i,t,1:3], wT[i,t])
  }
  # females
  prAge[1,t,1] <- N[2,1,t] / (N[2,1,t] + N[4,1,t] + N[5,1,t])
  prAge[1,t,2] <- N[4,1,t] / (N[2,1,t] + N[4,1,t] + N[5,1,t])
  prAge[1,t,3] <- N[5,1,t] / (N[2,1,t] + N[4,1,t] + N[5,1,t])	
  # males
  prAge[2,t,1] <- N[2,2,t] / (N[2,2,t] + N[4,2,t] + N[5,2,t] + I[t])
  prAge[2,t,2] <- N[4,2,t] / (N[2,2,t] + N[4,2,t] + N[5,2,t] + I[t])
  prAge[2,t,3] <- (N[5,2,t]  + I[t])/ (N[2,2,t] + N[4,2,t] + N[5,2,t] + I[t])	

}

# 2.3. Capture-recapture model (CJS model with multinomial likelihood)
for (k in 1:2){
  for (t in 1: (nyears-1)){
  marr[t,1:nyears,k] ~ dmulti(pr[t,1:nyears,k], rel[t,k])
  }
  # Define the cell probabilities of the m-arrays
  for (t in 1:(nyears-1)){
    # Main diagonal
    q[1,k,t] <- 1-p[1,k,t]
    q[2,k,t] <- 1-p[2,k,t]
    pr[t,t,k] <- phi[k,t] * p[1,k,t]
    # Further above main diagonal
    for (j in (t+2):(nyears-1)){
      pr[t,j,k] <- prod(phi[k,(t):j]) * q[1,k,t] * prod(q[2,k,(t+1):(j-1)]) * p[2,k,j] 
    } #j
    # Below main diagonal
    for (j in 1:(t-1)){
      pr[t,j,k] <- 0
    } #j
  } #t
  # One above main diagonal
  for (t in 1:(nyears-2)){
    pr[t,t+1,k] <- phi[k,t] * phi[k,t+1] * q[1,k,t] * p[2,k,t+1]
  }
  # Last column: probability of non-recapture
  for (t in 1:(nyears-1)){
    pr[t,nyears,k] <- 1-sum(pr[t,1:(nyears-1),k])
  } #t
} #k


# 2.4. Dead-recovery model
for (k in 1:2){
  for (t in 1:(nyears-1)){
    marrD[t,1:nyears,k] ~ dmulti(prD[t,1:nyears,k], relD[t,k])
  }
  # Define the cell probabilities of the m-array
  for (t in 1:(nyears-1)){
    # Main diagonal
    prD[t,t,k] <- (1-s[1,k,t]) * r[k,t]
    # Further than three above main diagonal
    for (j in (t+3):(nyears-1)){
      prD[t,j,k] <- s[1,k,t] * s[2,k,t+1] * prod(s[3,k,(t+2):(j-1)]) * (1-s[3,k,j]) * r[k,j]
    } #j
    # Below main diagonal
    for (j in 1:(t-1)){
      prD[t,j,k] <- 0
    } #j
  } #t
  # One above main diagonal
  for (t in 1:(nyears-2)){
    prD[t,t+1,k] <- s[1,k,t] * (1-s[2,k,t+1]) * r[k,t+1]
  } #t
  # Two above main diagonal
  for (t in 1:(nyears-3)){
    prD[t,t+2,k] <- s[1,k,t] * s[2,k,t+1] * (1-s[3,k,t+2]) * r[k,t+2]
  } #t
  # Last column: probability of non-recovery
  for (t in 1:(nyears-1)){
    prD[t,nyears,k] <- 1-sum(prD[t,1:(nyears-1),k])
  } #t
} #k


# 2.5. Probability of breeding success
# Define the multinomial likelihood
for (a in 1:3){
  for (t in 1:nyears){
    psy[a,t] ~ dbinom(eta[a,t], pst[a,t])
  } #t
} #a


# 2.6. Sex ratio data
for (i in 1:nsr){
  fchicks[i] ~ dbin(xi[agef.sr[i], year.sr[i]], tchicks[i])
}

# 2.7. Productivity data
# Define the normal likelihood
for (i in 1:nchicks){
  yc[i] ~ dnorm(rho[agef.chicks[i], year.chicks[i]], sd=sigma.chicks)
}
})



# Parameters to monitor
parameters <- c('N', 'F', 'M', 'I', 'FB', 's', 'alpha', 'phi', 'xi', 'eta', 'rho', 'p', 'r', 'mean.s', 'sigma.s', 'mean.alpha', 'sigma.alpha', 'mean.xi', 'sigma.xi', 'mean.eta', 'sigma.eta', 'mean.rho', 'sigma.rho', 'sigma.chicks', 'mean.omega', 'sigma.omega', 'mean.p', 'sigma.p', 'mean.r', 'sigma.r', 'sigma.obs', 'beta.s', 'beta.eta', 'beta.rho', 'beta.xi', 'beta.alpha', 'eps.s', 'eps.eta', 'eps.rho', 'eps.xi', 'eps.alpha')

# Build model
model_ipm2 <- nimbleModel(code_ipm2, data=gh.data, constants=gh.constants, inits=inits2(), calculate=FALSE)

# Compile the model
Cmodel_ipm2 <- compileNimble(model_ipm2)

# Build MCMC
conf_ipm2 <- configureMCMC(Cmodel_ipm2, enableWAIC=TRUE, useConjugacy=FALSE, monitors=parameters)
mcmc_ipm2 <- buildMCMC(conf_ipm2, useConjugacy=FALSE)

# Compile MCMC
Cmcmc_ipm2 <- compileNimble(mcmc_ipm2, project=model_ipm2)

# Run the MCMC
res_ipm2 <- runMCMC(Cmcmc_ipm2, niter=110000, nburnin=10000, thin=50, nchains=3, WAIC=TRUE, samplesAsCodaMCMC=TRUE)

# Inspect results
MCMCsummary(res_ipm2$samples, round=3)

# Save results
save(res_ipm2, code_ipm2, gh.data, gh.constants, file='Model2.Rdata')



#############################################

# Model 3
# s(s*a3*t), p(a2*t), r(t), eta(a3*t), rho(a3*t), xi(a3*t), alpha(s*a2*t), omega(f:-; m:t)
# - This model is the same as model 1, but it consideres density-dependence on juvenile female survival


# Write NIMBLE model file
code_ipm3 <- nimbleCode({

# 1. Priors and linear models
# 1.1. CMR data
for (t in 1:(nyears-1)){
  for (k in 1:2){
    phi[k,t] <- s[3,k,t]
    for (a in 1:2){
      p[a,k,t] <- ilogit(lp[a,t])
    } #a
    r[k,t] <- ilogit(lr[t])
  } #k

  # 1.2. Dead recovery data
  ls[1,1,t] <- mean.ls[1,1] + beta.s * (sum(N[,1,t])-50)/10 + eps.s[1,1,t]
  ls[2,1,t] <- mean.ls[2,1] + eps.s[2,1,t]
  ls[3,1,t] <- mean.ls[3,1] + eps.s[3,1,t]
  ls[1,2,t] <- mean.ls[1,2] + eps.s[1,2,t]
  ls[2,2,t] <- mean.ls[2,2] + eps.s[2,2,t]
  ls[3,2,t] <- mean.ls[3,2] + eps.s[3,2,t]
  for (k in 1:2){
	for (a in 1:3){
	  eps.s[a,k,t] ~ dnorm(0, sd=sigma.s[a,k])
	  s[a,k,t] <- plogis(ls[a,k,t])
    } #a
  } #k
  for (a in 1:2){
    lp[a,t] ~ dnorm(mean.lp[a], sd=sigma.p[a])
  } #a
  lr[t] ~ dnorm(mean.lr, sd=sigma.r)
} #t
for (k in 1:2){
  for (a in 1:3){
    mean.s[a,k] ~ dunif(0, 1)
    mean.ls[a,k] <- logit(mean.s[a,k])
    sigma.s[a,k] ~ dunif(0, 2)
  } #a
} #k
beta.s ~ dnorm(0, sd=1)
for (a in 1:2){
  mean.p[a] ~ dunif(0, 1)
  mean.lp[a] <- logit(mean.p[a])
  sigma.p[a] ~ dunif(0, 3)
}
mean.r ~ dunif(0, 1)
mean.lr <- logit(mean.r)
sigma.r ~ dunif(0, 3)

# 1.3. Probability of breeding success
for (a in 1:3){
  for (t in 1:nyears){
    leta[a,t] ~ dnorm(mean.leta[a], sd=sigma.eta[a])
  } #t
  eta[a,1:nyears] <- ilogit(leta[a,1:nyears])
  mean.eta[a] ~ dunif(0, 1)
  mean.leta[a] <- logit(mean.eta[a])
  sigma.eta[a] ~ dunif(0, 2)
} #a

# 1.4. Number of chicks given success
for (a in 1:3){
  for (t in 1:nyears){
    rho[a,t] ~ dnorm(mean.rho[a], sd=sigma.rho[a])
  } #t
  mean.rho[a] ~ dnorm(2, 0.01)
  sigma.rho[a] ~ dunif(0, 2)
} #a
sigma.chicks ~ dunif(0, 2)

# 1.5. Sex ratio of chicks
for (a in 1:3){
  for (t in 1:nyears){
    lxi[a,t] ~ dnorm(mean.lxi[a], sd=sigma.xi[a])
  } #t
  xi[a,1:nyears] <- ilogit(lxi[a,1:nyears])
  mean.xi[a] ~ dunif(0, 1)
  mean.lxi[a] <- logit(mean.xi[a])
  sigma.xi[a] ~ dunif(0, 2)
} #a

# 1.6. Recruitment probability (hidden parameter)
for (k in 1:2){
  for (a in 1:2){
    for (t in 1:(nyears-1)){
      lalpha[a,k,t] ~ dnorm(mean.lalpha[a,k], sd=sigma.alpha[a,k])
      alpha[a,k,t] <- ilogit(lalpha[a,k,t])
    } #t
    mean.alpha[a,k] ~ dunif(0, 1)
    mean.lalpha[a,k] <- logit(mean.alpha[a,k])
    sigma.alpha[a,k] ~ dunif(0, 2)
  } #a
} #k

# 1.7. Male immigration
for (t in 1:nyears){
  log.omega[t] ~ dnorm(mean.lomega, sd=sigma.omega)
  omega[t] <- exp(log.omega[t])
}
sigma.omega ~ dunif(0, 2)
mean.omega ~ dunif(0, 20)
mean.lomega <- log(mean.omega)

# 1.8. Residual / observation error
sigma.obs ~ dunif(0.02, 0.3)

# 1.9. Priors for the initial population size: discrete uniform distributions
for (a in 1:5){
  N[a,1,1] ~ dcat(pNinit[a,])
  N[a,2,1] ~ dcat(pNinit[a,])
}
I[1] ~ dpois(omega[1])


# 2. Likelihoods
# 2.1 State-space model for population counts
# Process model of the state-space model: our model of population dynamics
for (t in 1:(nyears-1)){
  # Females
  F[1,t] ~ dpois(N[2,1,t] * eta[1,t] * rho[1,t] * xi[1,t])              # Total number of female fledglings produced by 1y
  F[2,t] ~ dpois(N[4,1,t] * eta[2,t] * rho[2,t] * xi[2,t])              # Total number of female fledglings produced by 2y
  F[3,t] ~ dpois(N[5,1,t] * eta[3,t] * rho[3,t] * xi[3,t])              # Total number of female fledglings produced by adults
  
  N[1,1,t+1] ~ dbin(s[1,1,t] * (1-alpha[1,1,t]), (F[1,t] + F[2,t] + F[3,t]))   # 1y NB
  N[2,1,t+1] ~ dbin(s[1,1,t] * alpha[1,1,t], (F[1,t] + F[2,t] + F[3,t]))       # 1y B 
  N[3,1,t+1] ~ dbin(s[2,1,t] * (1-alpha[2,1,t]), N[1,1,t])
  n[1,1,t] ~ dbin(s[2,1,t] * alpha[2,1,t], N[1,1,t])
  n[2,1,t] ~ dbin(s[2,1,t], N[2,1,t]) 
  N[4,1,t+1] <- n[1,1,t] + n[2,1,t]  
  N[5,1,t+1] ~ dbin(s[3,1,t], (N[3,1,t] + N[4,1,t] + N[5,1,t]))
  
  # First-time breeders of different ages
  FB[1,1,t] <- N[2,1,t+1]
  FB[2,1,t] <- n[1,1,t]
  FB[3,1,t] ~ dbin(s[3,1,t], N[3,1,t])
  
  # Males
  M[1,t] ~ dpois(N[2,1,t] * eta[1,t] * rho[1,t] * (1-xi[1,t]))            # Total number of male fledglings produced by 1y
  M[2,t] ~ dpois(N[4,1,t] * eta[2,t] * rho[2,t] * (1-xi[2,t]))            # Total number of male fledglings produced by 2y
  M[3,t] ~ dpois(N[5,1,t] * eta[3,t] * rho[3,t] * (1-xi[3,t]))            # Total number of male fledglings produced by adults

  N[1,2,t+1] ~ dbin(s[1,2,t] * (1-alpha[1,2,t]), (M[1,t] + M[2,t] + M[3,t]))   # 1y NB
  N[2,2,t+1] ~ dbin(s[1,2,t] * alpha[1,2,t], (M[1,t] + M[2,t] + M[3,t]))       # 1y B 
  N[3,2,t+1] ~ dbin(s[2,2,t] * (1-alpha[2,2,t]), N[1,2,t])
  n[1,2,t] ~ dbin(s[2,2,t] * alpha[2,2,t], N[1,2,t])
  n[2,2,t] ~ dbin(s[2,2,t], N[2,2,t])
  N[4,2,t+1]  <- n[1,2,t] + n[2,2,t]
  N[5,2,t+1] ~ dbin(s[3,2,t], (N[3,2,t] + N[4,2,t] + N[5,2,t] + I[t]))
  I[t+1] ~ dpois(omega[t+1])

  # First-time breeders of different ages
  FB[1,2,t] <- N[2,2,t+1]
  FB[2,2,t] <- n[1,2,t]
  FB[3,2,t] ~ dbin(s[3,2,t], N[3,2,t])
}

# Number of fledglings produced in the last study year
F[1,nyears] ~ dpois(N[2,1,nyears] * eta[1,nyears] * rho[1,nyears] * xi[1,nyears])              
F[2,nyears] ~ dpois(N[4,1,nyears] * eta[2,nyears] * rho[2,nyears] * xi[2,nyears])              
F[3,nyears] ~ dpois(N[5,1,nyears] * eta[3,nyears] * rho[3,nyears] * xi[3,nyears])              
M[1,nyears] ~ dpois(N[2,1,nyears] * eta[1,nyears] * rho[1,nyears] * (1-xi[1,nyears]))        
M[2,nyears] ~ dpois(N[4,1,nyears] * eta[2,nyears] * rho[2,nyears] * (1-xi[2,nyears]))
M[3,nyears] ~ dpois(N[5,1,nyears] * eta[3,nyears] * rho[3,nyears] * (1-xi[3,nyears]))   

# Observation process of the state-space model
for (t in 1:nyears){
  lcountf[t] ~ dnorm(log(N[2,1,t] + N[4,1,t] + N[5,1,t]), sd=sigma.obs)
  lcountm[t] ~ dnorm(log(N[2,2,t] + N[4,2,t] + N[5,2,t] + I[t]), sd=sigma.obs)
}


# 2.2. Observed age distribution of breeding individuals
for (t in 1:nyears){
  for (i in 1:2){       # sex
    w[i,t,1:3] ~ dmulti(prAge[i,t,1:3], wT[i,t])
  }
  # females
  prAge[1,t,1] <- N[2,1,t] / (N[2,1,t] + N[4,1,t] + N[5,1,t])
  prAge[1,t,2] <- N[4,1,t] / (N[2,1,t] + N[4,1,t] + N[5,1,t])
  prAge[1,t,3] <- N[5,1,t] / (N[2,1,t] + N[4,1,t] + N[5,1,t])	
  # males
  prAge[2,t,1] <- N[2,2,t] / (N[2,2,t] + N[4,2,t] + N[5,2,t] + I[t])
  prAge[2,t,2] <- N[4,2,t] / (N[2,2,t] + N[4,2,t] + N[5,2,t] + I[t])
  prAge[2,t,3] <- (N[5,2,t]  + I[t])/ (N[2,2,t] + N[4,2,t] + N[5,2,t] + I[t])	
}

# 2.3. Capture-recapture model (CJS model with multinomial likelihood)
for (k in 1:2){
  for (t in 1: (nyears-1)){
  marr[t,1:nyears,k] ~ dmulti(pr[t,1:nyears,k], rel[t,k])
  }
  # Define the cell probabilities of the m-arrays
  for (t in 1:(nyears-1)){
    # Main diagonal
    q[1,k,t] <- 1-p[1,k,t]
    q[2,k,t] <- 1-p[2,k,t]
    pr[t,t,k] <- phi[k,t] * p[1,k,t]
    # Further above main diagonal
    for (j in (t+2):(nyears-1)){
      pr[t,j,k] <- prod(phi[k,(t):j]) * q[1,k,t] * prod(q[2,k,(t+1):(j-1)]) * p[2,k,j] 
    } #j
    # Below main diagonal
    for (j in 1:(t-1)){
      pr[t,j,k] <- 0
    } #j
  } #t
  # One above main diagonal
  for (t in 1:(nyears-2)){
    pr[t,t+1,k] <- phi[k,t] * phi[k,t+1] * q[1,k,t] * p[2,k,t+1]
  }
  # Last column: probability of non-recapture
  for (t in 1:(nyears-1)){
    pr[t,nyears,k] <- 1-sum(pr[t,1:(nyears-1),k])
  } #t
} #k


# 2.4. Dead-recovery model
for (k in 1:2){
  for (t in 1:(nyears-1)){
    marrD[t,1:nyears,k] ~ dmulti(prD[t,1:nyears,k], relD[t,k])
  }
  # Define the cell probabilities of the m-array
  for (t in 1:(nyears-1)){
    # Main diagonal
    prD[t,t,k] <- (1-s[1,k,t]) * r[k,t]
    # Further than three above main diagonal
    for (j in (t+3):(nyears-1)){
      prD[t,j,k] <- s[1,k,t] * s[2,k,t+1] * prod(s[3,k,(t+2):(j-1)]) * (1-s[3,k,j]) * r[k,j]
    } #j
    # Below main diagonal
    for (j in 1:(t-1)){
      prD[t,j,k] <- 0
    } #j
  } #t
  # One above main diagonal
  for (t in 1:(nyears-2)){
    prD[t,t+1,k] <- s[1,k,t] * (1-s[2,k,t+1]) * r[k,t+1]
  } #t
  # Two above main diagonal
  for (t in 1:(nyears-3)){
    prD[t,t+2,k] <- s[1,k,t] * s[2,k,t+1] * (1-s[3,k,t+2]) * r[k,t+2]
  } #t
  # Last column: probability of non-recovery
  for (t in 1:(nyears-1)){
    prD[t,nyears,k] <- 1-sum(prD[t,1:(nyears-1),k])
  } #t
} #k


# 2.5. Probability of breeding success
# Define the multinomial likelihood
for (a in 1:3){
  for (t in 1:nyears){
    psy[a,t] ~ dbinom(eta[a,t], pst[a,t])
  } #t
} #a


# 2.6. Sex ratio data
for (i in 1:nsr){
  fchicks[i] ~ dbin(xi[agef.sr[i], year.sr[i]], tchicks[i])
}


# 2.7. Productivity data
# Define the normal likelihood
for (i in 1:nchicks){
  yc[i] ~ dnorm(rho[agef.chicks[i], year.chicks[i]], sd=sigma.chicks)
}
})



# Parameters to monitor
parameters <- c('N', 'F', 'M', 'I', 'FB', 's', 'alpha', 'phi', 'xi', 'eta', 'rho', 'p', 'r', 'mean.s', 'sigma.s', 'mean.alpha', 'sigma.alpha', 'mean.xi', 'sigma.xi', 'mean.eta', 'sigma.eta', 'mean.rho', 'sigma.rho', 'sigma.chicks', 'mean.omega', 'sigma.omega', 'mean.p', 'sigma.p', 'mean.r', 'sigma.r', 'sigma.obs', 'beta.s', 'eps.s')

# Build model
model_ipm3 <- nimbleModel(code_ipm3, data=gh.data, constants=gh.constants, inits=inits3(), calculate=FALSE)

# Compile the model, together with the functions
Cmodel_ipm3 <- compileNimble(model_ipm3)

# Build MCMC
conf_ipm3 <- configureMCMC(Cmodel_ipm3, enableWAIC=TRUE, useConjugacy=FALSE, monitors=parameters)
mcmc_ipm3 <- buildMCMC(conf_ipm3, useConjugacy=FALSE)

# Compile MCMC
Cmcmc_ipm3 <- compileNimble(mcmc_ipm3, project=model_ipm3)

# Run the MCMC
res_ipm3 <- runMCMC(Cmcmc_ipm3, niter=110000, nburnin=10000, thin=50, nchains=3, WAIC=TRUE, samplesAsCodaMCMC=TRUE)

# Inspect results
MCMCsummary(res_ipm3$samples, round=3)

# Save results
save(res_ipm3, code_ipm3, gh.data, gh.constants, file='Model3.Rdata')



#############################################

# Model 4
# s(s*a3*t), p(a2*t), r(t), eta(a3*t), rho(a3*t), xi(a3*t), alpha(s*a2*t), omega(f:t; m:t)
# - This model is the same as model 1, but it allows female immigration


# Write NIMBLE model file
code_ipm4 <- nimbleCode({

# 1. Priors and linear models
# 1.1. CMR data
for (t in 1:(nyears-1)){
  for (k in 1:2){
    phi[k,t] <- s[3,k,t]
    for (a in 1:2){
      p[a,k,t] <- ilogit(lp[a,t])
    } #a

  # 1.2. Dead recovery data
    for (a in 1:3){
      ls[a,k,t] ~ dnorm(mean.ls[a,k], sd=sigma.s[a,k])
      s[a,k,t] <- ilogit(ls[a,k,t])
    } #a
    r[k,t] <- ilogit(lr[t])
  } #k
  for (a in 1:2){
    lp[a,t] ~ dnorm(mean.lp[a], sd=sigma.p[a])
  } #a
  lr[t] ~ dnorm(mean.lr, sd=sigma.r)
} #t
for (k in 1:2){
  for (a in 1:3){
    mean.s[a,k] ~ dunif(0, 1)
    mean.ls[a,k] <- logit(mean.s[a,k])
    sigma.s[a,k] ~ dunif(0, 2)
  } #a
} #k
for (a in 1:2){
  mean.p[a] ~ dunif(0, 1)
  mean.lp[a] <- logit(mean.p[a])
  sigma.p[a] ~ dunif(0, 3)
}
mean.r ~ dunif(0, 1)
mean.lr <- logit(mean.r)
sigma.r ~ dunif(0, 3)

# 1.3. Probability of breeding success
for (a in 1:3){
  for (t in 1:nyears){
    leta[a,t] ~ dnorm(mean.leta[a], sd=sigma.eta[a])
  } #t
  eta[a,1:nyears] <- ilogit(leta[a,1:nyears])
  mean.eta[a] ~ dunif(0, 1)
  mean.leta[a] <- logit(mean.eta[a])
  sigma.eta[a] ~ dunif(0, 2)
} #a

# 1.4. Number of chicks given success
for (a in 1:3){
  for (t in 1:nyears){
    rho[a,t] ~ dnorm(mean.rho[a], sd=sigma.rho[a])
  } #t
  mean.rho[a] ~ dnorm(2, 0.01)
  sigma.rho[a] ~ dunif(0, 2)
} #a
sigma.chicks ~ dunif(0, 2)

# 1.5. Sex ratio of chicks
for (a in 1:3){
  for (t in 1:nyears){
    lxi[a,t] ~ dnorm(mean.lxi[a], sd=sigma.xi[a])
  } #t
  xi[a,1:nyears] <- ilogit(lxi[a,1:nyears])
  mean.xi[a] ~ dunif(0, 1)
  mean.lxi[a] <- logit(mean.xi[a])
  sigma.xi[a] ~ dunif(0, 2)
} #a

# 1.6. Recruitment probability (hidden parameter)
for (k in 1:2){
  for (a in 1:2){
    for (t in 1:(nyears-1)){
      lalpha[a,k,t] ~ dnorm(mean.lalpha[a,k], sd=sigma.alpha[a,k])
      alpha[a,k,t] <- ilogit(lalpha[a,k,t])
    } #t
    mean.alpha[a,k] ~ dunif(0, 1)
    mean.lalpha[a,k] <- logit(mean.alpha[a,k])
    sigma.alpha[a,k] ~ dunif(0, 2)
  } #a
} #k

# 1.7. Immigration
for (k in 1:2){
  for (t in 1:nyears){
    log.omega[k,t] ~ dnorm(mean.lomega[k], sd=sigma.omega[k])
    omega[k,t] <- exp(log.omega[k,t])
  }
  sigma.omega[k] ~ dunif(0, 2)
  mean.omega[k] ~ dunif(0, 20)
  mean.lomega[k] <- log(mean.omega[k])
}

# 1.8. Residual / observation error
sigma.obs ~ dunif(0.02, 0.3)

# 1.9. Priors for the initial population size: discrete uniform distributions
for (a in 1:5){
  N[a,1,1] ~ dcat(pNinit[a,])
  N[a,2,1] ~ dcat(pNinit[a,])
}
I[1,1] ~ dpois(omega[1,1])
I[2,1] ~ dpois(omega[2,1])


# 2. Likelihoods
# 2.1 State-space model for population counts
# Process model of the state-space model: our model of population dynamics
for (t in 1:(nyears-1)){
  # Females
  F[1,t] ~ dpois(N[2,1,t] * eta[1,t] * rho[1,t] * xi[1,t])              # Total number of female fledglings produced by 1y
  F[2,t] ~ dpois(N[4,1,t] * eta[2,t] * rho[2,t] * xi[2,t])              # Total number of female fledglings produced by 2y
  F[3,t] ~ dpois((N[5,1,t] + I[1,t]) * eta[3,t] * rho[3,t] * xi[3,t])   # Total number of female fledglings produced by adults
  
  N[1,1,t+1] ~ dbin(s[1,1,t] * (1-alpha[1,1,t]), (F[1,t] + F[2,t] + F[3,t]))   # 1y NB
  N[2,1,t+1] ~ dbin(s[1,1,t] * alpha[1,1,t], (F[1,t] + F[2,t] + F[3,t]))       # 1y B 
  N[3,1,t+1] ~ dbin(s[2,1,t] * (1-alpha[2,1,t]), N[1,1,t])
  n[1,1,t] ~ dbin(s[2,1,t] * alpha[2,1,t], N[1,1,t])
  n[2,1,t] ~ dbin(s[2,1,t], N[2,1,t])  
  N[4,1,t+1] <- n[1,1,t] + n[2,1,t]  
  N[5,1,t+1] ~ dbin(s[3,1,t], (N[3,1,t] + N[4,1,t] + N[5,1,t] + I[1,t]))
  I[1,t+1] ~ dpois(omega[1,t+1])

  # First-time breeders of different ages
  FB[1,1,t] <- N[2,1,t+1]
  FB[2,1,t] <- n[1,1,t]
  FB[3,1,t] ~ dbin(s[3,1,t], N[3,1,t])
  
  # Males
  M[1,t] ~ dpois(N[2,1,t] * eta[1,t] * rho[1,t] * (1-xi[1,t]))              # Total number of male fledglings produced by 1y
  M[2,t] ~ dpois(N[4,1,t] * eta[2,t] * rho[2,t] * (1-xi[2,t]))              # Total number of male fledglings produced by 2y
  M[3,t] ~ dpois((N[5,1,t] + I[1,t]) * eta[3,t] * rho[3,t] * (1-xi[3,t]))   # Total number of male fledglings produced by adults

  N[1,2,t+1] ~ dbin(s[1,2,t] * (1-alpha[1,2,t]), (M[1,t] + M[2,t] + M[3,t]))   # 1y NB
  N[2,2,t+1] ~ dbin(s[1,2,t] * alpha[1,2,t], (M[1,t] + M[2,t] + M[3,t]))       # 1y B  
  N[3,2,t+1] ~ dbin(s[2,2,t] * (1-alpha[2,2,t]), N[1,2,t])
  n[1,2,t] ~ dbin(s[2,2,t] * alpha[2,2,t], N[1,2,t])
  n[2,2,t] ~ dbin(s[2,2,t], N[2,2,t])
  N[4,2,t+1]  <- n[1,2,t] + n[2,2,t]
  N[5,2,t+1] ~ dbin(s[3,2,t], (N[3,2,t] + N[4,2,t] + N[5,2,t] + I[2,t]))
  I[2,t+1] ~ dpois(omega[2,t+1])

  # First-time breeders of different ages
  FB[1,2,t] <- N[2,2,t+1]
  FB[2,2,t] <- n[1,2,t]
  FB[3,2,t] ~ dbin(s[3,2,t], N[3,2,t])
}

# Number of fledglings produced in the last study year
F[1,nyears] ~ dpois(N[2,1,nyears] * eta[1,nyears] * rho[1,nyears] * xi[1,nyears])              
F[2,nyears] ~ dpois(N[4,1,nyears] * eta[2,nyears] * rho[2,nyears] * xi[2,nyears])              
F[3,nyears] ~ dpois((N[5,1,nyears] + I[1,nyears]) * eta[3,nyears] * rho[3,nyears] * xi[3,nyears])              
M[1,nyears] ~ dpois(N[2,1,nyears] * eta[1,nyears] * rho[1,nyears] * (1-xi[1,nyears]))        
M[2,nyears] ~ dpois(N[4,1,nyears] * eta[2,nyears] * rho[2,nyears] * (1-xi[2,nyears]))
M[3,nyears] ~ dpois((N[5,1,nyears] + I[1,nyears]) * eta[3,nyears] * rho[3,nyears] * (1-xi[3,nyears]))   

# Observation model of the state-space model
for (t in 1:nyears){
  lcountf[t] ~ dnorm(log(N[2,1,t] + N[4,1,t] + N[5,1,t] + I[1,t]), sd=sigma.obs)
  lcountm[t] ~ dnorm(log(N[2,2,t] + N[4,2,t] + N[5,2,t] + I[2,t]), sd=sigma.obs)
}


# 2.2. Observed age distribution of breeding individuals
for (t in 1:nyears){
  for (i in 1:2){       # sex
    w[i,t,1:3] ~ dmulti(prAge[i,t,1:3], wT[i,t])
  }
  # females
  prAge[1,t,1] <- N[2,1,t] / (N[2,1,t] + N[4,1,t] + N[5,1,t] + I[1,t])
  prAge[1,t,2] <- N[4,1,t] / (N[2,1,t] + N[4,1,t] + N[5,1,t] + I[1,t])
  prAge[1,t,3] <- (N[5,1,t] + I[1,t])/ (N[2,1,t] + N[4,1,t] + N[5,1,t] + I[1,t])	
  # males
  prAge[2,t,1] <- N[2,2,t] / (N[2,2,t] + N[4,2,t] + N[5,2,t] + I[2,t])
  prAge[2,t,2] <- N[4,2,t] / (N[2,2,t] + N[4,2,t] + N[5,2,t] + I[2,t])
  prAge[2,t,3] <- (N[5,2,t] + I[2,t])/ (N[2,2,t] + N[4,2,t] + N[5,2,t] + I[2,t])	
}

# 2.3. Capture-recapture model (CJS model with multinomial likelihood)
for (k in 1:2){
  for (t in 1: (nyears-1)){
  marr[t,1:nyears,k] ~ dmulti(pr[t,1:nyears,k], rel[t,k])
  }
  # Define the cell probabilities of the m-arrays
  for (t in 1:(nyears-1)){
    # Main diagonal
    q[1,k,t] <- 1-p[1,k,t]
    q[2,k,t] <- 1-p[2,k,t]
    pr[t,t,k] <- phi[k,t] * p[1,k,t]
    # Further above main diagonal
    for (j in (t+2):(nyears-1)){
      pr[t,j,k] <- prod(phi[k,(t):j]) * q[1,k,t] * prod(q[2,k,(t+1):(j-1)]) * p[2,k,j] 
    } #j
    # Below main diagonal
    for (j in 1:(t-1)){
      pr[t,j,k] <- 0
    } #j
  } #t
  # One above main diagonal
  for (t in 1:(nyears-2)){
    pr[t,t+1,k] <- phi[k,t] * phi[k,t+1] * q[1,k,t] * p[2,k,t+1]
  } #t
  # Last column: probability of non-recapture
  for (t in 1:(nyears-1)){
    pr[t,nyears,k] <- 1-sum(pr[t,1:(nyears-1),k])
  } #t
} #k


# 2.4. Dead-recovery model
for (k in 1:2){
  for (t in 1:(nyears-1)){
    marrD[t,1:nyears,k] ~ dmulti(prD[t,1:nyears,k], relD[t,k])
  }
  # Define the cell probabilities of the m-array
  for (t in 1:(nyears-1)){
    # Main diagonal
    prD[t,t,k] <- (1-s[1,k,t]) * r[k,t]
    # Further than three above main diagonal
    for (j in (t+3):(nyears-1)){
      prD[t,j,k] <- s[1,k,t] * s[2,k,t+1] * prod(s[3,k,(t+2):(j-1)]) * (1-s[3,k,j]) * r[k,j]
    } #j
    # Below main diagonal
    for (j in 1:(t-1)){
      prD[t,j,k] <- 0
    } #j
  } #t
  # One above main diagonal
  for (t in 1:(nyears-2)){
    prD[t,t+1,k] <- s[1,k,t] * (1-s[2,k,t+1]) * r[k,t+1]
  } #t
  # Two above main diagonal
  for (t in 1:(nyears-3)){
    prD[t,t+2,k] <- s[1,k,t] * s[2,k,t+1] * (1-s[3,k,t+2]) * r[k,t+2]
  } #t
  # Last column: probability of non-recovery
  for (t in 1:(nyears-1)){
    prD[t,nyears,k] <- 1-sum(prD[t,1:(nyears-1),k])
  } #t
} #k


# 2.5. Probability of breeding success
# Define the multinomial likelihood
for (a in 1:3){
  for (t in 1:nyears){
    psy[a,t] ~ dbinom(eta[a,t], pst[a,t])
  } #t
} #a


# 2.6. Sex ratio data
for (i in 1:nsr){
  fchicks[i] ~ dbin(xi[agef.sr[i], year.sr[i]], tchicks[i])
}


# 2.7. Productivity data
# Define the normal likelihood
for (i in 1:nchicks){
  yc[i] ~ dnorm(rho[agef.chicks[i], year.chicks[i]], sd=sigma.chicks)
}
})



# Parameters to monitor
parameters <- c('N', 'F', 'M', 'I', 'FB', 's', 'alpha', 'phi', 'xi', 'eta', 'rho', 'p', 'r', 'mean.s', 'sigma.s', 'mean.alpha', 'sigma.alpha', 'mean.xi', 'sigma.xi', 'mean.eta', 'sigma.eta', 'mean.rho', 'sigma.rho', 'sigma.chicks', 'mean.omega', 'sigma.omega', 'mean.p', 'sigma.p', 'mean.r', 'sigma.r', 'sigma.obs')

# Build model
model_ipm4 <- nimbleModel(code_ipm4, data=gh.data, constants=gh.constants, inits=inits4(), calculate=FALSE)

# Compile the model
Cmodel_ipm4 <- compileNimble(model_ipm4)

# Build MCMC
conf_ipm4 <- configureMCMC(Cmodel_ipm4, enableWAIC=TRUE, useConjugacy=FALSE, monitors=parameters)
mcmc_ipm4 <- buildMCMC(conf_ipm4, useConjugacy=FALSE)

# Compile MCMC
Cmcmc_ipm4 <- compileNimble(mcmc_ipm4, project=model_ipm4)

# Run the MCMC
res_ipm4 <- runMCMC(Cmcmc_ipm4, niter=110000, nburnin=10000, thin=50, nchains=3, WAIC=TRUE, samplesAsCodaMCMC=TRUE)

# Inspect results
MCMCsummary(res_ipm4$samples, round=3)

# Save results
save(res_ipm4, code_ipm4, gh.data, gh.constants, file='Model4.Rdata')