#######################################################################################
#
# Productivity drives the dynamics of a red kite source population that depends on immigration 
#
# Thomas Pfeiffer & Michael Schaub
#
# To be published in Journal of Avian Biology
#
# Article DOI: 10.1111/jav.02984
#
#######################################################################################


# Code for all the analyses performed and for all figures produced in the paper
#
# The complete code is written in R, but some analyses are performed in JAGS
#
# MS, October 22, 2022



########################################################
#
# 1. Load libraries
#
########################################################

library(jagsUI)
library(maps)
library(RColorBrewer)
library(IPMbook)
library(scales)
library(sp)

########################################################
#
# 2. Load custom functions
#
########################################################

setwd('...')
source('Custom_functions.txt')


########################################################
#
#  3. Load data and do some manipulations
#
########################################################

setwd('...')

# 1. Capture-recovery-resighting data from Weimar (1985-2018)
# - Birds marked with "plaste" wing tag are excluded
# - Resightings only during breeding and inside the study area

# 1.1. CMR of indv. marked as juveniles with a ring only (will be joint with the ring-recovery data from Hiddensee)
# - Capture-histories of 1365 individuals (rows) across 34 years (columns, 1985-2018)
# - Codes of the capture-histories
#   - 0: not captured
#   - 1: marked as juvenile
#   - 2: recovered dead in the first year of life
#   - 3: recovered dead after the first year of life

JuvRing <- read.table("CMRJuvRing.csv", header = TRUE, sep = ";")
CHJuvRing <- as.matrix(JuvRing[,2:35])

# Compute vector with occasion of first capture
f1 <- apply(CHJuvRing, 1, getFirst)

# Recode CH
CMR1 <- CHJuvRing 
CMR1[CMR1==2] <- 1
CMR1[CMR1==3] <- 1

  
# 1.2. CMR of indv. marked as juveniles with wing tag ("Draht")
# - Capture-histories of 568 individuals (rows) across 12 years (columns, 2007-2018)
# - Codes of the capture-histories
#   - 0: not captured
#   - 1: marked as juvenile
#   - 2: resighted as 1y non-breeder in the study area
#   - 3: resighted as 2y non-breeder in the study area
#   - 4: resighted as 3y non-breeder in the study area
#   - 5: resighted as adult (older than 3 years) non-breeder in the study area
#   - 6: resighted as  adult breeder in the study area
#   - 7: recovered dead in the first year of life
#   - 8: recovered dead after the first year of life

JuvDraht <- read.table("CMRJuvDraht.csv", header = TRUE, sep = ";")
CHJuvDraht <- as.matrix(JuvDraht[,2:ncol(JuvDraht)])

# Create multistate m-array
marr1 <- marray(CHJuvDraht, 4)          

# Years
year1 <- 2007:2018

# Create vector with the resighting periods
# Period 1: 2008-2014
# Period 2: 2015-2018
period <- c(rep(1, 8), rep(2, 4))


# 1.3. CMR of indv. marked as adults with wing tag ("Draht")
# - Capture-histories of 18 individuals (rows) across 12 years (columns, 2007-2018)
# - Codes of the capture-histories
#   - 0: not captured
#   - 1: marked and resighted as adult breeder (in the study area)
#   - 2: resighted as adult non-breeder in the study area
#   - 3: recovered dead
# - The variable 'fr' is a counter variable. A value of '-1' indicates that the corresponding individual was censored at last observation. This was needed because one individual was actually recaptured, the wing tag removed and a satellite tag was fitted instead. Hence for the wing tag data set it needs to be censored.

AdDraht <- read.table("CMRAdDraht.csv", header = TRUE, sep = ";")
CHAdDraht <- as.matrix(AdDraht[,2:(ncol(AdDraht)-1)])
fr <- AdDraht[,'fr']

# Create m-array
marr2 <- marray(CHAdDraht, 3, freq=fr)

# Years
year2 <- 2007:2018


# 1.4. CMR of indv. marked as adults with satellite transmitter 
# - Capture-histories of 62 individuals (rows) across 17 years (columns, 2002-2018)
# - Codes of the capture-histories
#   - 0: not captured
#   - 1: marked and detected as adult breeder (in the study area)
#   - 2: detected as adult non-breeder in the study area
#   - 3: found dead
# - The variable 'fr' is a counter variable. A value of '-1' indicates that the corresponding individual was censored at last observation. 

AdSat <- read.table("CMRAdSat.csv", header = TRUE, sep = ";")
CHAdSat <- as.matrix(AdSat[,2:(ncol(AdSat)-1)])
fr <- AdSat[,'fr']

# Create m-array
marr3 <- marray(CHAdSat, 1, freq=fr)

# Years
year3 <- 2002:2018


# 1.5. CMR of indv. marked as juveniles with satellite transmitter
# - A total of 13 juveniles were marked in 2017 with a transmitter
# - In year 2018:
#   - 1 marked individual was detected as non-breeder in the study area
#   - 6 marked individuals were alive, but outside the study area 
#   - 6 marked individuals were found dead
#   - 0 marked individual was not encountered
# We store this information in a scalar and a vector:

totalJuvSat <- 13
totalRecov <- c(1, 6, 6, 0)
   
# Years
year4 <- 2017:2018


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

# 2. Ring-recovery data from the Beringungszentrale Hiddensee (1985-2018) (data curtosy: Christof Herrmann)
# - Capture-histories of 25615 individuals (rows) across 34 years (columns, 1985-2018)
# - Codes of the capture-histories
#   - 0: not captured
#   - 1: marked 
#   - 2: recovered dead
# - Additional variables:
#   - 'Age': age at marking (Juv (=nestling), 1y, Ad)
#   - 'WingTag': whether or not the individual had a wing tag (ja, nein)
#   - 'Weimar': whether or not the individual belonged to the Weimar study population (ja, nein)
#   - 'X': longitude of ringing location
#   - 'Y': latitude of ringing location

ReReData <- read.table('RingRecoveries.csv', header = TRUE, sep = ';')
ch.all <- as.matrix(ReReData[,2:35])
age.all <- ReReData[,36]
tag.all <- ReReData[,37]
Weimar <- ReReData[,38]

# Exclude the Weimar birds
ch <- ch.all[Weimar == 'nein',]
tag <- tag.all[Weimar == 'nein']
age <- age.all[Weimar == 'nein']
first <- apply(ch, 1, getFirst)

# Prepare data for juveniles
z.j <- which(age == 'Juv')
ch.j <- ch[z.j,]
# record dead recoveries
ch.j[ch.j==2] <- 1
tag.j <- tag[z.j]

# Ringed first year birds 
z.1 <- which(age == '1y')
ch.1 <- ch[z.1,]
# record dead recoveries
ch.1[ch.1==2] <- 1
tag.1 <- tag[z.1]

# Ringed adult year birds 
z.a <- which(age == 'Ad')
ch.a <- ch[z.a,]
# record dead recoveries
ch.a[ch.a==2] <- 1
tag.a <- tag[z.a]


# Compute sample sizes
# Number of juveniles without wing tags
dim(ch.j[tag.j=='nein', ])[1]                      # Ringed
length(which(rowSums(ch.j[tag.j=='nein', ])==2))   # Recovered dead

# Number of juveniles with wing tags
dim(ch.j[tag.j=='ja', ])[1]                        # Ringed
length(which(rowSums(ch.j[tag.j=='ja', ])==2))     # Recovered dead

# Number of 1y without wing tags
dim(ch.1[tag.1=='nein', ])[1]                      # Ringed
length(which(rowSums(ch.1[tag.1=='nein', ])==2))   # Recovered dead

# Number of 1y with wing tags
dim(ch.1[tag.1=='ja', ])[1]                        # Ringed
length(which(rowSums(ch.1[tag.1=='ja', ])==2))     # Recovered dead

# Number of adults without wing tags
dim(ch.a[tag.a=='nein', ])[1]                      # Ringed
length(which(rowSums(ch.a[tag.a=='nein', ])==2))   # Recovered dead

# Number of adults with wing tags
dim(ch.a[tag.a=='ja', ])[1]                        # Ringed
length(which(rowSums(ch.a[tag.a=='ja', ])==2))     # Recovered dead


# Join the ring-recovery data of individuals marked as juvenile and with a ring only of Hiddensee and Weimar: they can be analysed with the same submodel. 
join.CMR <- rbind(ch.j[tag.j=='nein', ], CMR1)

# Create the m-arrays
marr.jr <- marrayDead(join.CMR)               # only ringed
marr.jw <- marrayDead(ch.j[tag.j=='ja', ])    # with wing tags
marr.1r <- marrayDead(ch.1[tag.1=='nein', ])
marr.1w <- marrayDead(ch.1[tag.1=='ja', ])
marr.ar <- marrayDead(ch.a[tag.a=='nein', ])
marr.aw <- marrayDead(ch.a[tag.a=='ja', ])


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

# 3. Count data from Weimar (1985-2018)
# Data from 1482 broods in the study area (1985-2018)
# The table has the following columns:
# - Jahr: year
# - Brutstatus: indicates whether the brood was successful (at least 1 fledgling) [eBP], or whether the brood failed [neBP]
# - Juv: brood size at the time of fledging
# - Mindest: whether the brood size as accurately recorded (0), or whether the number refers to a minimum (1)

bruten <- read.table("Bruten.csv", header = TRUE, sep = ";")

# Get the annual number of breeding pairs
bp <- table(bruten$Brutstatus, bruten$Jahr)
counts <- colSums(bp)


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

# 4. Data on productivity from Weimar (1985-2018)
# - The are from the same data table as annual breeding pairs

u <- which(is.na(bruten$Juv) & bruten$Brutstatus == "eBP")

incl <- which(!is.na(bruten$Juv))
Juv <- bruten$Juv[incl]
d <- bruten$Mindest[incl]
Jahr <- bruten$Jahr[incl]

# all these 'manipulations and calculations' are because we model also the 'minimum' data
Juv <- c(Juv, rep(1,length(u)))
d <- c(d, rep(1, length(u)))
Jahr <- c(Jahr, bruten$Jahr[u])

nyears <- length(unique(Jahr))
censor <- Juv+1
censor[d==1] <- Juv[d==1]

init.Juv <- Juv+1
init.Juv[d==0] <- NA
Juv[d==1] <- NA

year.total <- 1985:2018



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

# 5. Data on the body condition of the nestlings
# Data on body condition from 1756 nestlings originating from 806 broods conducted in the study area (1989-2018)
# See main text of the paper for the definition of 'condition'
# Each row corresponds to a nestling
# The table has the following columns:
# - Ezu: Body condition of each nestling
# - broodID: identifier of the brood
# - Year:
# - Juv: number of nestlings in the brood

Cond <- read.table("Condition.csv", header = TRUE, sep = ";")




########################################################
#
#  4. Define the IPM in BUGS language
#
########################################################

setwd('...')


# Specify model in BUGS language
cat(file = "M1.txt", "
model {

# Linear models for parameters
# Survival
for (t in 1:(n.occasions1-1)){
  for (j in 1:3){
    s1[j,t] <- s[j,(t+start1-1)]
  } #j
} #t
for (t in 1:(n.occasions2-1)){
  s2[3,t] <- s[3,(t+start2-1)]
}
for (t in 1:(n.occasions3-1)){
  s3[3,t] <- s[3,(t+start3-1)]
}
s4 <- s[1,33]     # Survival from 2017 - 2018   
for (t in 1:n.occasions5){
  for (j in 1:3){
    s5[j,t] <- s[j,(t+start5-1)]
  } #j
} #t
   
# Apply linear models for the different survival rates
for (t in 1:(nyears-1)){
  for (age in 1:3){
    s[age,t] <- ilogit(logit.s[t,age])
  } #age
} #t 
for (t in 1:(nyears-1)){
  logit.s[t,1:3] ~ dmnorm(mus[], taus[,])
}
for (i in 1:3){
  mean.s[i] ~ dunif(0, 1)
  mus[i] <- logit(mean.s[i])
}
taus[1:3,1:3] ~ dwish(Q[,], 4)
sigmas[1:3,1:3] <- inverse(taus[1:3,1:3])
rhos[1] <- sigmas[1,2] / sqrt(sigmas[1,1] * sigmas[2,2])   # Correlation s1 vs s2
rhos[2] <- sigmas[1,3] / sqrt(sigmas[1,1] * sigmas[3,3])   # Correlation s1 vs s3
rhos[3] <- sigmas[2,3] / sqrt(sigmas[2,2] * sigmas[3,3])   # Correlation s1 vs s3


# Age at first reproduction
for (age in 1:2){
  for (t in 1:(n.occasions1-1)){
    a1[age,t] <- a[age,t+start1-1]
  } #t
  # Apply linear model on age at first reproduction
  for (t in 1:(nyears-1)){
    a[age,t] <- ilogit(eta[t,age])
  } #t
} #age


# Breeding probability
for (t in 1:(n.occasions1-1)){
  b1[t] <- b[t+start1]
}
for (t in 1:(n.occasions2-1)){
  b2[t] <- b[t+start2]
}
for (t in 1:(n.occasions3-1)){
  b3[t] <- b[t+start3]
}

# Apply linear model for breeding probability
for (t in 2:nyears){
  b[t] <- ilogit(eta[t,3])
}
b[1] <- mean.b        # Model for the first year


# Site fidelity
for (t in 1:(n.occasions1-1)){
  fi1[t] <- fi[t+start1-1]
}
fi4 <- fi[33]   
   
# Apply linear model for site fidelity
for (t in 1:(nyears-1)){
  fi[t] <- mean.fi
}
mean.fi ~ dunif(0, 1)


# Immigration
for (t in 1:nyears){
  im[t] <- exp(log.im[t])
  log.im[t] ~ dnorm(mean.log.im, tau.im)
}
mean.im ~ dunif(0.0001, 1)
mean.log.im <- log(mean.im)
sigma.im ~ dunif(0.0001, 3)
tau.im <- pow(sigma.im, -2)


# Productivity
for (t in 2:nyears){
  f[t] <- exp(eta[t-1,4])      # Construction to allow that the correlation between the rates refers to the same years
}
log.f[1] ~ dunif(-10, 10)    # Prior for productivity in the first year. Is a fixed time effect, but this is justified, because there are many data
f[1] <- exp(log.f[1])

# Multivariate Normal for a, b and f with the inverse Wishart prior
for (t in 1:nyears){
  eta[t,1:4] ~ dmnorm(mu[], tau[,])
}
tau[1:4,1:4] ~ dwish(W[,], 5)
sigma[1:4,1:4] <- inverse(tau[1:4,1:4])
rho[1] <- sigma[1,2] / sqrt(sigma[1,1] * sigma[2,2])   # Correlation a1 vs a2
rho[2] <- sigma[1,3] / sqrt(sigma[1,1] * sigma[3,3])   # Correlation a1 vs b
rho[3] <- sigma[1,4] / sqrt(sigma[1,1] * sigma[4,4])   # Correlation a1 vs f
rho[4] <- sigma[2,3] / sqrt(sigma[2,2] * sigma[3,3])   # Correlation a2 vs b
rho[5] <- sigma[2,4] / sqrt(sigma[2,2] * sigma[4,4])   # Correlation a2 vs f
rho[6] <- sigma[3,4] / sqrt(sigma[3,3] * sigma[4,4])   # Correlation b vs f
   
mean.a[1] ~ dunif(0, 1)
mean.a[2] ~ dunif(0, 1)
mean.b ~ dunif(0, 0.98)  
mean.f ~ dunif(0, 5)
mu[1] <- logit(mean.a[1])
mu[2] <- logit(mean.a[2])
mu[3] <- logit(mean.b)
mu[4] <- log(mean.f)


# Resighting probabilities
# Resighting probabilities of birds with wing tags 
for (t in 1:(n.occasions1-1)){
  pN1[t] <- period.pN[period[t]]
  pN2[t] <- period.pN[period[t]]
  pF1[t] <- period.pF[period[t]]
  pBss1[t] <- period.pBss[period[t]]
  pBns1[t] <- period.pBns[period[t]]
  pBss2[t] <- period.pBss[period[t]]
  pBns2[t] <- period.pBns[period[t]]
}
for (pe in 1:2){
  period.pN[pe] ~ dunif(0, 1)
  period.pF[pe] ~ dunif(0, 1)
  period.pBss[pe] ~ dunif(0, 1)
  period.pBns[pe] <- period.pF[pe]
}

# Resighting probabilities of birds with satellite transmitters
# Relocation probability of non-breeders
for (t in 1:(n.occasions3-1)){
  pN3[t] <- pNs[t+start3-1]
}
pN4 <- pNs[33]

# Relocation probability of breeders
for (t in 1:(n.occasions3-1)){
  pB3[t] <- pBs[t+start3-1]
}

# Apply linear models for relocation probabilities  
for (t in 1:(nyears-1)){
  pNs[t] <- mean.ps
  pBs[t] <- mean.ps
}
mean.ps ~ dunif(0, 1)


# Recovery probabilities
# Recovery probabilities of birds with satellite transmitters
for (t in 1:(n.occasions3-1)){
  ra3[t] <- ras[t+start3-1]
}
rj4 <- ras[33]   # Recovery in year 2018
for (t in 1:(nyears-1)){
  ras[t] <- mean.ras
}
mean.ras ~ dunif(0, 1)

# Recovery probabilities of birds with wing tags
# Adult ring recovery of birds with wing tag 
for (t in 1:(n.occasions1-1)){
  ra1[t] <- ilogit(logit.r[4,t+start1-1])
}
for (t in 1:(n.occasions2-1)){
  ra2[t] <- ilogit(logit.r[4,t+start2-1])
}
for (t in 1:n.occasions5){
  ra5[2,t] <- ilogit(logit.r[4,t+start5-1])
}

# Juvenile ring recovery of birds with wing tag 
for (t in 1:(n.occasions1-1)){
  rj1[t] <- ilogit(logit.r[3,t+start1-1])
}
for (t in 1:n.occasions5){
  rj5[2,t] <- ilogit(logit.r[3,t+start5-1])
}

# Recovery probabilities of birds with rings only
for (t in 1:n.occasions5){
  rj5[1,t] <- ilogit(logit.r[1,t+start5-1])
  ra5[1,t] <- ilogit(logit.r[2,t+start5-1])
}

# Apply linear models for the different recovery probabilities
for (t in 1:(nyears-1)){
  logit.r[1:4,t] ~ dmnorm(mur[], taur[,])
  for (j in 1:4){
    r[j,t] <- ilogit(logit.r[j,t])
  } #j
} #t
for (i in 1:4){
  mean.r[i] ~ dunif(0, 0.5)
  mur[i] <- logit(mean.r[i])
}
taur[1:4,1:4] ~ dwish(W[,], 5)
sigmar[1:4,1:4] <- inverse(taur[1:4,1:4])
rhor[1] <- sigmar[1,2] / sqrt(sigmar[1,1] * sigmar[2,2])   # Correlation r1 vs r2
rhor[2] <- sigmar[1,3] / sqrt(sigmar[1,1] * sigmar[3,3])   # Correlation r1 vs r3
rhor[3] <- sigmar[1,4] / sqrt(sigmar[1,1] * sigmar[4,4])   # Correlation r1 vs r4
rhor[4] <- sigmar[2,3] / sqrt(sigmar[2,2] * sigmar[3,3])   # Correlation r2 vs r3
rhor[5] <- sigmar[2,4] / sqrt(sigmar[2,2] * sigmar[4,4])   # Correlation r2 vs r4
rhor[6] <- sigmar[3,4] / sqrt(sigmar[3,3] * sigmar[4,4])   # Correlation r3 vs r4


# Observation error for the count data
sigma.obs ~ dunif(0.001, 10)
tau.obs <- pow(sigma.obs, -2)


# Initial population sizes  (means according to stable age distribution)
N[1,1] ~ dcat(Ninit1)     # Total number of local 1y
N[2,1] ~ dcat(Ninit2)     # 2-years first time breeders 
N[3,1] ~ dcat(Ninit2)     # 2-years non-breeder
N[4,1] ~ dcat(Ninit2)     # 3-years first time breeders
N[5,1] ~ dcat(Ninit2)     # 3-years non-breeder
N[6,1] ~ dcat(Ninit2)     # 4-year first time breeders
N[7,1] ~ dcat(Ninit2)     # Experienced breeders
N[8,1] ~ dcat(Ninit2)     # Immigrants

N.pred[1,1] ~ dcat(Ninit1)     # Total number of local 1y
N.pred[2,1] ~ dcat(Ninit2)     # 2-years first time breeders 
N.pred[3,1] ~ dcat(Ninit2)     # 2-years non-breeder
N.pred[4,1] ~ dcat(Ninit2)     # 3-years first time breeders
N.pred[5,1] ~ dcat(Ninit2)     # 3-years non-breeder
N.pred[6,1] ~ dcat(Ninit2)     # 4-year first time breeders
N.pred[7,1] ~ dcat(Ninit2)     # Experienced breeders
N.pred[8,1] ~ dcat(Ninit2)     # Immigrants


#-------------------------------------------------
# 2. The likelihoods of the single data sets
#-------------------------------------------------
# 2.1. Likelihood for population population count data (state-space model)
# 2.1.1 System process
for (t in 2:nyears){
  N[1,t] ~ dpois(s[1,t-1] * fi[t-1] * f[t-1]/2 * (N[2,t-1] + N[4,t-1] + N[6,t-1] + N[8,t-1]) + s[1,t-1] * fi[t-1] * f[t-1]/2 * b[t] * N[7,t-1])         
  N[2,t] ~ dbin(s[2,t-1] * a[1,t-1], N[1,t-1])                                
  N[3,t] ~ dbin(s[2,t-1] * (1-a[1,t-1]), N[1,t-1])                                
  N[4,t] ~ dbin(s[3,t-1] * a[2,t-1], N[3,t-1])                                
  N[5,t] ~ dbin(s[3,t-1] * (1-a[2,t-1]), N[3,t-1])                                
  N[6,t] ~ dbin(s[3,t-1], N[5,t-1])                                
  N[7,t] ~ dbin(s[3,t-1], (N[2,t-1] + N[4,t-1] + N[6,t-1] + N[7,t-1] + N[8,t-1]))       
  N[8,t] ~ dpois((N[2,t-1] + N[4,t-1] + N[6,t-1] + b[t-1]*N[7,t-1] + N[8,t-1]) * im[t-1])     # Immigration rate with respect to local number of breeders in previous year
}      

for (t in 1:nyears){
  Np[t] <- N[2,t] + N[4,t] + N[6,t] + N[7,t] * b[t] + N[8,t]   # Number of breeders in a given year
}

# 2.1.2 Observation process
for (t in 1:nyears){
  counts[t] ~ dnorm(Np[t], tau.obs)
  counts.pred[t] ~ dnorm(Np[t], tau.obs)
}

# 2.1.3 GOF for count data: mean absolute percentage error
for (t in 1:nyears){
  disc.C[t] <- pow(((counts[t] - Np[t]) / counts[t]) * ((counts[t] - Np[t]) / counts[t]), 0.5)
  discN.C[t] <- pow(((counts.pred[t] - Np[t]) / (counts.pred[t] + 0.001)) * ((counts.pred[t] - Np[t]) / (counts.pred[t] + 0.001)), 0.5)
}
fit.C <- 100 / nyears * sum(disc.C)
fitN.C <- 100 / nyears * sum(discN.C)


# 2.2. Likelihood for data on reproductive success
for (i in 1:nbr){
  Juv[i] ~ dpois(f[Jahr[i]])   
  delta[i] <- step(Juv[i] - censor[i])
  d[i] ~ dbern(delta[i])
   
  # Prediction for GOF
  # Simulate new data
  Juv.pred[i] ~ dpois(f[Jahr[i]])
  # Expected values
  Juv.exp[i] <- f[Jahr[i]]
   
  # Discrepancy meassure
  dev[i] <- Juv[i] * log(Juv[i]/Juv.exp[i])
  devN[i] <- Juv.pred[i] * log(Juv.pred[i]/Juv.exp[i])      
}

fit.dev <- 2*sum(dev)
fitN.dev <- 2*sum(devN)


# 2.3. Ringing data from Weimar
# 2.3.1. Ringed as juv with 'Draht' wing tag
# -------------------------------------------------
# Parameters:
# s: survival probability
# r: recovery probability
# a: probability of first time reproduction
# fi: site fidelity
# b: breeding probability of experienced breeder
# pF: probability to resight a first-time breeder
# pBss: probability to resight a breeder that has been seen at the previous occasion
# pBns: probability to resight a breeder that has not been seen at the previous occasion
# pN: probability to resight a non-breeder
# -------------------------------------------------
# States:
# 1 juv 
# 2 1y non-breeder inside
# 3 2y non-breeder inside
# 4 3y non-breeder inside
# 5 Ad non-breeder inside
# 6 Ad breeder inside, seen
# 7 recently dead Juv (and recovered)
# 8 recently dead Ad (and recovered)
# 9 1y indv outside
# 10 2y+ indv outside
# 11 Ad breeder inside, not seen
# 12 dead (not recovered)
# -------------------------------------------------

# Define state-transition and reencounter probabilities
for (t in 1:(n.occasions1-1)){
  psi1[1,t,1] <- 0
  psi1[1,t,2] <- s1[1,t]*fi1[t]
  psi1[1,t,3] <- 0
  psi1[1,t,4] <- 0
  psi1[1,t,5] <- 0
  psi1[1,t,6] <- 0
  psi1[1,t,7] <- 1-s1[1,t]
  psi1[1,t,8] <- 0
  psi1[1,t,9] <- s1[1,t]*(1-fi1[t])
  psi1[1,t,10] <- 0
  psi1[1,t,11] <- 0
  psi1[1,t,12] <- 0

  psi1[2,t,1] <- 0
  psi1[2,t,2] <- 0
  psi1[2,t,3] <- s1[2,t]*(1-a1[1,t])
  psi1[2,t,4] <- 0
  psi1[2,t,5] <- 0
  psi1[2,t,6] <- s1[2,t]*a1[1,t]*pBns1[t]
  psi1[2,t,7] <- 0
  psi1[2,t,8] <- 1-s1[2,t]
  psi1[2,t,9] <- 0
  psi1[2,t,10] <- 0
  psi1[2,t,11] <- s1[2,t]*a1[1,t]*(1-pBns1[t])
  psi1[2,t,12] <- 0

  psi1[3,t,1] <- 0
  psi1[3,t,2] <- 0
  psi1[3,t,3] <- 0
  psi1[3,t,4] <- s1[3,t]*(1-a1[2,t])
  psi1[3,t,5] <- 0
  psi1[3,t,6] <- s1[3,t]*a1[2,t]*pBns1[t]
  psi1[3,t,7] <- 0
  psi1[3,t,8] <- 1-s1[3,t]
  psi1[3,t,9] <- 0
  psi1[3,t,10] <- 0
  psi1[3,t,11] <- s1[3,t]*a1[2,t]*(1-pBns1[t])
  psi1[3,t,12] <- 0

  psi1[4,t,1] <- 0
  psi1[4,t,2] <- 0
  psi1[4,t,3] <- 0
  psi1[4,t,4] <- 0
  psi1[4,t,5] <- 0
  psi1[4,t,6] <- s1[3,t]*pBns1[t]
  psi1[4,t,7] <- 0
  psi1[4,t,8] <- 1-s1[3,t]
  psi1[4,t,9] <- 0
  psi1[4,t,10] <- 0
  psi1[4,t,11] <- s1[3,t]*(1-pBns1[t])
  psi1[4,t,12] <- 0

  psi1[5,t,1] <- 0
  psi1[5,t,2] <- 0
  psi1[5,t,3] <- 0
  psi1[5,t,4] <- 0
  psi1[5,t,5] <- s1[3,t]*(1-b1[t])
  psi1[5,t,6] <- s1[3,t]*b1[t]*pBss1[t]
  psi1[5,t,7] <- 0
  psi1[5,t,8] <- 1-s1[3,t]
  psi1[5,t,9] <- 0
  psi1[5,t,10] <- 0
  psi1[5,t,11] <- s1[3,t]*b1[t]*(1-pBss1[t])
  psi1[5,t,12] <- 0

  psi1[6,t,1] <- 0
  psi1[6,t,2] <- 0
  psi1[6,t,3] <- 0
  psi1[6,t,4] <- 0
  psi1[6,t,5] <- s1[3,t]*(1-b1[t])
  psi1[6,t,6] <- s1[3,t]*b1[t]*pBss1[t]
  psi1[6,t,7] <- 0
  psi1[6,t,8] <- 1-s1[3,t]
  psi1[6,t,9] <- 0
  psi1[6,t,10] <- 0
  psi1[6,t,11] <- s1[3,t]*b1[t]*(1-pBss1[t])
  psi1[6,t,12] <- 0

  psi1[7,t,1] <- 0
  psi1[7,t,2] <- 0
  psi1[7,t,3] <- 0
  psi1[7,t,4] <- 0
  psi1[7,t,5] <- 0
  psi1[7,t,6] <- 0
  psi1[7,t,7] <- 0
  psi1[7,t,8] <- 0
  psi1[7,t,9] <- 0
  psi1[7,t,10] <- 0
  psi1[7,t,11] <- 0
  psi1[7,t,12] <- 1

  psi1[8,t,1] <- 0
  psi1[8,t,2] <- 0
  psi1[8,t,3] <- 0
  psi1[8,t,4] <- 0
  psi1[8,t,5] <- 0
  psi1[8,t,6] <- 0
  psi1[8,t,7] <- 0
  psi1[8,t,8] <- 0
  psi1[8,t,9] <- 0
  psi1[8,t,10] <- 0
  psi1[8,t,11] <- 0
  psi1[8,t,12] <- 1

  psi1[9,t,1] <- 0
  psi1[9,t,2] <- 0
  psi1[9,t,3] <- 0
  psi1[9,t,4] <- 0
  psi1[9,t,5] <- 0
  psi1[9,t,6] <- 0
  psi1[9,t,7] <- 0
  psi1[9,t,8] <- 1-s1[2,t]
  psi1[9,t,9] <- 0
  psi1[9,t,10] <- s1[2,t]
  psi1[9,t,11] <- 0
  psi1[9,t,12] <- 0

  psi1[10,t,1] <- 0
  psi1[10,t,2] <- 0
  psi1[10,t,3] <- 0
  psi1[10,t,4] <- 0
  psi1[10,t,5] <- 0
  psi1[10,t,6] <- 0
  psi1[10,t,7] <- 0
  psi1[10,t,8] <- 1-s1[3,t]
  psi1[10,t,9] <- 0
  psi1[10,t,10] <- s1[3,t]
  psi1[10,t,11] <- 0
  psi1[10,t,12] <- 0

  psi1[11,t,1] <- 0
  psi1[11,t,2] <- 0
  psi1[11,t,3] <- 0
  psi1[11,t,4] <- 0
  psi1[11,t,5] <- s1[3,t]*(1-b1[t])
  psi1[11,t,6] <- s1[3,t]*b1[t]*pBns1[t]
  psi1[11,t,7] <- 0
  psi1[11,t,8] <- 0
  psi1[11,t,9] <- 0
  psi1[11,t,10] <- 0
  psi1[11,t,11] <- s1[3,t]*b1[t]*(1-pBns1[t])
  psi1[11,t,12] <- 0

  psi1[12,t,1] <- 0
  psi1[12,t,2] <- 0
  psi1[12,t,3] <- 0
  psi1[12,t,4] <- 0
  psi1[12,t,5] <- 0
  psi1[12,t,6] <- 0
  psi1[12,t,7] <- 0
  psi1[12,t,8] <- 0
  psi1[12,t,9] <- 0
  psi1[12,t,10] <- 0
  psi1[12,t,11] <- 0
  psi1[12,t,12] <- 1
 
  po1[1,t] <- 0
  po1[2,t] <- pF1[t]
  po1[3,t] <- pF1[t]
  po1[4,t] <- pF1[t]
  po1[5,t] <- pN1[t]
  po1[6,t] <- 1
  po1[7,t] <- rj1[t]
  po1[8,t] <- ra1[t]
  po1[9,t] <- 0
  po1[10,t] <- 0
  po1[11,t] <- 0
  po1[12,t] <- 0


  # Calculate probability of non-encounter (dq) and reshape the array for the encounter probabilities      
  for (s in 1:ns1){
    dp1[s,t,s] <- po1[s,t]
    dq1[s,t,s] <- 1-po1[s,t]
  } #s
  for (s in 1:(ns1-1)){
    for (m in (s+1):ns1){
      dp1[s,t,m] <- 0
      dq1[s,t,m] <- 0
    } #s
  } #m
  for (s in 2:ns1){
    for (m in 1:(s-1)){
      dp1[s,t,m] <- 0
      dq1[s,t,m] <- 0
    } #s
  } #m
} #t

# Define the multinomial likelihood
for (t in 1:((n.occasions1-1)*ns1)){
  marr1[t,1:(n.occasions1*ns1-(ns1-1))] ~ dmulti(pr1[t,], rel1[t])
  marr1N[t,1:(n.occasions1*ns1-(ns1-1))] ~ dmulti(pr1[t,], rel1[t])
}

# Define the cell probabilities of the multistate m-array   
# Define matrix U: product of probabilities of state-transition and non-encounter (this is just done because there is no product function for matrix multiplication in JAGS)
for (t in 1:(n.occasions1-2)){
  U1[(t-1)*ns1+(1:ns1), (t-1)*ns1+(1:ns1)] <- ones1
  for (j in (t+1):(n.occasions1-1)){
    U1[(t-1)*ns1+(1:ns1), (j-1)*ns1+(1:ns1)] <- U1[(t-1)*ns1+(1:ns1), (j-2)*ns1+(1:ns1)] %*% psi1[,t,] %*% dq1[,t,]
  } #j
} #t
U1[(n.occasions1-2)*ns1+(1:ns1), (n.occasions1-2)*ns1+(1:ns1)] <- ones1

# Diagonal
for (t in 1:(n.occasions1-2)){
  pr1[(t-1)*ns1+(1:ns1),(t-1)*ns1+(1:ns1)] <- U1[(t-1)*ns1+(1:ns1),(t-1)*ns1+(1:ns1)] %*% psi1[,t,] %*% dp1[,t,]
  # Above main diagonal
  for (j in (t+1):(n.occasions1-1)){
    pr1[(t-1)*ns1+(1:ns1), (j-1)*ns1+(1:ns1)] <- U1[(t-1)*ns1+(1:ns1), (j-1)*ns1+(1:ns1)] %*% psi1[,j,] %*% dp1[,j,]
  } #j
} #t
pr1[(n.occasions1-2)*ns1+(1:ns1), (n.occasions1-2)*ns1+(1:ns1)] <- psi1[,n.occasions1-1,] %*% dp1[,n.occasions1-1,]

# Below main diagonal
for (t in 2:(n.occasions1-1)){
  for (j in 1:(t-1)){
    pr1[(t-1)*ns1+(1:ns1),(j-1)*ns1+(1:ns1)] <- zero1
  } #j
} #t

# Last column: probability of non-recapture
for (t in 1:((n.occasions1-1)*ns1)){
  pr1[t,(n.occasions1*ns1-(ns1-1))] <- 1-sum(pr1[t,1:((n.occasions1-1)*ns1)])
}


# Generate expected m-arrays and Freeman-Tukey test statistics for GOF
for (t in 1:((n.occasions1-1)*ns1)){
  for (j in 1:(n.occasions1*ns1-(ns1-1))){
    marr1.E[t,j] <- pr1[t,j] * rel1[t]
    E.org1[t,j] <- pow((pow(marr1[t,j], 0.5) - pow(marr1.E[t,j], 0.5)), 2)
    E.new1[t,j] <- pow((pow(marr1N[t,j], 0.5) - pow(marr1.E[t,j], 0.5)), 2)
  } #j
} #t

fit.1 <- sum(E.org1[,])
fitN.1 <- sum(E.new1[,])


# 2.3.2. Marked as adult with wing tag 'Draht'
# -------------------------------------------------
# Parameters:
# s: survival probability
# r: recovery probability
# b: breeding probability of experienced breeder
# pBss: probability to resight a breeder that has been seen at the previous occasion
# pBns: probability to resight a breeder that has not been seen at the previous occasion
# pN: probability to resight a non-breeder
# -------------------------------------------------
# States:
# 1 Breeding, seen 
# 2 Non-breeding, seen
# 3 recently dead (and recovered)
# 4 Breeding, not seen
# 5 Non-breeding, not seen
# 6 dead (not recovered)
# -------------------------------------------------

# Define state-transition and observation matrices
# Define probabilities of state S(t+1) given S(t)
for (t in 1:(n.occasions2-1)){
  psi2[1,t,1] <- b2[t]*s2[3,t]*pBss2[t]
  psi2[1,t,2] <- (1-b2[t])*s2[3,t]*pN2[t]
  psi2[1,t,3] <- 1-s2[3,t]
  psi2[1,t,4] <- b2[t]*s2[3,t]*(1-pBss2[t])
  psi2[1,t,5] <- (1-b2[t])*s2[3,t]*(1-pN2[t])
  psi2[1,t,6] <- 0

  psi2[2,t,1] <- b2[t]*s2[3,t]*pBss2[t]
  psi2[2,t,2] <- (1-b2[t])*s2[3,t]*pN2[t]
  psi2[2,t,3] <- 1-s2[3,t]
  psi2[2,t,4] <- b2[t]*s2[3,t]*(1-pBss2[t])
  psi2[2,t,5] <- (1-b2[t])*s2[3,t]*(1-pN2[t])
  psi2[2,t,6] <- 0

  psi2[3,t,1] <- 0
  psi2[3,t,2] <- 0
  psi2[3,t,3] <- 0
  psi2[3,t,4] <- 0
  psi2[3,t,5] <- 0
  psi2[3,t,6] <- 1

  psi2[4,t,1] <- b2[t]*s2[3,t]*pBns2[t]
  psi2[4,t,2] <- (1-b2[t])*s2[3,t]*pN2[t]
  psi2[4,t,3] <- 1-s2[3,t]
  psi2[4,t,4] <- b2[t]*s2[3,t]*(1-pBns2[t])
  psi2[4,t,5] <- (1-b2[t])*s2[3,t]*(1-pN2[t])
  psi2[4,t,6] <- 0

  psi2[5,t,1] <- b2[t]*s2[3,t]*pBns2[t]
  psi2[5,t,2] <- (1-b2[t])*s2[3,t]*pN2[t]
  psi2[5,t,3] <- 1-s2[3,t]
  psi2[5,t,4] <- b2[t]*s2[3,t]*(1-pBns2[t])
  psi2[5,t,5] <- (1-b2[t])*s2[3,t]*(1-pN2[t])
  psi2[5,t,6] <- 0

  psi2[6,t,1] <- 0
  psi2[6,t,2] <- 0
  psi2[6,t,3] <- 0
  psi2[6,t,4] <- 0
  psi2[6,t,5] <- 0
  psi2[6,t,6] <- 1

  po2[1,t] <- 1
  po2[2,t] <- 1
  po2[3,t] <- ra2[t]
  po2[4,t] <- 0
  po2[5,t] <- 0
  po2[6,t] <- 0


  # Calculate probability of non-encounter (dq) and reshape the array for the encounter probabilities      
  for (s in 1:ns2){
    dp2[s,t,s] <- po2[s,t]
    dq2[s,t,s] <- 1-po2[s,t]
  } #s
  for (s in 1:(ns2-1)){
    for (m in (s+1):ns2){
      dp2[s,t,m] <- 0
      dq2[s,t,m] <- 0
    } #s
  } #m
  for (s in 2:ns2){
    for (m in 1:(s-1)){
      dp2[s,t,m] <- 0
      dq2[s,t,m] <- 0
    } #s
  } #m
} #t

# Define the multinomial likelihood
for (t in 1:((n.occasions2-1)*ns2)){
  marr2[t,1:(n.occasions2*ns2-(ns2-1))] ~ dmulti(pr2[t,], rel2[t])
  marr2N[t,1:(n.occasions2*ns2-(ns2-1))] ~ dmulti(pr2[t,], rel2[t])
}

# Define the cell probabilities of the multistate m-array   
# Define matrix U: product of probabilities of state-transition and non-encounter (this is just done because there is no product function for matrix multiplication in JAGS)
for (t in 1:(n.occasions2-2)){
  U2[(t-1)*ns2+(1:ns2), (t-1)*ns2+(1:ns2)] <- ones2
  for (j in (t+1):(n.occasions2-1)){
    U2[(t-1)*ns2+(1:ns2), (j-1)*ns2+(1:ns2)] <- U2[(t-1)*ns2+(1:ns2), (j-2)*ns2+(1:ns2)] %*% psi2[,t,] %*% dq2[,t,]
  } #j
} #t
U2[(n.occasions2-2)*ns2+(1:ns2), (n.occasions2-2)*ns2+(1:ns2)] <- ones2

# Diagonal
for (t in 1:(n.occasions2-2)){
  pr2[(t-1)*ns2+(1:ns2),(t-1)*ns2+(1:ns2)] <- U2[(t-1)*ns2+(1:ns2),(t-1)*ns2+(1:ns2)] %*% psi2[,t,] %*% dp2[,t,]
  # Above main diagonal
  for (j in (t+1):(n.occasions2-1)){
    pr2[(t-1)*ns2+(1:ns2), (j-1)*ns2+(1:ns2)] <- U2[(t-1)*ns2+(1:ns2), (j-1)*ns2+(1:ns2)] %*% psi2[,j,] %*% dp2[,j,]
  } #j
} #t
pr2[(n.occasions2-2)*ns2+(1:ns2), (n.occasions2-2)*ns2+(1:ns2)] <- psi2[,n.occasions2-1,] %*% dp2[,n.occasions2-1,]

# Below main diagonal
for (t in 2:(n.occasions2-1)){
  for (j in 1:(t-1)){
    pr2[(t-1)*ns2+(1:ns2),(j-1)*ns2+(1:ns2)] <- zero2
  } #j
} #t

# Last column: probability of non-recapture
for (t in 1:((n.occasions2-1)*ns2)){
  pr2[t,(n.occasions2*ns2-(ns2-1))] <- 1-sum(pr2[t,1:((n.occasions2-1)*ns2)])
}


# Generate expected m-arrays and Freeman-Tukey test statistics for GOF
for (t in 1:((n.occasions2-1)*ns2)){
  for (j in 1:(n.occasions2*ns2-(ns2-1))){
    marr2.E[t,j] <- pr2[t,j] * rel2[t]
    E.org2[t,j] <- pow((pow(marr2[t,j], 0.5) - pow(marr2.E[t,j], 0.5)), 2)
    E.new2[t,j] <- pow((pow(marr2N[t,j], 0.5) - pow(marr2.E[t,j], 0.5)), 2)
  } #j
} #t

fit.2 <- sum(E.org2[,])
fitN.2 <- sum(E.new2[,])


# 2.3.3. Marked as adult with satellite transmitter 
# -------------------------------------------------
# Parameters:
# s: survival probability
# r: recovery probability
# b: breeding probability of experienced breeder
# pB: probability to resight a breeder 
# pN: probability to resight a non-breeder 
# -------------------------------------------------
# States:
# 1 Breeding 
# 2 Non-breeding
# 3 recently dead (and recovered)
# 4 dead (not recovered)
# -------------------------------------------------

# Define state-transition and observation matrices
# Define probabilities of state S(t+1) given S(t)
for (t in 1:(n.occasions3-1)){
  psi3[1,t,1] <- b3[t]*s3[3,t]
  psi3[1,t,2] <- (1-b3[t])*s3[3,t]
  psi3[1,t,3] <- 1-s3[3,t]
  psi3[1,t,4] <- 0

  psi3[2,t,1] <- b3[t]*s3[3,t]
  psi3[2,t,2] <- (1-b3[t])*s3[3,t]
  psi3[2,t,3] <- 1-s3[3,t]
  psi3[2,t,4] <- 0

  psi3[3,t,1] <- 0
  psi3[3,t,2] <- 0
  psi3[3,t,3] <- 0
  psi3[3,t,4] <- 1

  psi3[4,t,1] <- 0
  psi3[4,t,2] <- 0
  psi3[4,t,3] <- 0
  psi3[4,t,4] <- 1

  po3[1,t] <- pB3[t]
  po3[2,t] <- pN3[t]
  po3[3,t] <- ra3[t]
  po3[4,t] <- 0

  # Calculate probability of non-encounter (dq) and reshape the array for the encounter probabilities      
  for (s in 1:ns3){
    dp3[s,t,s] <- po3[s,t]
    dq3[s,t,s] <- 1-po3[s,t]
  } #s
  for (s in 1:(ns3-1)){
    for (m in (s+1):ns3){
      dp3[s,t,m] <- 0
      dq3[s,t,m] <- 0
    } #s
  } #m
  for (s in 2:ns3){
    for (m in 1:(s-1)){
      dp3[s,t,m] <- 0
      dq3[s,t,m] <- 0
    } #s
  } #m
} #t

# Define the multinomial likelihood
for (t in 1:((n.occasions3-1)*ns3)){
  marr3[t,1:(n.occasions3*ns3-(ns3-1))] ~ dmulti(pr3[t,], rel3[t])
  marr3N[t,1:(n.occasions3*ns3-(ns3-1))] ~ dmulti(pr3[t,], rel3[t])
}

# Define the cell probabilities of the multistate m-array   
# Define matrix U: product of probabilities of state-transition and non-encounter (this is just done because there is no product function for matrix multiplication in JAGS)
for (t in 1:(n.occasions3-2)){
  U3[(t-1)*ns3+(1:ns3), (t-1)*ns3+(1:ns3)] <- ones3
  for (j in (t+1):(n.occasions3-1)){
    U3[(t-1)*ns3+(1:ns3), (j-1)*ns3+(1:ns3)] <- U3[(t-1)*ns3+(1:ns3), (j-2)*ns3+(1:ns3)] %*% psi3[,t,] %*% dq3[,t,]
  } #j
} #t
U3[(n.occasions3-2)*ns3+(1:ns3), (n.occasions3-2)*ns3+(1:ns3)] <- ones3

# Diagonal
for (t in 1:(n.occasions3-2)){
  pr3[(t-1)*ns3+(1:ns3),(t-1)*ns3+(1:ns3)] <- U3[(t-1)*ns3+(1:ns3),(t-1)*ns3+(1:ns3)] %*% psi3[,t,] %*% dp3[,t,]
  # Above main diagonal
  for (j in (t+1):(n.occasions3-1)){
    pr3[(t-1)*ns3+(1:ns3), (j-1)*ns3+(1:ns3)] <- U3[(t-1)*ns3+(1:ns3), (j-1)*ns3+(1:ns3)] %*% psi3[,j,] %*% dp3[,j,]
  } #j
} #t
pr3[(n.occasions3-2)*ns3+(1:ns3), (n.occasions3-2)*ns3+(1:ns3)] <- psi3[,n.occasions3-1,] %*% dp3[,n.occasions3-1,]

# Below main diagonal
for (t in 2:(n.occasions3-1)){
  for (j in 1:(t-1)){
    pr3[(t-1)*ns3+(1:ns3),(j-1)*ns3+(1:ns3)] <- zero3
  } #j
} #t

# Last column: probability of non-recapture
for (t in 1:((n.occasions3-1)*ns3)){
  pr3[t,(n.occasions3*ns3-(ns3-1))] <- 1-sum(pr3[t,1:((n.occasions3-1)*ns3)])
}


# Generate expected m-arrays and Freeman-Tukey test statistics for GOF
for (t in 1:((n.occasions3-1)*ns3)){
  for (j in 1:(n.occasions3*ns3-(ns3-1))){
    marr3.E[t,j] <- pr3[t,j] * rel3[t]
    E.org3[t,j] <- pow((pow(marr3[t,j], 0.5) - pow(marr3.E[t,j], 0.5)), 2)
    E.new3[t,j] <- pow((pow(marr3N[t,j], 0.5) - pow(marr3.E[t,j], 0.5)), 2)
  } #j
} #t

fit.3 <- sum(E.org3[,])
fitN.3 <- sum(E.new3[,])


# 2.3.4. Marked as juvenile with satellite transmitter 
# -------------------------------------------------
# Parameters:
# s: survival probability
# r: recovery probability
# -------------------------------------------------

# Since all the marked birds are identified in the next year, we use a simple logistic regression model

totalRecov[1:4] ~ dmulti(pJuvSat[1:4], totalJuvSat)
totalRecovN[1:4] ~ dmulti(pJuvSat[1:4], totalJuvSat)

pJuvSat[1] <- s4 * fi4 * pN4
pJuvSat[2] <- s4 * (1-fi4) * pN4
pJuvSat[3] <- (1-s4) * rj4
pJuvSat[4] <- (1-s4) * (1-rj4) + s4 * (1-pN4)

# Generate expected values and Freeman-Tukey test statistics for GOF
for (i in 1:4){
  d4.E[i] <- pJuvSat[i] * totalJuvSat
  E.org4[i] <- pow((pow(totalRecov[i], 0.5) - pow(d4.E[i], 0.5)), 2)
  E.new4[i] <- pow((pow(totalRecovN[i], 0.5) - pow(d4.E[i], 0.5)), 2)
}

fit.4 <- sum(E.org4)
fitN.4 <- sum(E.new4)



# 2.3.5. Ring-recovery data from Hiddensee
# Define the multinomial likelihoods
for (t in 1:n.occasions5){
  marr.jr[t,1:(n.occasions5+1)] ~ dmulti(pr.jr[t,], rel.jr[t])
  marr.jw[t,1:(n.occasions5+1)] ~ dmulti(pr.jw[t,], rel.jw[t])
  marr.1r[t,1:(n.occasions5+1)] ~ dmulti(pr.1r[t,], rel.1r[t])
  marr.1w[t,1:(n.occasions5+1)] ~ dmulti(pr.1w[t,], rel.1w[t])
  marr.ar[t,1:(n.occasions5+1)] ~ dmulti(pr.ar[t,], rel.ar[t])
  marr.aw[t,1:(n.occasions5+1)] ~ dmulti(pr.aw[t,], rel.aw[t])
  # For GOF
  marr.jrN[t,1:(n.occasions5+1)] ~ dmulti(pr.jr[t,], rel.jr[t])
  marr.jwN[t,1:(n.occasions5+1)] ~ dmulti(pr.jw[t,], rel.jw[t])
  marr.1rN[t,1:(n.occasions5+1)] ~ dmulti(pr.1r[t,], rel.1r[t])
  marr.1wN[t,1:(n.occasions5+1)] ~ dmulti(pr.1w[t,], rel.1w[t])
  marr.arN[t,1:(n.occasions5+1)] ~ dmulti(pr.ar[t,], rel.ar[t])
  marr.awN[t,1:(n.occasions5+1)] ~ dmulti(pr.aw[t,], rel.aw[t])   
}

# Define the cell probabilities of the m-array
# Main diagonal
for (t in 1:n.occasions5){
  pr.jr[t,t] <- (1-s5[1,t])*rj5[1,t]
  pr.jw[t,t] <- (1-s5[1,t])*rj5[2,t]
  pr.1r[t,t] <- (1-s5[2,t])*ra5[1,t]
  pr.1w[t,t] <- (1-s5[2,t])*ra5[2,t]
  pr.ar[t,t] <- (1-s5[3,t])*ra5[1,t]
  pr.aw[t,t] <- (1-s5[3,t])*ra5[2,t]
  
  # Below main diagonal
  for (j in 1:(t-1)){
    pr.jr[t,j] <- 0
    pr.jw[t,j] <- 0
    pr.1r[t,j] <- 0
    pr.1w[t,j] <- 0
    pr.ar[t,j] <- 0
    pr.aw[t,j] <- 0
  } #j
} #t
for (t in 1:(n.occasions5-1)){
  # One above main diagonal
  pr.jr[t,t+1] <- s5[1,t]*(1-s5[2,t+1])*ra5[1,t+1] 
  pr.jw[t,t+1] <- s5[1,t]*(1-s5[2,t+1])*ra5[2,t+1] 
  pr.1r[t,t+1] <- s5[2,t]*(1-s5[3,t+1])*ra5[1,t+1] 
  pr.1w[t,t+1] <- s5[2,t]*(1-s5[3,t+1])*ra5[2,t+1] 
  pr.ar[t,t+1] <- s5[3,t]*(1-s5[3,t+1])*ra5[1,t+1] 
  pr.aw[t,t+1] <- s5[3,t]*(1-s5[3,t+1])*ra5[2,t+1] 
}
for (t in 1:(n.occasions5-2)){
  # Two above main diagonal
  pr.jr[t,t+2] <- s5[1,t]*s5[2,t+1]*(1-s5[3,t+2])*ra5[1,t+2] 
  pr.jw[t,t+2] <- s5[1,t]*s5[2,t+1]*(1-s5[3,t+2])*ra5[2,t+2] 
  pr.1r[t,t+2] <- s5[2,t]*s5[3,t+1]*(1-s5[3,t+2])*ra5[1,t+2] 
  pr.1w[t,t+2] <- s5[2,t]*s5[3,t+1]*(1-s5[3,t+2])*ra5[2,t+2] 
  pr.ar[t,t+2] <- s5[3,t]*s5[3,t+1]*(1-s5[3,t+2])*ra5[1,t+2] 
  pr.aw[t,t+2] <- s5[3,t]*s5[3,t+1]*(1-s5[3,t+2])*ra5[2,t+2] 
}

# Further above main diagonal
for (t in 1:(n.occasions5-3)){
  for (j in (t+3):n.occasions5){
    pr.jr[t,j] <- s5[1,t]*s5[2,t+1]*prod(s5[3,(t+2):(j-1)])*(1-s5[3,j])*ra5[1,j]
    pr.jw[t,j] <- s5[1,t]*s5[2,t+1]*prod(s5[3,(t+2):(j-1)])*(1-s5[3,j])*ra5[2,j]      
    pr.1r[t,j] <- s5[2,t]*s5[3,t+1]*prod(s5[3,(t+2):(j-1)])*(1-s5[3,j])*ra5[1,j]
    pr.1w[t,j] <- s5[2,t]*s5[3,t+1]*prod(s5[3,(t+2):(j-1)])*(1-s5[3,j])*ra5[2,j]
    pr.ar[t,j] <- s5[3,t]*s5[3,t+1]*prod(s5[3,(t+2):(j-1)])*(1-s5[3,j])*ra5[1,j]
    pr.aw[t,j] <- s5[3,t]*s5[3,t+1]*prod(s5[3,(t+2):(j-1)])*(1-s5[3,j])*ra5[2,j]
  } #j
} #t

# Last column: probability of non-recovery
for (t in 1:n.occasions5){
  pr.jr[t,n.occasions5+1] <- 1-sum(pr.jr[t,1:n.occasions5])
  pr.jw[t,n.occasions5+1] <- 1-sum(pr.jw[t,1:n.occasions5])
  pr.1r[t,n.occasions5+1] <- 1-sum(pr.1r[t,1:n.occasions5])
  pr.1w[t,n.occasions5+1] <- 1-sum(pr.1w[t,1:n.occasions5])
  pr.ar[t,n.occasions5+1] <- 1-sum(pr.ar[t,1:n.occasions5])
  pr.aw[t,n.occasions5+1] <- 1-sum(pr.aw[t,1:n.occasions5])
}

# Generate expected m-arrays and Freeman-Tukey test statistics for GOF
for (t in 1:n.occasions5){
  for (j in 1:(n.occasions5+1)){
    marr.jrE[t,j] <- pr.jr[t,j] * rel.jr[t]
    marr.jwE[t,j] <- pr.jw[t,j] * rel.jw[t]
    marr.1rE[t,j] <- pr.1r[t,j] * rel.1r[t]
    marr.1wE[t,j] <- pr.1w[t,j] * rel.1w[t]
    marr.arE[t,j] <- pr.ar[t,j] * rel.ar[t]
    marr.awE[t,j] <- pr.aw[t,j] * rel.aw[t]
    E.org.jr[t,j] <- pow((pow(marr.jr[t,j], 0.5) - pow(marr.jrE[t,j], 0.5)), 2)
    E.org.jw[t,j] <- pow((pow(marr.jw[t,j], 0.5) - pow(marr.jwE[t,j], 0.5)), 2)
    E.org.1r[t,j] <- pow((pow(marr.1r[t,j], 0.5) - pow(marr.1rE[t,j], 0.5)), 2)
    E.org.1w[t,j] <- pow((pow(marr.1w[t,j], 0.5) - pow(marr.1wE[t,j], 0.5)), 2)
    E.org.ar[t,j] <- pow((pow(marr.ar[t,j], 0.5) - pow(marr.arE[t,j], 0.5)), 2)
    E.org.aw[t,j] <- pow((pow(marr.aw[t,j], 0.5) - pow(marr.awE[t,j], 0.5)), 2)
    E.new.jr[t,j] <- pow((pow(marr.jrN[t,j], 0.5) - pow(marr.jrE[t,j], 0.5)), 2)
    E.new.jw[t,j] <- pow((pow(marr.jwN[t,j], 0.5) - pow(marr.jwE[t,j], 0.5)), 2)
    E.new.1r[t,j] <- pow((pow(marr.1rN[t,j], 0.5) - pow(marr.1rE[t,j], 0.5)), 2)
    E.new.1w[t,j] <- pow((pow(marr.1wN[t,j], 0.5) - pow(marr.1wE[t,j], 0.5)), 2)
    E.new.ar[t,j] <- pow((pow(marr.arN[t,j], 0.5) - pow(marr.arE[t,j], 0.5)), 2)
    E.new.aw[t,j] <- pow((pow(marr.awN[t,j], 0.5) - pow(marr.awE[t,j], 0.5)), 2)
  } #j
} #t

fit.jr <- sum(E.org.jr[,])
fit.jw <- sum(E.org.jw[,])
fit.1r <- sum(E.org.1r[,])
fit.1w <- sum(E.org.1w[,])
fit.ar <- sum(E.org.ar[,])
fit.aw <- sum(E.org.aw[,])
fitN.jr <- sum(E.new.jr[,])
fitN.jw <- sum(E.new.jw[,])
fitN.1r <- sum(E.new.1r[,])
fitN.1w <- sum(E.new.1w[,])
fitN.ar <- sum(E.new.ar[,])
fitN.aw <- sum(E.new.aw[,])

fit.5 <- fit.jr + fit.jw + fit.1r + fit.1w + fit.ar + fit.aw  
fitN.5 <- fitN.jr + fitN.jw + fitN.1r + fitN.1w + fitN.ar + fitN.aw 
}
")




########################################################
#
#  5. Fit the IPM with JAGS
#
########################################################

# 5.1. Define initial values
inits <- function(){list(Juv = init.Juv, mean.s = c(runif(1, 0.05, 0.1), runif(1, 0.4, 0.6), runif(1, 0.4, 0.6)), mean.fi = 0.8, mean.im = 0.001,  mean.a = c(runif(1, 0.05, 0.1), runif(1, 0.4, 0.55)), mean.b = runif(1, 0.9, 0.96))}  

# 5.2 Bundle the data
jags.data <- list(marr1 = marr1, rel1 = rowSums(marr1), ns1 = 12, zero1 = diag(12)*0, ones1 = diag(12), 
                  marr2 = marr2, rel2 = rowSums(marr2), ns2 = 6, zero2 = diag(6)*0, ones2 = diag(6),
                  marr3 = marr3, rel3 = rowSums(marr3), ns3 = 4, zero3 = diag(4)*0, ones3 = diag(4),
                  marr.jr = marr.jr, marr.1r = marr.1r, marr.ar = marr.ar, marr.jw = marr.jw, marr.1w = marr.1w, marr.aw = marr.aw, rel.jr = rowSums(marr.jr), rel.1r = rowSums(marr.1r), rel.ar = rowSums(marr.ar), rel.jw = rowSums(marr.jw), rel.1w = rowSums(marr.1w), rel.aw = rowSums(marr.aw), totalJuvSat = totalJuvSat, totalRecov = totalRecov,
                  n.occasions1 = ncol(CHJuvDraht), n.occasions2 = ncol(CHAdDraht), n.occasions3 = ncol(CHAdSat), n.occasions5 = nrow(marr.ar),
                  start1 = min(match(year1, year.total)), start2 = min(match(year2, year.total)), start3 = min(match(year3, year.total)), start5 = 1, 
                  Juv = Juv, Jahr = Jahr-min(Jahr)+1, nyears = nyears, nbr = length(Juv), censor = censor, d = d, counts = counts, Ninit1 = dUnif(1, 60), Ninit2 = dUnif(1, 30), period = period, W = diag(4), Q = diag(3))

# 5.3. Define parameters to be monitored
parameters <- c("mean.s", "mean.b", "mean.a", "mean.fi", "mean.im", "mean.f", "mean.r", "mean.ras", "sigmas", "sigmar", "sigma.im", "sigma", "rhos", "rho", "rhor", "period.pN", "period.pBss", "period.pBns", "mean.ps", "sigma.obs",
"s", "f", "a", "b", "fi", "r", "im", "N", "Np", "mu", "mur", "fit.jr", "fit.jw", "fit.1r", "fit.1w", "fit.ar", "fit.aw", "fitN.jr", "fitN.jw", "fitN.1r", "fitN.1w", "fitN.ar", "fitN.aw", "fit.C", "fitN.C", "fit.dev", "fitN.dev", "fit.1", "fitN.1", "fit.2", "fitN.2", "fit.3", "fitN.3", "fit.4", "fitN.4", "fit.5", "fitN.5")


# 5.4. Define MCMC settings
ni <- 600000; nt <- 200; nb <- 100000; nc <- 4; na <- 5000

# 5.5. Run the model 
mod <- jags(jags.data, inits, parameters, "M1.txt", n.chains=nc, n.thin=nt, n.iter=ni, n.burnin=nb, n.adapt=na, parallel=T)

# 5.6. Save the results (and the used data) 
save(mod, jags.data, file="IPM_Results.Rdata")



########################################################
#
#  6. Define linear mixed model for nestling condition in BUGS language
#
########################################################

# Feature of model 1:
# - Needs to account for brood ID (as some nestlings originate from the same nest)
# - No linear trend

# Specify model in BUGS language
cat(file = "Ezu1.txt", "
model {

for (t in 1:nyears){
  beta[t] ~ dnorm(0, tau.beta)
  pred[t] <- alpha + beta[t]
}
alpha ~ dnorm(0, 0.001)

sigma ~ dunif(0, 100)
tau <- pow(sigma, -2)

sigma.beta ~ dunif(0, 10)
tau.beta <- pow(sigma.beta, -2)

sigma.gamma ~ dunif(0, 20)
tau.gamma <- pow(sigma.gamma, -2)

# Likelihood
for (i in 1:length(y)){
  y[i] ~ dnorm(mu[i], tau)
  mu[i] <- alpha + beta[year[i]] + gamma[broodID[i]]
}

for (i in 1:nID){
  gamma[i] ~ dnorm(0, tau.gamma)   
}
}
")


# Feature of model 2:
# - Needs to account for brood ID (as some nestlings originate from the same nest)
# - Linear trend over time

# Specify model in BUGS language
cat(file = "Ezu2.txt", "
model {

for (t in 1:nyears){
  beta[t] ~ dnorm(0, tau.beta)
  pred[t] <- alpha[1] + alpha[2] * (t-nyears/2) + beta[t]
}
for (i in 1:2){
  alpha[i] ~ dnorm(0, 0.001)
}

sigma ~ dunif(0, 100)
tau <- pow(sigma, -2)

sigma.beta ~ dunif(0, 10)
tau.beta <- pow(sigma.beta, -2)

sigma.gamma ~ dunif(0, 20)
tau.gamma <- pow(sigma.gamma, -2)

# Likelihood
for (i in 1:length(y)){
  y[i] ~ dnorm(mu[i], tau)
  mu[i] <- alpha[1] + alpha[2] * (year[i]-nyears/2) + beta[year[i]] + gamma[broodID[i]]
}

for (i in 1:nID){
  gamma[i] ~ dnorm(0, tau.gamma)   
}
}
")


########################################################
#
#  7. Fit the linear mixed model for nestling condition with JAGS
#
########################################################

# 7.1. Define initial values
inits <- function(){list(sigma.gamma = runif(1, 0, 20))}  

# 7.2. Bundle the data
jags.data <- list(y = Cond$Ezu, broodID = as.numeric(factor(Cond$broodID)), juv = Cond$Juv, year = Cond$Year-min(Cond$Year)+1, nID = length(unique(Cond$broodID)), nyears = length(unique(Cond$Year)))

# 7.3. Define the parameters monitored
parameters <- c("alpha", "beta", "pred", "sigma", "sigma.beta", "sigma.gamma")

# 7.4. Define the MCMC settings (such that we get the same number of samples as in the main analyses)
ni <- 10000; nt <- 2; nb <- 5000; nc <- 4; na <- 2000

# 7.5. Run the models
# Model 1
e1 <- jags(jags.data, inits, parameters, "Ezu1.txt", n.chains = nc, n.thin=nt, n.iter=ni, n.burnin=nb, n.adapt=na, parallel=T)

# Model 2
e2 <- jags(jags.data, inits, parameters, "Ezu2.txt", n.chains = nc, n.thin=nt, n.iter=ni, n.burnin=nb, n.adapt=na, parallel=T)


# 7.6. Save the model results
save(jags.data, e1, e2, file = "NutCond_Results.Rdata")



########################################################
#
#  8. Prospective and retrospective perturbation analyses
#
########################################################

setwd('...')

# 8.1. Load the 'data' (results from the IPM) 
load("IPM_Results.Rdata")
draws <- mod$sims.list     # store MCMC draws in a new object to simplify calculation
n.mcmc <- mod$mcmc.info$n.samples
ny <- dim(mod$mean$N)[2]

# 8.2. Perform tLTRE calculations
# Step 1: Compute realized population growth rates for female red kites. 
# Mean population growth rate (geometric mean)
Ntot <- apply(mod$sims.list$N, c(1,3), sum)           # Total population size
log_lam <- log(Ntot[,2:ny]) - log(Ntot[,1:(ny-1)])    # Realized growth rates (log)

# Mean growth
mlog_lam <- apply(log_lam, 1, mean)
m_lam <- exp(mlog_lam)                 # on natural scale
mean(m_lam)
quantile(m_lam, c(0.025, 0.975))


# Step 2: Calculate normalized abundances for each female age class at each time step and for each of the saved MCMC samples. 
n1 <- draws$N[, 1,] / Ntot
n2 <- draws$N[, 2,] / Ntot
n3 <- draws$N[, 3,] / Ntot
n4 <- draws$N[, 4,] / Ntot
n5 <- draws$N[, 5,] / Ntot
n6 <- draws$N[, 6,] / Ntot
n7 <- draws$N[, 7,] / Ntot
n8 <- draws$N[, 8,] / Ntot


# Step 3: Calculate the transient sensitivities (and elasticities) for each demographic parameter, evaluated at temporal means of each parameter. 

# Compute symbolic sensitivities of observed population growth rate to changes in vital rates and components of population structure.                       
obsgr <- expression(((n2+n4+n6+n8)*0.5*f*s1*fi + n7*0.5*f*s1*fi*b + n1*s2*a1 + n1*s2*(1-a1) + n3*s3*a2 + n3*s3*(1-a2) + n5*s3 + s3*(n2+n4+n6+n7+n8) + im*(n2+n4+n6+b*n7+n8)) / (n1 + n2 + n3 + n4 + n5 + n6 + n7 + n8))

# Extract the mean demographic rates and population sizes and store them in a list
mu <- list(f=draws$mean.f, s1=draws$mean.s[,1], s2=draws$mean.s[,2], s3=draws$mean.s[,3], fi=draws$mean.fi, a1=draws$mean.a[,1], a2=draws$mean.a[,2], b=draws$mean.b, im=draws$mean.im, n1=apply(n1, 1, mean), n2=apply(n2, 1, mean), n3=apply(n3, 1, mean), n4=apply(n4, 1, mean), n5=apply(n5, 1, mean), n6=apply(n6, 1, mean), n7=apply(n7, 1, mean), n8=apply(n8, 1, mean))

# Calculate sensitivities
n.draws <- mod$mcmc.info$n.samples
sens <- matrix(NA, n.draws, 17)  
colnames(sens) <- c("s1", "s2", "s3", "f", "fi", "a1", "a2", "b", "im", "n1", "n2", "n3", "n4", "n5", "n6", "n7", "n8")
sens[,"s1"] <- eval(D(obsgr, "s1"), envir=mu)
sens[,"s2"] <- eval(D(obsgr, "s2"), envir=mu)
sens[,"s3"] <- eval(D(obsgr, "s3"), envir=mu)
sens[,"f"] <- eval(D(obsgr, "f"), envir=mu)
sens[,"fi"] <- eval(D(obsgr, "fi"), envir=mu)
sens[,"a1"] <- eval(D(obsgr, "a1"), envir=mu)
sens[,"a2"] <- eval(D(obsgr, "a2"), envir=mu)
sens[,"b"] <- eval(D(obsgr, "b"), envir=mu)
sens[,"im"] <- eval(D(obsgr, "im"), envir=mu)
sens[,"n1"] <- eval(D(obsgr, "n1"), envir=mu)
sens[,"n2"] <- eval(D(obsgr, "n2"), envir=mu)
sens[,"n3"] <- eval(D(obsgr, "n3"), envir=mu)
sens[,"n4"] <- eval(D(obsgr, "n4"), envir=mu)
sens[,"n5"] <- eval(D(obsgr, "n5"), envir=mu)
sens[,"n6"] <- eval(D(obsgr, "n6"), envir=mu)
sens[,"n7"] <- eval(D(obsgr, "n7"), envir=mu)
sens[,"n8"] <- eval(D(obsgr, "n8"), envir=mu)

# Calculate elasticities
elas <- matrix(NA, n.draws, 17)  
colnames(elas) <- c("s1", "s2", "s3", "f", "fi", "a1", "a2", "b", "im", "n1", "n2", "n3", "n4", "n5", "n6", "n7", "n8")
elas[,"s1"] <- sens[,"s1"] * mu$s1 / m_lam
elas[,"s2"] <- sens[,"s2"] * mu$s2 / m_lam
elas[,"s3"] <- sens[,"s3"] * mu$s3 / m_lam
elas[,"f"] <- sens[,"f"] * mu$f / m_lam
elas[,"fi"] <- sens[,"fi"] * mu$fi / m_lam
elas[,"a1"] <- sens[,"a1"] * mu$a1 / m_lam
elas[,"a2"] <- sens[,"a2"] * mu$a2 / m_lam
elas[,"b"] <- sens[,"b"] * mu$b / m_lam
elas[,"im"] <- sens[,"im"] * mu$im / m_lam
elas[,"n1"] <- sens[,"n1"] * mu$s1 / m_lam
elas[,"n2"] <- sens[,"n2"] * mu$n2 / m_lam
elas[,"n3"] <- sens[,"n3"] * mu$n3 / m_lam
elas[,"n4"] <- sens[,"n4"] * mu$n4 / m_lam
elas[,"n5"] <- sens[,"n5"] * mu$n5 / m_lam
elas[,"n6"] <- sens[,"n6"] * mu$n6 / m_lam
elas[,"n7"] <- sens[,"n7"] * mu$n7 / m_lam
elas[,"n8"] <- sens[,"n8"] * mu$n8 / m_lam


# Define matrix to store results
cont <- matrix(NA, nrow=n.draws, ncol=17) 
colnames(cont) <- c("s1", "s2", "s3", "f", "fi", "a1", "a2", "b", "im", "n1", "n2", "n3", "n4", "n5", "n6", "n7", "n8")

# Calculate contributions for each demographic rate and stage-structured population size at each MCMC draw
for (s in 1:n.draws){
  dp_stoch <- cbind(draws$s[s,1,], draws$s[s,2,], draws$s[s,3,], draws$f[s,1:(ny-1)], draws$fi[s,], draws$a[s,1,], draws$a[s,2,], draws$b[s,1:(ny-1)], draws$im[s,1:(ny-1)], n1[s,1:(ny-1)], n2[s,1:(ny-1)], n3[s,1:(ny-1)], n4[s,1:(ny-1)], n5[s,1:(ny-1)], n6[s,1:(ny-1)], n7[s,1:(ny-1)], n8[s,1:(ny-1)])
  # Derive process variance and covariance among demographic parameters using shrunken estimates of demographic rates and proportional pop. sizes
  dp_varcov <- var(dp_stoch)
  sensvec <- sens[s, ]
  # Calculate demographic contributions
  contmatrix <- dp_varcov * outer(sensvec, sensvec)
  cont[s, ] <- rowSums(contmatrix)
}

# 8.3. Save results
save(cont, elas, sens, log_lam, file='Perturbation_Results.Rdata')



########################################################
#
#  9. Source-sink dynamics: calculations of C and R metrics (Runge et al. 2006, Am. Nat.)
#
########################################################

# 9.1. load data (IPM result file)
setwd('...')
load("IPM_Results.Rdata")
n.mcmc <- mod$mcmc.info$n.samples

# 9.2. Calculate stable stage distribution (based on eigen analysis from the Leslie matrix (A) without immigration)
stand <- matrix(NA, nrow=n.mcmc, ncol=7)
for (i in 1:n.mcmc){
  A <- matrix(c(0, mod$sims.list$mean.f[i]/2*mod$sims.list$mean.s[i,1]*mod$sims.list$mean.fi[i], 0, mod$sims.list$mean.f[i]/2*mod$sims.list$mean.s[i,1]*mod$sims.list$mean.fi[i], 0, mod$sims.list$mean.f[i]/2*mod$sims.list$mean.s[i,1]*mod$sims.list$mean.fi[i], mod$sims.list$mean.f[i]/2*mod$sims.list$mean.s[i,1]*mod$sims.list$mean.b[i]*mod$sims.list$mean.fi[i],  
  mod$sims.list$mean.a[i,1]*mod$sims.list$mean.s[i,2], 0, 0, 0, 0, 0, 0,
  (1-mod$sims.list$mean.a[i,1])*mod$sims.list$mean.s[i,2], 0, 0, 0, 0, 0, 0,
  0, 0, mod$sims.list$mean.a[i,2]*mod$sims.list$mean.s[i,3], 0, 0, 0, 0,
  0, 0, (1-mod$sims.list$mean.a[i,2])*mod$sims.list$mean.s[i,3], 0, 0, 0, 0,
  0, 0, 0, 0, mod$sims.list$mean.s[i,3], 0, 0,
  0, mod$sims.list$mean.s[i,3], 0, mod$sims.list$mean.s[i,3], 0, mod$sims.list$mean.s[i,3], mod$sims.list$mean.s[i,3]), 
  ncol=7, nrow=7, byrow=TRUE)

  z <- which(Re(eigen(A)$values)==max(Re(eigen(A)$values)))	
  revec <- Re(eigen(A)$vectors[,z])
  stand[i,1:7] <- as.numeric(matrix(revec/sum(revec))) # standardised eigenvector
}

# 9.3. Calculate the C metrics
Cmet <- mod$sims.list$mean.s[,2] * stand[,1] + (mod$sims.list$mean.f/2*mod$sims.list$mean.s[,1] + mod$sims.list$mean.s[,3]) * stand[,2] +
mod$sims.list$mean.s[,3] * stand[,3] + (mod$sims.list$mean.f/2*mod$sims.list$mean.s[,1] + mod$sims.list$mean.s[,3]) * stand[,4] + 
mod$sims.list$mean.s[,3] * stand[,5] + (mod$sims.list$mean.f/2*mod$sims.list$mean.s[,1] + mod$sims.list$mean.s[,3]) * stand[,6] + 
(mod$sims.list$mean.f/2*mod$sims.list$mean.s[,1]*mod$sims.list$mean.b + mod$sims.list$mean.s[,3]) * stand[,7]

# 9.4. Calculate the R metrics
Rmet <- mod$sims.list$mean.s[,2] * stand[,1] + (mod$sims.list$mean.f/2*mod$sims.list$mean.s[,1]*mod$sims.list$mean.fi + mod$sims.list$mean.s[,3]) * stand[,2] +
mod$sims.list$mean.s[,3] * stand[,3] + (mod$sims.list$mean.f/2*mod$sims.list$mean.s[,1]*mod$sims.list$mean.fi + mod$sims.list$mean.s[,3]) * stand[,4] + 
mod$sims.list$mean.s[,3] * stand[,5] + (mod$sims.list$mean.f/2*mod$sims.list$mean.s[,1]*mod$sims.list$mean.fi + mod$sims.list$mean.s[,3]) * stand[,6] + 
(mod$sims.list$mean.f/2*mod$sims.list$mean.s[,1]*mod$sims.list$mean.b*mod$sims.list$mean.fi + mod$sims.list$mean.s[,3]) * stand[,7] 

# 9.5. Summarize the calculations (means and 95& CRI)
mean(Cmet)
quantile(Cmet, c(0.025, 0.975))

mean(Rmet)
quantile(Rmet, c(0.025, 0.975))


########################################################
#
#  10. Misc stats
#
########################################################

setwd('...')
load("IPM_Results.Rdata")

# Proportion of floaters among new breeders   
n.mcmc <- mod$mcmc.info$n.samples
ny <- dim(mod$mean$N)[2]

# Compute total population size (Ntot), number of non-breeders (NB), of breeders (Breeders) and first-time breeders (FB)
Ntot <- apply(mod$sims.list$N, c(1,3), sum)     
NB <- apply(mod$sims.list$N[,c(1,3,5),], c(1,3), sum)
Breeders <- apply(mod$sims.list$N[,c(2,4,6,7,8),], c(1,3), sum)
FB <- apply(mod$sims.list$N[,c(2,4,6),], c(1,3), sum)   

# Range of the number of breeding pairs
round(range(apply(Breeders, 2, mean)))

# Range of the total number of individuals
round(range(apply(Ntot, 2, mean)))
  
# Mean population growth rate (geometric mean)
lam_t <- log_lam  <- matrix(NA, n.mcmc, ny-1)
for (t in 1:(ny-1)){
  lam_t[,t] <- Ntot[,t+1] / Ntot[,t]                           # Growth rate based on total popualtion size
  log_lam[,t] <- log(apply(mod$sims.list$N[,,t+1], 1, sum)) - log(apply(mod$sims.list$N[,,t], 1, sum))      # for mean geometric mean   
}

# Compute geometric mean and SD of growth rate   
exp(mean(apply(log_lam, 1, mean)))
quantile(exp(apply(log_lam, 1, mean)), prob = c(0.025, 0.975))

# Proportion of non-breeders vs. breeders 
NBratio <- matrix(NA, nrow = n.mcmc, ncol = ny)
for (t in 1:ny){
  NBratio[,t] <- NB[,t] / Breeders[,t]
}
range(round(apply(NBratio, 2, mean)[2:ny], 2))
mean(apply(NBratio[,2:ny], 1, mean))
quantile(apply(NBratio[,2:ny], 1, mean), c(0.025, 0.975))

# Proportion of immigrants among new breeders 
Iratio <- matrix(NA, nrow = n.mcmc, ncol = ny)
for (t in 1:ny){
  ifelse((mod$sims.list$N[,2,t] + mod$sims.list$N[,4,t] + mod$sims.list$N[,6,t] + mod$sims.list$N[,8,t]) > 0, 
  Iratio[,t] <- mod$sims.list$N[,8,t] / (mod$sims.list$N[,2,t] + mod$sims.list$N[,4,t] + mod$sims.list$N[,6,t] + mod$sims.list$N[,8,t]), Iratio[,t] <- 0) 
}
mean(apply(Iratio[,2:ny], 1, mean))
quantile(apply(Iratio[,2:ny], 1, mean), c(0.025, 0.975))


# By how many time exceed experienced breeders the first time breeders?
fbratio <- rowSums(mod$sims.list$N[,7,2:34]) / rowSums(FB[,2:34])
mean(fbratio)
quantile(fbratio, c(0.025, 0.975))

# By how many time exceed experienced breeders the immigrants?
z <- which(rowSums(mod$sims.list$N[,8,2:34])==0) # which MCMC sample is 0: must be excluded
imratio <- rowSums(mod$sims.list$N[-z,7,2:34]) / rowSums(mod$sims.list$N[-z,8,2:34])
mean(imratio)
quantile(imratio, c(0.025, 0.975))


# Temporal variability of the realized population growth rate
load("Perturbation_Results.Rdata")
p <- c(0.025, 0.975)
m_lam <- exp(log_lam)                 
var_lam <- apply(m_lam, 1, var)
mean(var_lam)
quantile(var_lam, p)     


# Total contribution of temporal variability of rates and population structure
mean(rowSums(cont))      
quantile(rowSums(cont), p)     

# Calculate how much has demography compared to population structure contributed?
# Demographic contributions
mean(rowSums(cont[,1:9]))        
quantile(rowSums(cont[,1:9]), p) 

# Contribution of population structure
mean(rowSums(cont[,10:17]))      
quantile(rowSums(cont[,10:17]), p)     


# Correlation between index of food supply and productivity 
load('IPM_Results.Rdata')
load("NutCond_Results.Rdata")
n.mcmc <- mod$mcmc.info$n.samples

cor_f <- numeric()
for (i in 1:n.mcmc){
  cor_f[i] <- cor(mod$sims.list$f[i,5:34], e1$sims.list$pred[i,1:30])
}

round(mean(cor_f), 3)
round(quantile(cor_f, c(0.025, 0.975)), 3)
round(mean(cor_f>0), 3)



########################################################
#
#  11. Figures for paper
#
########################################################
  
# Figure 1: Population composition

# a. Load data
setwd('...')
load("IPM_Results.Rdata")

# b. Prepare some variables
ny <- dim(mod$mean$N)[2]
year <- 1985:2018
n.mcmc <- mod$mcmc.info$n.samples

# c. Compute total population size and population growth rate
Ntot <- apply(mod$sims.list$N, c(1,3), sum)     
NB <- apply(mod$sims.list$N[,c(1,3,5),], c(1,3), sum)       # locally born non-breeders
Breeders <- apply(mod$sims.list$N[,c(2,4,6,7,8),], c(1,3), sum)
FB <- apply(mod$sims.list$N[,c(2,4,6),], c(1,3), sum)       # local recruits (first time breeders)

# d. Define some figure settings
col.dot <- viridis_pal(option='E')(10)[c(5,8)]
mag <- 3
gs <- 1.62
cex.tif <- mag*1.3
lwd.tif <- 3*mag
d <- 0.25

# e. Choose folder to store figure
setwd('...')

# f. Produce figure
tiff("Fig.1.tif", width = gs*1000*mag, height = 2000*mag, units = "px", pointsize = 35, compression = "lzw")
layout(matrix(1:2, 2, 1, byrow = TRUE), widths = gs, heights = c(1, 1), TRUE)
par(las = 1, cex = cex.tif, mar = c(3, 4, 0.5, 2))

# panel 1: overall
plot(y = apply(Ntot, 2, mean)[2:ny], x = 2:ny, type = "b", pch = 16, ylim = c(0, max(apply(Ntot, 2, quant))), ylab = "Number", axes = F, xlab = NA, xlim = c(1,ny), lwd = lwd.tif)
axis(2, las = 1, lwd = lwd.tif)
axis(1, at = 1:ny, labels = NA, tcl = -0.25, lwd = lwd.tif)
axis(1, at = c(1, 6, 11, 16, 21, 26, 31), labels = NA, tcl = -0.5, lwd = lwd.tif)
axis(1, at = c(1, 6, 11, 16, 21, 26, 31), labels = NA, tcl = -0.5, lwd = lwd.tif)
segments(2:ny, apply(Ntot, 2, quant)[1,2:ny], 2:ny, apply(Ntot, 2, quant)[2,2:ny], lwd = lwd.tif)
points(x = (2:ny)+d, y = apply(Breeders, 2, mean)[2:ny], type = "b", pch = 15, col = col.dot[1], lwd = lwd.tif)
segments((2:ny)+d, apply(Breeders, 2, quant)[1,2:ny], (2:ny)+d, apply(Breeders, 2, quant)[2,2:ny], col = col.dot[1], lwd = lwd.tif)
points(x = (2:ny)-d, y = apply(NB, 2, mean)[2:ny], type = "b", pch = 17, col = col.dot[2], lwd = lwd.tif)
segments((2:ny)-d, apply(NB, 2, quant)[1,2:ny], (2:ny)-d, apply(NB, 2, quant)[2,2:ny], col = col.dot[2], lwd = lwd.tif)
corner.label('A', cex = 1.25)
legend(x = 1, y = max(apply(Ntot, 2, quant))-2, legend = c("Total", "Breeders", "Locally born non-breeders"), pch = c(16,15,17), bty = "n", col = c("black", col.dot))

# panel 2: composition of breeders       
plot(y = mod$mean$N[7,2:ny], x = 2:ny, type = "b", pch = 16, ylim = c(0, max(mod$q97.5$N[7,2:ny])), ylab = "Number", axes = F, xlab = NA, xlim = c(1,ny), lwd = lwd.tif)
axis(2, las = 1, lwd = lwd.tif)
axis(1, at = 1:34, tcl = -0.25, labels = NA, lwd = lwd.tif)
axis(1, at = c(1, 6, 11, 16, 21, 26, 31), labels = c("1985", "1990", "1995", "2000", "2005", "2010", "2015"), lwd = lwd.tif)
segments(2:ny, mod$q2.5$N[7,2:ny], 2:ny, mod$q97.5$N[7,2:ny], lwd = lwd.tif)
points(x = (2:ny)+d, y=apply(FB, 2, mean)[2:ny], type = "b", pch = 15, col = col.dot[1], lwd = lwd.tif)
segments((2:ny)+d, apply(FB, 2, quant)[1,2:ny], (2:ny)+d, apply(FB, 2, quant)[2,2:ny], col = col.dot[1], lwd = lwd.tif)
points(x=(2:ny)-d, y=mod$mean$N[8,2:ny], type = "b", pch = 17, col = col.dot[2], lwd = lwd.tif)
segments((2:ny)-d, mod$q2.5$N[8,2:ny], (2:ny)-d, mod$q97.5$N[8,2:ny], col = col.dot[2], lwd = lwd.tif)
corner.label('B', cex = 1.25)
legend(x = 1, y = max(mod$q97.5$N[7,2:ny])-2, legend = c("Experienced breeders", "Local recruits", "Immigrants"), pch = c(16,15,17), bty = "n", col = c("black", col.dot))

dev.off()


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


# Figure 2. Demographic rates and nestling condition

# a. Load data
setwd('...')
load("IPM_Results.Rdata")
load("NutCond_Results.Rdata")

# b. Prepare (define) some variables
ny <- dim(mod$mean$N)[2]

# c. Define some figure settings
d <- 0.35
co <- brewer.pal(n = 8, name = 'OrRd')[8]
co1 <- brewer.pal(n = 8, name = 'OrRd')[5]
co2 <- col2rgb(co1)
co.t <- rgb(co2[1], co2[2], co2[3], max = 255, alpha = 50)   # Make the colour transparent
mag <- 1.5
cex.tif <- mag*1.4
lwd.tif <- 2.5*mag
lwd.tif2 <- 4*mag
cex.font <- 1.5

# d. Choose folder to store figure
setwd('...')

# e. Produce figure
tiff("Fig.2.tif", width = 3.75*1000*mag, height = 3.1*1000*mag, units = "px", pointsize = 35, compression = "lzw")
layout(matrix(1:9, 3, 3, byrow = TRUE), width = c(1.25,1.25,1.25), height = c(1,1,1.1))

par(cex = cex.tif, mar = c(3, 4, 0.5, 2))
# 1. Survival rates of the three age classes
plot(y = mod$mean$s[1,], x = 1:(ny-1)+0.5, type = "b", pch = 16, xlim = c(1,ny), ylim = c(0.3, 0.9), ylab = expression(paste("Juvenile survival (", s[1], ")")), xlab = "", axes = FALSE, lwd = lwd.tif)
segments(1:(ny-1)+0.5, mod$q2.5$s[1,], 1:(ny-1)+0.5, mod$q97.5$s[1,], lwd = lwd.tif)
axis(1, at = (1:ny), labels = NA, tcl = -0.25, lwd = lwd.tif2)
axis(1, at = c(1, 6, 11, 16, 21, 26, 31), labels = NA, tcl = -0.5, lwd = lwd.tif2)
axis(2, las = 1, lwd = lwd.tif2)
segments(1, mod$mean$mean.s[1], ny+d-0.5, mod$mean$mean.s[1], lwd = lwd.tif2, lty = 2, col = co)
polygon(y = c(mod$q2.5$mean.s[1], mod$q2.5$mean.s[1], mod$q97.5$mean.s[1], mod$q97.5$mean.s[1]), x = c(1,ny+d-0.5,ny+d-0.5,1), border = FALSE, col = co.t)
corner.label('A', cex = 1.25)

plot(y = mod$mean$s[2,], x = 1:(ny-1)+0.5, type = "b", pch = 16, ylim = c(0.3, 0.9), ylab = expression(paste("Subadult survival (", s[2], ")")), xlab = "", axes = FALSE, lwd = lwd.tif)
segments(1:(ny-1)+0.5, mod$q2.5$s[2,], 1:(ny-1)+0.5, mod$q97.5$s[2,], lwd = lwd.tif)
axis(1, at = (1:ny), labels = NA, tcl = -0.25, lwd = lwd.tif2)
axis(1, at = c(1, 6, 11, 16, 21, 26, 31), labels = NA, tcl = -0.5, lwd = lwd.tif2)
axis(2, las = 1, lwd = lwd.tif2)
segments(1, mod$mean$mean.s[2], ny+d-0.5, mod$mean$mean.s[2], lwd = lwd.tif2, lty = 2, col = co)
polygon(y = c(mod$q2.5$mean.s[2], mod$q2.5$mean.s[2], mod$q97.5$mean.s[2], mod$q97.5$mean.s[2]), x = c(1,ny+d-0.5,ny+d-0.5,1), border = FALSE, col = co.t)
corner.label('B', cex = 1.25)

plot(y = mod$mean$s[3,], x = 1:(ny-1)+0.5, type = "b", pch = 16, ylim = c(0.3, 0.9), ylab = expression(paste("Adult survival (", s[3], ")")), xlab = "", axes = FALSE, lwd = lwd.tif)
segments(1:(ny-1)+0.5, mod$q2.5$s[3,], 1:(ny-1)+0.5, mod$q97.5$s[3,], lwd = lwd.tif)
axis(1, at = (1:ny), labels = NA, tcl = -0.25, lwd = lwd.tif2)
axis(1, at = c(1, 6, 11, 16, 21, 26, 31), labels = NA, tcl = -0.5, lwd = lwd.tif2)
axis(2, las = 1, lwd = lwd.tif2)
segments(1, mod$mean$mean.s[3], ny+d-0.5, mod$mean$mean.s[3], lwd = lwd.tif2, lty = 2, col = co)
polygon(y = c(mod$q2.5$mean.s[3], mod$q2.5$mean.s[3], mod$q97.5$mean.s[3], mod$q97.5$mean.s[3]), x = c(1,ny+d-0.5,ny+d-0.5,1), border = FALSE, col = co.t)
corner.label('C', cex = 1.25)

# 2. Productivity
plot(y = mod$mean$f, x = 1:ny, type = "b", pch = 16, ylim = range(c(mod$q2.5$f, mod$q97.5$f)), ylab = "Productivity (f)", xlab = "", axes = FALSE, lwd = lwd.tif)
segments(1:ny, mod$q2.5$f, 1:ny, mod$q97.5$f, lwd = lwd.tif)
axis(1, at = 1:ny, labels = NA, tcl = -0.25, lwd = lwd.tif2)
axis(1, at = c(1, 6, 11, 16, 21, 26, 31), labels = NA, tcl = -0.5, lwd = lwd.tif2)
axis(2, las = 1, lwd = lwd.tif2)
segments(1-d, mod$mean$mean.f, ny+d, mod$mean$mean.f, lwd = lwd.tif2, lty = 2, col = co)
polygon(y = c(mod$q2.5$mean.f, mod$q2.5$mean.f, mod$q97.5$mean.f, mod$q97.5$mean.f), x = c(1-d,ny+d,ny+d,1-d), border = FALSE, col = co.t)
corner.label('D', cex = 1.25)

# 3. Probability to start to reproduce
plot(y = mod$mean$a[1,], x = 2:ny, type = "b", pch = 16, ylim = c(0, 1), ylab = expression(paste("Prob. first breeding (", alpha[1], ")")), xlab = "", axes = FALSE, xlim = c(1,ny), lwd = lwd.tif)
segments(2:ny, mod$q2.5$a[1,], 2:ny, mod$q97.5$a[1,], lwd = lwd.tif)
axis(1, at = 1:ny, labels = NA, tcl = -0.25, lwd = lwd.tif2)
axis(1, at = c(1, 6, 11, 16, 21, 26, 31), labels = NA, tcl = -0.5, lwd = lwd.tif2)
axis(2, las = 1, lwd = lwd.tif2)
segments(2-d, mod$mean$mean.a[1], ny+d, mod$mean$mean.a[1], lwd = lwd.tif2, lty = 2, col = co)
polygon(y = c(mod$q2.5$mean.a[1], mod$q2.5$mean.a[1], mod$q97.5$mean.a[1], mod$q97.5$mean.a[1]), x = c(2-d,ny+d,ny+d,2-d), border = FALSE, col = co.t)
corner.label('E', cex = 1.25)

plot(y = mod$mean$a[2,], x = 2:ny, type = "b", pch = 16, ylim = c(0, 1), ylab = expression(paste("Prob. first breeding (", alpha[2], ")")), xlab = "", axes = FALSE, xlim = c(1,ny), lwd = lwd.tif)
segments(2:ny, mod$q2.5$a[2,], 2:ny, mod$q97.5$a[2,], lwd = lwd.tif)
axis(1, at = 1:ny, labels = NA, tcl = -0.25, lwd = lwd.tif2)
axis(1, at = c(1, 6, 11, 16, 21, 26, 31), labels = NA, tcl = -0.5, lwd = lwd.tif2)
axis(2, las = 1, lwd = lwd.tif2)
segments(2-d, mod$mean$mean.a[2], ny+d, mod$mean$mean.a[2], lwd = lwd.tif2, lty = 2, col = co)
polygon(y = c(mod$q2.5$mean.a[2], mod$q2.5$mean.a[2], mod$q97.5$mean.a[2], mod$q97.5$mean.a[2]), x = c(2-d,ny+d,ny+d,2-d), border = FALSE, col = co.t)
corner.label('F', cex = 1.25)

# 4. Breeding propensity
par(cex = cex.tif, mar = c(4, 4, 0.5, 2))
plot(y = mod$mean$b[2:ny], x = 2:ny, type = "b", pch = 16, ylim = c(0.7, 1), ylab = "Breeding propensity (b)", xlab = "", axes = FALSE)
segments(2:ny, mod$q2.5$b[2:ny], 2:ny, mod$q97.5$b[2:ny], lwd = lwd.tif)
axis(1, at = 1:ny, labels = NA, tcl = -0.25, lwd = lwd.tif2)
axis(1, at = c(1, 6, 11, 16, 21, 26, 31), labels = c("1985", "1990", "1995", "2000", "2005", "2010", "2015"), tcl = -0.5, lwd = lwd.tif2)
axis(2, las = 1, lwd = lwd.tif2)
segments(2-d, mod$mean$mean.b, ny+d, mod$mean$mean.b, lwd = lwd.tif2, lty = 2, col = co)
polygon(y = c(mod$q2.5$mean.b, mod$q2.5$mean.b, mod$q97.5$mean.b, mod$q97.5$mean.b), x = c(2-d,ny+d,ny+d,2-d), border = FALSE, col = co.t)
corner.label('G', cex = 1.25)

# 5. Immigration rate
plot(y = mod$mean$im, x = 1:ny, type = "b", pch = 16, ylim = c(0, 0.5), ylab = expression(paste("Immigration rate (", omega, ")")), xlab = "", axes = FALSE, lwd = lwd.tif)
segments(1:ny, mod$q2.5$im, 1:ny, mod$q97.5$im, lwd = lwd.tif)
axis(1, at = 1:ny, labels = NA, tcl = -0.25, lwd = lwd.tif2)
axis(1, at = c(1, 6, 11, 16, 21, 26, 31), labels = c("1985", "1990", "1995", "2000", "2005", "2010", "2015"), tcl = -0.5, lwd = lwd.tif2)
axis(2, las = 1, lwd = lwd.tif2)
segments(1-d, mod$mean$mean.im, ny+d, mod$mean$mean.im, lwd = lwd.tif2, lty = 2, col = co)
polygon(y = c(mod$q2.5$mean.im, mod$q2.5$mean.im, mod$q97.5$mean.im, mod$q97.5$mean.im), x = c(1-d,ny+d,ny+d,1-d), border = FALSE, col = co.t)
corner.label('H', cex = 1.25)

# 6. Nestling condition
plot(x = 5:34, y = e1$mean$pred, pch = 16, type = "b", ylim = range(c(e1$q2.5$pred, e1$q97.5$pred)), ylab = "Index of nestling condition", xlab = NA, axes = F, xlim = c(1,34), lwd = lwd.tif)
segments(5:34, e1$q2.5$pred, 5:34, e1$q97.5$pred, lwd = lwd.tif)
axis(2, las = 1, lwd = lwd.tif2)
axis(1, at = 1:34, tcl = -0.25, labels = NA, lwd = lwd.tif2)
axis(1, at = c(1, 6, 11, 16, 21, 26, 31), labels = c("1985", "1990", "1995", "2000", "2005", "2010", "2015"), tcl = -0.5, lwd = lwd.tif2)
segments(5-d, e1$mean$alpha, ny+d, e1$mean$alpha, lwd = lwd.tif2, lty = 2, col = co)
polygon(y = c(e1$q2.5$alpha, e1$q2.5$alpha, e1$q97.5$alpha, e1$q97.5$alpha), x = c(5-d,ny+d,ny+d,5-d), border = FALSE, col = co.t)
corner.label('I', cex = 1.25)

dev.off()


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

# Figure 3: tLTRE based on IPM

# a. Load data
setwd('...')
load("Perturbation_Results.Rdata")

# b. Calculate relative contributions
r.cont <- cont/rowSums(cont)

# c. Define some figure settings
names <- c(expression(s[1]), expression(s[2]), expression(s[3]), expression(f), expression(upsilon), expression(alpha[1]), expression(alpha[2]), expression(b), expression(omega), expression(N[1]), expression(N[2]), expression(N[3]), expression(N[4]), expression(N[5]), expression(N[6]), expression(N[7]), expression(N[8]))
col.dot <- c(brewer.pal(n = 8, name = 'OrRd')[7], brewer.pal(n = 8, name = 'Blues')[7]) 
colours <- c(rep(col.dot[1], 9), rep(col.dot[2], 8))
p <- c(0.025, 0.975)
mag <- 2
cex.tif <- mag*1.4
lwd.tif <- 3*mag
cex.font <- 1.5

# d. Choose folder to store figure
setwd('...')

# e. Produce figure
tiff("Fig.3.tif", width = 2*1000*mag, height = 2.1*1000*mag, units = "px", pointsize = 35, compression = "lzw")
layout(matrix(1:2, 2, 1, byrow = TRUE), width = 2, height = c(1,1.1))
par(cex = cex.tif, mar = c(3,4,1,1))

# Panel 1
CRI <- apply(elas, 2, quantile, probs=p)
a <- barplot(colMeans(elas), names.arg = NA, ylab = "Elasticity", las = 1, ylim = range(CRI), col = colours, border = NA, axes = FALSE)
segments(a, CRI[1,], a, CRI[2,], lwd = lwd.tif)
axis(2, las = 1, lwd = lwd.tif)
segments(min(a)-0.5, 0, max(a)+0.5, 0, lwd = lwd.tif)
corner.label('A', cex = 1.25)
legend("topright", pch = c(15, 15), col = col.dot, legend = c("Demographic rates", "Population stage structure"), bty = "n", cex = 1.25)

# Panel 2
par(mar = c(4,4,0,1))
CRI <- apply(r.cont, 2, quantile, probs=p)
a <- barplot(colMeans(r.cont), names.arg = names, ylab = expression(paste("Relative contribution to  ", lambda)), las = 1, ylim = range(CRI), col = colours, lwd = lwd.tif, border = NA)
segments(a, CRI[1,], a, CRI[2,], lwd = lwd.tif)
segments(min(a)-0.5, 0, max(a)+0.5, 0, lwd = lwd.tif)
corner.label('B', cex = 1.25)

dev.off()



############################################
#
# 12. Figures for the appendices
#
############################################


# Appendix 4, Fig. 1: Posterior predictive checks

# a. Load data
setwd('...')
load('IPM_Results.Rdata')

# b. Do some calculations

# c. Define some figure settings
mag <- 4
gs <- 1.62
cex.tif <- mag*1.3
lwd.tif <- 3*mag

# d. Choose folder to store figure
setwd('...')

# e. Make figure
tiff("Fig.A4-1.tif", width = 4*1000*mag, height = 2000*mag, units = "px", pointsize = 35, compression = "lzw")
layout(matrix(1:8, 2, 4, byrow = TRUE), widths = rep(1, 4), heights = c(1, 1), TRUE)
par(las = 1, cex = cex.tif, mar = c(4.1, 4.5, 3, 2))
   
z <- range(c(mod$sims.list$fit.C, mod$sims.list$fitN.C)) 
plot(x = mod$sims.list$fit.C, y = mod$sims.list$fitN.C, pch = 16, ylim = z, xlim = z, ylab = "Replicate data", xlab = "Actual data", las = 1, main = "Count data", col = brewer.pal(n = 8, name = 'OrRd')[7], axes=FALSE, cex=0.3)
axis(1, lwd=lwd.tif)
axis(2, lwd=lwd.tif)
segments(z[1], z[1], z[2], z[2], lwd=lwd.tif/2, lty=2)
text(x=z[1], y = z[2], paste("P = ", round(mean(mod$sims.list$fit.C > mod$sims.list$fitN.C), 3)), adj = c(-0.1,0.5))

z <- range(c(mod$sims.list$fit.dev, mod$sims.list$fitN.dev)) 
plot(x = mod$sims.list$fit.dev, y = mod$sims.list$fitN.dev, pch = 16, ylim = z, xlim = z, ylab = "Replicate data", xlab = "Actual data", las = 1, main = "Productivity data", col = brewer.pal(n = 8, name = 'OrRd')[7], axes=FALSE, cex=0.3)
axis(1, lwd=lwd.tif)
axis(2, lwd=lwd.tif)
segments(z[1], z[1], z[2], z[2], lwd=lwd.tif/2, lty=2)
text(x=z[1], y = z[2], paste("P = ", round(mean(mod$sims.list$fit.dev > mod$sims.list$fitN.dev), 3)), adj = c(-0.1,0.5))

z <- range(c(mod$sims.list$fit.1, mod$sims.list$fitN.1)) 
plot(x = mod$sims.list$fit.1, y = mod$sims.list$fitN.1, pch = 16, ylim = z, xlim = z, ylab = "Replicate data", xlab = "Actual data", las = 1, main = "Resightings, juv. (wing tags)", col = brewer.pal(n = 8, name = 'OrRd')[7], axes=FALSE, cex=0.3)
axis(1, lwd=lwd.tif)
axis(2, lwd=lwd.tif)
segments(z[1], z[1], z[2], z[2], lwd=lwd.tif/2, lty=2)
text(x=z[1], y = z[2], paste("P = ", round(mean(mod$sims.list$fit.1 > mod$sims.list$fitN.1), 3)), adj = c(-0.1,0.5))

z <- range(c(mod$sims.list$fit.4, mod$sims.list$fitN.4)) 
plot(x = mod$sims.list$fit.4, y = mod$sims.list$fitN.4, pch = 16, ylim = z, xlim = z, ylab = "Replicate data", xlab = "Actual data", las = 1, main = "Resightings, juv. (transmitter)", col = brewer.pal(n = 8, name = 'OrRd')[7], axes=FALSE, cex=0.3)
axis(1, lwd=lwd.tif)
axis(2, lwd=lwd.tif)
segments(z[1], z[1], z[2], z[2], lwd=lwd.tif/2, lty=2)
text(x=z[1], y = z[2], paste("P = ", round(mean(mod$sims.list$fit.4 > mod$sims.list$fitN.4), 3)), adj = c(-0.1,0.5))

z <- range(c(mod$sims.list$fit.2, mod$sims.list$fitN.2)) 
plot(x = mod$sims.list$fit.2, y = mod$sims.list$fitN.2, pch = 16, ylim = z, xlim = z, ylab = "Replicate data", xlab = "Actual data", las = 1, main = "Resightings, ad. (wing tags)", col = brewer.pal(n = 8, name = 'OrRd')[7], axes=FALSE, cex=0.3)
axis(1, lwd=lwd.tif)
axis(2, lwd=lwd.tif)
segments(z[1], z[1], z[2], z[2], lwd=lwd.tif/2, lty=2)
text(x=z[1], y = z[2], paste("P = ", round(mean(mod$sims.list$fit.2 > mod$sims.list$fitN.2), 3)), adj = c(-0.1,0.5))


z <- range(c(mod$sims.list$fit.3, mod$sims.list$fitN.3)) 
plot(x = mod$sims.list$fit.3, y = mod$sims.list$fitN.3, pch = 16, ylim = z, xlim = z, ylab = "Replicate data", xlab = "Actual data", las = 1, main = "Resightings, ad. (transmitter)", col = brewer.pal(n = 8, name = 'OrRd')[7], axes=FALSE, cex=0.3)
axis(1, lwd=lwd.tif)
axis(2, lwd=lwd.tif)
segments(z[1], z[1], z[2], z[2], lwd=lwd.tif/2, lty=2)
text(x=z[1], y = z[2], paste("P = ", round(mean(mod$sims.list$fit.3 > mod$sims.list$fitN.3), 3)), adj = c(-0.1,0.5))

z <- range(c(mod$sims.list$fit.5, mod$sims.list$fitN.5)) 
plot(x = mod$sims.list$fit.5, y = mod$sims.list$fitN.5, pch = 16, ylim = z, xlim = z, ylab = "Replicate data", xlab = "Actual data", las = 1, main = "Recovery data Hiddensee", col = brewer.pal(n = 8, name = 'OrRd')[7], axes=FALSE, cex=0.3)
axis(1, lwd=lwd.tif)
axis(2, lwd=lwd.tif)
segments(z[1], z[1], z[2], z[2], lwd=lwd.tif/2, lty=2)
text(x=z[1], y = z[2], paste("P = ", round(mean(mod$sims.list$fit.5 > mod$sims.list$fitN.5), 3)), adj = c(-0.1,0.5))

dev.off()

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

# Appendix 4, Fig. 2: Histogram of the R-hat values

# a. Load data
setwd('...')
load("IPM_Results.Rdata")

# b. Do some calculations

# c. Define some figure settings
mag <- 4
cex.tif <- mag*1
lwd.tif <- 3*mag

# d. Choose folder to store figure
setwd('...')

# e. Make figure
tiff("Fig.A4-2.tif", width = 1000*mag, height = 1000*mag, units = "px", pointsize = 35, compression = "lzw")
par(las = 1, cex = cex.tif, mar = c(4.1, 4.5, 3, 2))

# Boxplot of the structural parameters (means and SD of demographic rates and of population sizes)
boxplot(mod$summary[c(1:65, 512:783),8], ylab = "R-hat", lwd=lwd.tif, axes=FALSE, cex=0.5)
axis(2, lwd=lwd.tif)

dev.off()



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

# Appendix 5, Fig. 1: visualisation of ring recoveries

# a. Load data
setwd('...')
ReReData <- read.csv('RingRecoveries.csv', header = TRUE, sep = ';', dec = ",")
load("MapsData.Rdata")

# b. Information about the MapsData.Rdata
# The data are from this webpage: https://gadm.org/
# de0 <- getData('GADM', country='DEU', level=0)
# de1 <- getData('GADM', country='DEU', level=1)
# save(de0, de1, file = "MapsData.Rdata")
# Koordinaten: https://www.koordinaten-umrechner.de/Dezimal/51.062867,11.361269?karte=OpenStreetMap&zoom=8

# c. Do some calculations
X <- ReReData[,39]
Y <- ReReData[,40]
W <- ReReData[,38]

# d. Figure settings
mag <- 5
cex.tif <- mag/1.5
lwd.tif <- 2*mag

# e. Choose folder to store figure
setwd('...')

# f. Make figure
tiff("Fig.A5-1.tif", width = 1000*mag, height = 1000*mag, units = "px", pointsize = 35, compression = "lzw")
par(mai=c(0,0,0,0), cex=cex.tif)

c1 <- brewer.pal(n = 8, name = 'OrRd')[2]
co <- c(NA, NA, c1, c1, NA, NA, NA, c1, NA, NA, NA, NA, c1, c1, NA, c1)
plot(de1, col = co, axes = F, las = 1, border = co)
plot(de0, lwd = lwd.tif, add=TRUE)
points(X[W=="nein"], Y[W=="nein"], pch = 16, cex = 0.5, col = brewer.pal(n = 8, name = 'OrRd')[7])
points(X[W=="ja"], Y[W=="ja"], pch = 16, col = brewer.pal(n = 8, name = 'Blues')[8], cex = 0.5)

scale.coord <- c(5, 54.5)
d <- 2.36
segments(scale.coord[1], scale.coord[2], scale.coord[1]+d, scale.coord[2], lwd=lwd.tif)
segments(scale.coord[1], scale.coord[2], scale.coord[1], scale.coord[2]-0.1, lwd=lwd.tif)
segments(scale.coord[1]+d/3, scale.coord[2], scale.coord[1]+d/3, scale.coord[2]-0.1, lwd=lwd.tif)
segments(scale.coord[1]+2*d/3, scale.coord[2], scale.coord[1]+2*d/3, scale.coord[2]-0.1, lwd=lwd.tif)
segments(scale.coord[1]+d, scale.coord[2], scale.coord[1]+d, scale.coord[2]-0.1, lwd=lwd.tif)
text(x = scale.coord[1], y = scale.coord[2] - 0.08, labels = "0", pos = 1)
text(x = scale.coord[1] + 0.8, y = scale.coord[2] - 0.08, labels = "50", pos = 1)
text(x = scale.coord[1] + 1.55, y = scale.coord[2] - 0.08, labels = "100", pos = 1)
text(x = scale.coord[1] + 2.35, y = scale.coord[2] - 0.08, labels = "150", pos = 1)

dev.off()

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


# Appendix 5, Fig. 2: Ring recovery probabilities

# a. Load data
setwd('...')
load('IPM_Results.Rdata')

# b. Do some calculations
ny <- dim(mod$mean$N)[2]

# c. Figure settings
co <- brewer.pal(n = 8, name = 'OrRd')[8]
co1 <- brewer.pal(n = 8, name = 'OrRd')[5]
co2 <- col2rgb(co1)
co.t <- rgb(co2[1], co2[2], co2[3], max = 255, alpha = 50)   # Make the colour transparent
mag <- 4
cex.tif <- mag*1.4
lwd.tif <- 3*mag
cex.font <- 1.5
gs <- 1.6

# d. Choose folder to store figure
setwd('...')

# e. Make figure
tiff("Fig.A5-2.tif", width = 2*1000*gs*mag, height = 2.1*1000*mag, units = "px", pointsize = 35, compression = "lzw")
layout(matrix(1:4, 2, 2, byrow = TRUE), width = c(gs, gs), height = c(1,1.1))
par(cex = cex.tif, mar = c(2.5,4.2,1,1), las=1)

plot(y = mod$mean$r[1,], x = 1:(ny-1)+0.5, type = "b", pch = 16, xlim = c(1,ny), ylim = c(0, 0.35), ylab = expression(paste("Juvenile dead recovery (ring, ", r[1], ")")), xlab = "", axes = FALSE, lwd=lwd.tif)
segments(1:(ny-1)+0.5, mod$q2.5$r[1,], 1:(ny-1)+0.5, mod$q97.5$r[1,], lwd=lwd.tif)
axis(1, at = (1:ny), labels = NA, tcl = -0.25, lwd=lwd.tif)
axis(1, at = c(1, 6, 11, 16, 21, 26, 31), labels = NA, tcl = -0.5, lwd=lwd.tif)
axis(2, las = 1, lwd=lwd.tif)
segments(1, mod$mean$mean.r[1], ny, mod$mean$mean.r[1], lty = 2, col = co, lwd=lwd.tif)
polygon(y = c(mod$q2.5$mean.r[1], mod$q2.5$mean.r[1], mod$q97.5$mean.r[1], mod$q97.5$mean.r[1]), x = c(1,ny,ny,1), border = FALSE, col = co.t)

plot(y = mod$mean$r[2,], x = 1:(ny-1)+0.5, type = "b", pch = 16, xlim = c(1,ny), ylim = c(0, 0.35), ylab = expression(paste("Adult dead recovery (ring, ", r[a], ")")), xlab = "", axes = FALSE, lwd=lwd.tif)
segments(1:(ny-1)+0.5, mod$q2.5$r[2,], 1:(ny-1)+0.5, mod$q97.5$r[2,], lwd=lwd.tif)
axis(1, at = (1:ny), labels = NA, tcl = -0.25, lwd=lwd.tif)
axis(1, at = c(1, 6, 11, 16, 21, 26, 31), labels = NA, tcl = -0.5, lwd=lwd.tif)
axis(2, las = 1, lwd=lwd.tif)
segments(1, mod$mean$mean.r[2], ny, mod$mean$mean.r[2], lwd=lwd.tif, lty = 2, col = co)
polygon(y = c(mod$q2.5$mean.r[2], mod$q2.5$mean.r[2], mod$q97.5$mean.r[2], mod$q97.5$mean.r[2]), x = c(1,ny,ny,1), border = FALSE, col = co.t)

# Note: nestlings with wing tags were released in years 1:4, 9, 11, 12, 14 until the end. Therefore there are no annual estimates in some years.
plot(y = mod$mean$r[3,c(1:4,9,11,12,14:(ny-1))], x = c(1:4,9,11,12,14:(ny-1))+0.5, type = "b", pch = 16, xlim = c(1,ny), ylim = c(0, 0.35), ylab = expression(paste("Juvenile dead recovery (wing tag, ", r^w, "" [1], ")")), xlab = "", axes = FALSE, lwd=lwd.tif)
segments(c(1:4,9,11,12,14:(ny-1))+0.5, mod$q2.5$r[3,c(1:4,9,11,12,14:(ny-1))], c(1:4,9,11,12,14:(ny-1))+0.5, mod$q97.5$r[3,c(1:4,9,11,12,14:(ny-1))], lwd=lwd.tif)
axis(1, at = (1:ny), labels = NA, tcl = -0.25, lwd=lwd.tif)
axis(1, at = c(1, 6, 11, 16, 21, 26, 31), labels = c("1985", "1990", "1995", "2000", "2005", "2010", "2015"), tcl = -0.5, lwd=lwd.tif)
axis(2, las = 1, lwd=lwd.tif)
segments(1, mod$mean$mean.r[3], ny, mod$mean$mean.r[3], lwd=lwd.tif, lty = 2, col = co)
polygon(y = c(mod$q2.5$mean.r[3], mod$q2.5$mean.r[3], mod$q97.5$mean.r[3], mod$q97.5$mean.r[3]), x = c(1,ny,ny,1), border = FALSE, col = co.t)

Note: no adult with a wing tag was released in year 1, hence there is no annual estimate in that year.
plot(y = mod$mean$r[4,2:(ny-1)], x = 2:(ny-1)+0.5, type = "b", pch = 16, xlim = c(1,ny), ylim = c(0, 0.35), ylab = expression(paste("Adult dead recovery (wing tag, ", r^w, "" [a], ")")), xlab = "", axes = FALSE, lwd=lwd.tif)
segments(2:(ny-1)+0.5, mod$q2.5$r[4,2:(ny-1)], 2:(ny-1)+0.5, mod$q97.5$r[4,2:(ny-1)], lwd=lwd.tif)
axis(1, at = (1:ny), labels = NA, tcl = -0.25, lwd=lwd.tif)
axis(1, at = c(1, 6, 11, 16, 21, 26, 31), labels = c("1985", "1990", "1995", "2000", "2005", "2010", "2015"), tcl = -0.5, lwd=lwd.tif)
axis(2, las = 1, lwd=lwd.tif)
segments(2, mod$mean$mean.r[4], ny, mod$mean$mean.r[4], lwd=lwd.tif, lty = 2, col = co)
polygon(y = c(mod$q2.5$mean.r[4], mod$q2.5$mean.r[4], mod$q97.5$mean.r[4], mod$q97.5$mean.r[4]), x = c(2,ny,ny,2), border = FALSE, col = co.t)

dev.off()