###############################################################
#
# 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 results provided in the text, for the prospective (elasticities) and retrospective (tLTRE) analyses, and for testing density-dependence at the population level
#
###############################################################


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

library(IPMbook)
library(nimble)
library(MCMCvis)
library(scales)
library(RColorBrewer)


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

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


#########################################
#
# Load data: population counts
#
#########################################

pop <- read.csv('Population-count_data.csv', header=TRUE, sep=';')
count <- pop[,2]


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

source('Functions.txt')


#########################################
#
# Load result file
#
#########################################

load('Model1.Rdata')

res_ipm <- res_ipm1
out <- do.call(rbind, res_ipm$samples)


# Find indices for parameters
parms <- colnames(out)
nyears <- 47

i_N0F1 <- seq(from=1, by=3, length.out=nyears)          # female fledglings produced by 1y females
i_N0F2 <- seq(from=2, by=3, length.out=nyears)          # female fledglings produced by 2y females
i_N0F3 <- seq(from=3, by=3, length.out=nyears)          # female fledglings produced by ad females
i_N1F <- seq(from=606, by=10, length.out=nyears)        # N[1, females]
i_N2F <- seq(from=607, by=10, length.out=nyears)        # N[2, females]
i_N3F <- seq(from=608, by=10, length.out=nyears)        # N[3, females]
i_N4F <- seq(from=609, by=10, length.out=nyears)        # N[4, females]
i_N5F <- seq(from=610, by=10, length.out=nyears)        # N[5, females]
i_NFF1 <- seq(from=142, by=6, length.out=(nyears-1))    # female locally born 1y first-time breeder
i_NFF2 <- seq(from=143, by=6, length.out=(nyears-1))    # female locally born 2y first-time breeder
i_NFF3 <- seq(from=144, by=6, length.out=(nyears-1))    # female locally born 3y first-time breeder

i_N0M1 <- seq(from=465, by=3, length.out=nyears)        # male fledglings produced by 1y females
i_N0M2 <- seq(from=466, by=3, length.out=nyears)        # male fledglings produced by 2y females
i_N0M3 <- seq(from=467, by=3, length.out=nyears)        # male fledglings produced by ad females
i_N1M <- seq(from=611, by=10, length.out=nyears)        # N[1, males]
i_N2M <- seq(from=612, by=10, length.out=nyears)        # N[2, males]
i_N3M <- seq(from=613, by=10, length.out=nyears)        # N[3, males]
i_N4M <- seq(from=614, by=10, length.out=nyears)        # N[4, males]
i_N5M <- seq(from=615, by=10, length.out=nyears)        # N[5, males]
i_NI <- 418:464                                         # Number of male immigrants
i_NFM1 <- seq(from=145, by=6, length.out=(nyears-1))    # male locally born 1y first-time breeder
i_NFM2 <- seq(from=146, by=6, length.out=(nyears-1))    # male locally born 2y first-time breeder
i_NFM3 <- seq(from=147, by=6, length.out=(nyears-1))    # male locally born 3y first-time breeder

i_phiF <- seq(from=1608, by=2, length.out=(nyears-1))   # female adult survival
i_phiM <- seq(from=1609, by=2, length.out=(nyears-1))   # male adult survival
i_p1F <- seq(from=1424, by=4, length.out=(nyears-1))    # p (seen before, females)
i_p2F <- seq(from=1425, by=4, length.out=(nyears-1))    # p (not seen before, females)
i_p1M <- seq(from=1426, by=4, length.out=(nyears-1))    # p (seen before, males)
i_p2M <- seq(from=1427, by=4, length.out=(nyears-1))    # p (not seen before, males)
i_mean.p <- 1409:1410

i_s1F <- seq(from=1933, by=6, length.out=nyears-1)
i_s2F <- seq(from=1934, by=6, length.out=nyears-1)
i_s3F <- seq(from=1935, by=6, length.out=nyears-1)
i_s1M <- seq(from=1936, by=6, length.out=nyears-1)
i_s2M <- seq(from=1937, by=6, length.out=nyears-1)
i_s3M <- seq(from=1938, by=6, length.out=nyears-1)
i_mean.SF <- 1415:1417
i_mean.SM <- 1418:1420
i_sigma.SF <- 2225:2227
i_sigma.SM <- 2228:2230

i_rF <- seq(from=1700, by=2, length.out=nyears-1)
i_rM <- seq(from=1701, by=2, length.out=nyears-1)
i_mean.r <- 1411

i_eta1 <- seq(from=1260, by=3, length.out=nyears)
i_eta2 <- seq(from=1261, by=3, length.out=nyears)
i_eta3 <- seq(from=1262, by=3, length.out=nyears)
i_mean.eta <- 1405:1407
i_sigma.eta <- 2214:2216

i_rho1 <- seq(from=1792, by=3, length.out=nyears)
i_rho2 <- seq(from=1793, by=3, length.out=nyears)
i_rho3 <- seq(from=1794, by=3, length.out=nyears)
i_mean.rho <- 1412:1414
i_sigma.rho <- 2222:2224

i_xi1 <- seq(from=2234, by=3, length.out=nyears)
i_xi2 <- seq(from=2235, by=3, length.out=nyears)
i_xi3 <- seq(from=2236, by=3, length.out=nyears)
i_mean.xi <- 1421:1423
i_sigma.xi <- 2231:2233

i_alpha1F <- seq(from=1076, by=4, length.out=nyears-1)
i_alpha2F <- seq(from=1077, by=4, length.out=nyears-1)
i_alpha1M <- seq(from=1078, by=4, length.out=nyears-1)
i_alpha2M <- seq(from=1079, by=4, length.out=nyears-1)
i_mean.alphaF <- 1401:1402
i_mean.alphaM <- 1403:1404
i_sigma.alphaF <- 2209:2210
i_sigma.alphaM <- 2211:2212

i_mean.omega <- 1408


#########################
#
# Misc stats that is used in the text
#
#########################

# 1. Identify the year with the highest and lowest number of breeding pairs
year <- 1968:2014
BreederF <- out[,i_N2F] + out[,i_N4F] + out[,i_N5F]
BreederM <- out[,i_N2M] + out[,i_N4M] + out[,i_N5M] + out[,i_NI] 
BR <- BreederF + BreederM
lo <- which.min(apply(BR, 2, mean))
hi <- which.max(apply(BR, 2, mean))
year[c(lo, hi)]
round(apply(BR, 2, mean)[c(lo, hi)]/2)
round(apply(BR/2, 2, quantile, c(0.025, 0.975))[,c(lo, hi)])


# 2. CV of the number of estimated breeding pairs
sd(apply(BR, 2, mean)/2)/mean(apply(BR, 2, mean)/2)*100


# 3. Geometric population growth rate
BreederF <- out[,i_N2F] + out[,i_N4F] + out[,i_N5F]
FloaterF <- out[,i_N1F] + out[,i_N3F]
TotalF <- BreederF + FloaterF
BreederM <- out[,i_N2M] + out[,i_N4M] + out[,i_N5M] + out[,i_NI] 
FloaterM <- out[,i_N1M] + out[,i_N3M]
TotalM <- BreederM + FloaterM
Ntot <- TotalF + TotalM

# Compute the realized population growth rate
log_lam <- log(Ntot[,-1]) - log(Ntot[,-47])    # Realized growth rates (log)

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

# Decline over the 46 years
round(mean(1-exp(apply(log_lam, 1, sum))), 3)
round(quantile(1-exp(apply(log_lam, 1, sum)), c(0.025, 0.975)), 3)


# 4. Probability that female survival is higher than male survival
apply(out[,i_mean.SF] > out[,i_mean.SM], 2, mean)


# 5. Probability that the recruitment probability of females is higher at age=2 than at age=1
mean(out[,i_mean.alphaF[1]] < out[,i_mean.alphaF[2]])


# 6. Probability that the recruitment probability of males is higher at age=2 than at age=1
mean(out[,i_mean.alphaM[1]] < out[,i_mean.alphaM[2]])


# 7. Probability that the recruitment probability of females at age 1 is higher than that of males
mean(out[,i_mean.alphaF[1]] > out[,i_mean.alphaM[1]])


# 8. Productivity per initiated brood
# Calculate the output per female and year
op <- out[,i_mean.eta] * out[,i_mean.rho]
round(apply(op, 2, mean), 2)
round(apply(op, 2, quantile, c(0.025, 0.975)), 2)


# 8. Fledgling sex ratio
mean(out[,i_mean.xi[1]] < out[,i_mean.xi[2]])
mean(out[,i_mean.xi[1]] < out[,i_mean.xi[3]])
mean(out[,i_mean.xi[2]] < out[,i_mean.xi[3]])

# Mean proportion of female fledglings
FledglingF <- out[,i_N0F1] + out[,i_N0F2] + out[,i_N0F3]
FledglingM <- out[,i_N0M1] + out[,i_N0M2] + out[,i_N0M3]
round(mean(FledglingF / (FledglingF + FledglingM)), 2)
round(quantile(FledglingF / (FledglingF + FledglingM), c(0.025, 0.975)), 2)
range(round(apply(FledglingF / (FledglingF + FledglingM), 2, mean), 2))


# 9. Adult sex ratio
Nf <- out[,i_N1F] + out[,i_N2F] + out[,i_N3F] + out[,i_N4F] + out[,i_N5F]
Nm <- out[,i_N1M] + out[,i_N2M] + out[,i_N3M] + out[,i_N4M] + out[,i_N5M] + out[,i_NI]
ASR <- Nm / (Nm + Nf)
round(mean(ASR), 3)
quantile(ASR, c(0.025, 0.975))
mean(ASR<0.5)


# 10. Average proportion of fledglings, non-breeders and breeders
# both sexes combined
FledglingF <- out[,i_N0F1] + out[,i_N0F2] + out[,i_N0F3]
BreederF <- out[,i_N2F] + out[,i_N4F] + out[,i_N5F]
FloaterF <- out[,i_N1F] + out[,i_N3F]
TotalF <- BreederF + FloaterF + FledglingF
FledglingM <- out[,i_N0M1] + out[,i_N0M2] + out[,i_N0M3]
BreederM <- out[,i_N2M] + out[,i_N4M] + out[,i_N5M] + out[,i_NI]     
FloaterM <- out[,i_N1M] + out[,i_N3M]
TotalM <- BreederM + FloaterM + FledglingM
Total <- TotalF + TotalM

round(mean((FledglingF + FledglingM) / Total), 2)
round(quantile((FledglingF + FledglingM) / Total, c(0.025, 0.975)), 2)
round(mean((FloaterF + FloaterM)/ Total), 2)
round(quantile((FloaterF + FloaterM)/ Total, c(0.025, 0.975)), 2)
round(mean((BreederF + BreederM) / Total), 2)
round(quantile((BreederF + BreederM) / Total, c(0.025, 0.975)), 2)


# 11. Average proportion of experienced and first time breeders among breeders
# Females
FirstF <- out[,i_NFF1] + out[,i_NFF2] + out[,i_NFF3]
ExpBF <- BreederF[,2:47] - FirstF
round(mean(ExpBF / BreederF[,2:47]), 2)
round(quantile(ExpBF / BreederF[,2:47], c(0.025, 0.975)), 2)
round(mean(FirstF / BreederF[,2:47]), 2)
round(quantile(FirstF / BreederF[,2:47], c(0.025, 0.975)), 2)

# Males
FirstM <- out[,i_NFM1] + out[,i_NFM2] + out[,i_NFM3]
ImM <- out[,i_NI[2:47]]
ExpBM <- BreederM[,2:47] - FirstM - ImM
round(mean(ExpBM / BreederM[,2:47]), 2)
round(quantile(ExpBM / BreederM[,2:47], c(0.025, 0.975)), 2)
round(mean(FirstM / BreederM[,2:47]), 2)
round(quantile(FirstM / BreederM[,2:47], c(0.025, 0.975)), 2)
round(mean(ImM / BreederM[,2:47]), 2)
round(quantile(ImM / BreederM[,2:47], c(0.025, 0.975)), 2)


# 12. Stength of density dependence on female juvenile survival as obtained from model 3
load('Model3.Rdata')

res_ipmDD <- res_ipm3
outDD <- do.call(rbind, res_ipmDD$samples)
round(mean(outDD[,'beta.s']), 3)
round(quantile(outDD[,'beta.s'], c(0.025, 0.975)), 3)


# 13. Number of female immigrants
# Load model with female immigration
load('Model4.Rdata')
res_ipmI <- res_ipm4
outI <- do.call(rbind, res_ipmI$samples)
median(outI[,'mean.omega[1]'])



#########################################################
#
# Values for Table A-2
#
#########################################################

load('Model1.Rdata')
res_ipm <- res_ipm1
out <- do.call(rbind, res_ipm$samples)

ro <- 3   # digits for rounding

# Means
# Female survival
round(apply(out[,i_mean.SF], 2, mean), ro)
round(apply(out[,i_mean.SF], 2, quantile, c(0.025, 0.975)), ro)

# Male survival
round(apply(out[,i_mean.SM], 2, mean), ro)
round(apply(out[,i_mean.SM], 2, quantile, c(0.025, 0.975)), ro)

# Breeding success
round(apply(out[,i_mean.eta], 2, mean), ro)
round(apply(out[,i_mean.eta], 2, quantile, c(0.025, 0.975)), ro)

# Productivity
round(apply(out[,i_mean.rho], 2, mean), ro)
round(apply(out[,i_mean.rho], 2, quantile, c(0.025, 0.975)), ro)

# Fledgling sex ratio
round(apply(out[,i_mean.xi], 2, mean), ro)
round(apply(out[,i_mean.xi], 2, quantile, c(0.025, 0.975)), ro)

# Recruitment females
round(apply(out[,i_mean.alphaF], 2, mean), ro)
round(apply(out[,i_mean.alphaF], 2, quantile, c(0.025, 0.975)), ro)

# Recruitment males
round(apply(out[,i_mean.alphaM], 2, mean), ro)
round(apply(out[,i_mean.alphaM], 2, quantile, c(0.025, 0.975)), ro)


# Temporal variability
# Female survival
round(apply(out[,i_sigma.SF], 2, mean), ro)
round(apply(out[,i_sigma.SF], 2, quantile, c(0.025, 0.975)), ro)

# Male survival
round(apply(out[,i_sigma.SM], 2, mean), ro)
round(apply(out[,i_sigma.SM], 2, quantile, c(0.025, 0.975)), ro)

# Breeding success
round(apply(out[,i_sigma.eta], 2, mean), ro)
round(apply(out[,i_sigma.eta], 2, quantile, c(0.025, 0.975)), ro)

# Productivity
round(apply(out[,i_sigma.rho], 2, mean), ro)
round(apply(out[,i_sigma.rho], 2, quantile, c(0.025, 0.975)), ro)

# Fledgling sex ratio
round(apply(out[,i_sigma.xi], 2, mean), ro)
round(apply(out[,i_sigma.xi], 2, quantile, c(0.025, 0.975)), ro)

# Recruitment females
round(apply(out[,i_sigma.alphaF], 2, mean), ro)
round(apply(out[,i_sigma.alphaF], 2, quantile, c(0.025, 0.975)), ro)

# Recruitment males
round(apply(out[,i_sigma.alphaM], 2, mean), ro)
round(apply(out[,i_sigma.alphaM], 2, quantile, c(0.025, 0.975)), ro)




#########################################################
#
# tLTRE: decomposition of the temporal variances
#
#########################################################

# The LTRE is performed for model 1

load('Model1.Rdata')
res_ipm <- res_ipm1
out <- do.call(rbind, res_ipm$samples)

year <- 1968:2014
nyears <- length(year)
n.draws <- nrow(out)

BreederF <- out[,i_N2F] + out[,i_N4F] + out[,i_N5F]
FloaterF <- out[,i_N1F] + out[,i_N3F]
TotalF <- BreederF + FloaterF

BreederM <- out[,i_N2M] + out[,i_N4M] + out[,i_N5M] + out[,i_NI] 
immRate <- out[,i_NI[-2]] / (BreederM[,-47])
FloaterM <- out[,i_N1M] + out[,i_N3M]
TotalM <- BreederM + FloaterM

Ntot <- TotalF + TotalM

# Compute the realized population growth rate
log_lam <- log(Ntot[,-1]) - log(Ntot[,-47])    # Realized growth rates (log)

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


# Define lambda
lambda <- expression((sf1 * (n2*xi1*eta1*rho1 + n4*xi2*eta2*rho2 + n5*xi3*eta3*rho3) + sm1 * (n2*(1-xi1)*eta1*rho1 + n4*(1-xi2)*eta2*rho2 + n5*(1-xi3)*eta3*rho3) + (n1+n2)*sf2 + (n3+n4+n5)*sf3 + (n6+n7)*sm2 + (n8+n9+n10+n11)*(sm3 + im)) / (n1 + n2 + n3 + n4 + n5 + n6 + n7 + n8 + n9 + n10 + n11))   

# --> Note: alpha is not part of this equation, because the transition to first breeding is captured in the stage structure

# Calculate proportional population sizes
n1 <- out[,i_N1F] / Ntot
n2 <- out[,i_N2F] / Ntot
n3 <- out[,i_N3F] / Ntot
n4 <- out[,i_N4F] / Ntot
n5 <- out[,i_N5F] / Ntot
n6 <- out[,i_N1M] / Ntot
n7 <- out[,i_N2M] / Ntot
n8 <- out[,i_N3M] / Ntot
n9 <- out[,i_N4M] / Ntot
n10 <- out[,i_N5M] / Ntot
n11 <- out[,i_NI] / Ntot


# Extract the mean demographic rates and population sizes and store them in a list
mu <- list(sf1=out[,i_mean.SF[1]], sf2=out[,i_mean.SF[2]], sf3=out[,i_mean.SF[3]], sm1=out[,i_mean.SM[1]], sm2=out[,i_mean.SM[2]], sm3=out[,i_mean.SM[3]], eta1=out[,i_mean.eta[1]], eta2=out[,i_mean.eta[2]], eta3=out[,i_mean.eta[3]], rho1=out[,i_mean.rho[1]], rho2=out[,i_mean.rho[2]], rho3=out[,i_mean.rho[3]], xi1=out[,i_mean.xi[1]], xi2=out[,i_mean.xi[2]], xi3=out[,i_mean.xi[3]], im = apply(immRate, 1, mean), n1=apply(n1, 1, mean, na.rm=TRUE), n2=apply(n2, 1, mean, na.rm=TRUE), n3=apply(n3, 1, mean, na.rm=TRUE), n4=apply(n4, 1, mean, na.rm=TRUE), n5=apply(n5, 1, mean, na.rm=TRUE), n6=apply(n6, 1, mean, na.rm=TRUE), n7=apply(n7, 1, mean, na.rm=TRUE), n8=apply(n8, 1, mean, na.rm=TRUE), n9=apply(n9, 1, mean, na.rm=TRUE), n10=apply(n10, 1, mean, na.rm=TRUE), n11=apply(n11, 1, mean, na.rm=TRUE))

# Calculate sensitivities
sens <- matrix(NA, n.draws, 27)  
colnames(sens) <- c("sf1", "sf2", "sf3", "sm1", "sm2", "sm3", "eta1", "eta2", "eta3", "rho1", "rho2", "rho3", "xi1", "xi2", "xi3", "im", "n1", "n2", "n3", "n4", "n5", "n6", "n7", "n8", "n9", "n10", "n11")
sens[,"sf1"] <- eval(D(lambda, "sf1"), envir=mu)
sens[,"sf2"] <- eval(D(lambda, "sf2"), envir=mu)
sens[,"sf3"] <- eval(D(lambda, "sf3"), envir=mu)
sens[,"sm1"] <- eval(D(lambda, "sm1"), envir=mu)
sens[,"sm2"] <- eval(D(lambda, "sm2"), envir=mu)
sens[,"sm3"] <- eval(D(lambda, "sm3"), envir=mu)
sens[,"eta1"] <- eval(D(lambda, "eta1"), envir=mu)
sens[,"eta2"] <- eval(D(lambda, "eta2"), envir=mu)
sens[,"eta3"] <- eval(D(lambda, "eta3"), envir=mu)
sens[,"rho1"] <- eval(D(lambda, "rho1"), envir=mu)
sens[,"rho2"] <- eval(D(lambda, "rho2"), envir=mu)
sens[,"rho3"] <- eval(D(lambda, "rho3"), envir=mu)
sens[,"xi1"] <- eval(D(lambda, "xi1"), envir=mu)
sens[,"xi2"] <- eval(D(lambda, "xi2"), envir=mu)
sens[,"xi3"] <- eval(D(lambda, "xi3"), envir=mu)
sens[,"im"] <- eval(D(lambda, "im"), envir=mu)
sens[,"n1"] <- eval(D(lambda, "n1"), envir=mu)
sens[,"n2"] <- eval(D(lambda, "n2"), envir=mu)
sens[,"n3"] <- eval(D(lambda, "n3"), envir=mu)
sens[,"n4"] <- eval(D(lambda, "n4"), envir=mu)
sens[,"n5"] <- eval(D(lambda, "n5"), envir=mu)
sens[,"n6"] <- eval(D(lambda, "n6"), envir=mu)
sens[,"n7"] <- eval(D(lambda, "n7"), envir=mu)
sens[,"n8"] <- eval(D(lambda, "n8"), envir=mu)
sens[,"n9"] <- eval(D(lambda, "n9"), envir=mu)
sens[,"n10"] <- eval(D(lambda, "n10"), envir=mu)
sens[,"n11"] <- eval(D(lambda, "n11"), envir=mu)

# Calculate elasticities
elas <- matrix(NA, n.draws, 27)  
colnames(elas) <- c("sf1", "sf2", "sf3", "sm1", "sm2", "sm3", "eta1", "eta2", "eta3", "rho1", "rho2", "rho3", "xi1", "xi2", "xi3", "im", "n1", "n2", "n3", "n4", "n5", "n6", "n7", "n8", "n9", "n10", "n11")
elas[,"sf1"] <- sens[,"sf1"] * mu$sf1 / m_lam
elas[,"sf2"] <- sens[,"sf2"] * mu$sf2 / m_lam
elas[,"sf3"] <- sens[,"sf3"] * mu$sf3 / m_lam
elas[,"sm1"] <- sens[,"sm1"] * mu$sm1 / m_lam
elas[,"sm2"] <- sens[,"sm2"] * mu$sm2 / m_lam
elas[,"sm3"] <- sens[,"sm3"] * mu$sm3 / m_lam
elas[,"eta1"] <- sens[,"eta1"] * mu$eta1 / m_lam
elas[,"eta2"] <- sens[,"eta2"] * mu$eta2 / m_lam
elas[,"eta3"] <- sens[,"eta3"] * mu$eta3 / m_lam
elas[,"rho1"] <- sens[,"rho1"] * mu$rho1 / m_lam
elas[,"rho2"] <- sens[,"rho2"] * mu$rho2 / m_lam
elas[,"rho3"] <- sens[,"rho3"] * mu$rho3 / m_lam
elas[,"xi1"] <- sens[,"xi1"] * mu$xi1 / m_lam
elas[,"xi2"] <- sens[,"xi2"] * mu$xi2 / m_lam
elas[,"xi3"] <- sens[,"xi3"] * mu$xi3 / m_lam
elas[,"im"] <- sens[,"im"] * mu$im / m_lam
elas[,"n1"] <- sens[,"n1"] * mu$n1 / 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
elas[,"n9"] <- sens[,"n9"] * mu$n9 / m_lam
elas[,"n10"] <- sens[,"n10"] * mu$n10 / m_lam
elas[,"n11"] <- sens[,"n11"] * mu$n11 / m_lam


# Define matrix to store results
cont <- matrix(NA, nrow=n.draws, ncol=27) 
colnames(cont) <- c("sf1", "sf2", "sf3", "sm1", "sm2", "sm3", "eta1", "eta2", "eta3", "rho1", "rho2", "rho3", "xi1", "xi2", "xi3", "im", "n1", "n2", "n3", "n4", "n5", "n6", "n7", "n8", "n9", "n10", "n11")

ra <- 1:(nyears-1)
# Calculate contributions for each demographic rate and stage-structured population size at each MCMC draw
for (s in 1:n.draws){
  dp_stoch <- cbind(out[s,i_s1F[ra]], out[s,i_s2F[ra]], out[s,i_s3F[ra]], out[s,i_s1M[ra]], out[s,i_s2M[ra]], out[s,i_s3M[ra]], out[s,i_eta1[ra]], out[s,i_eta2[ra]], out[s,i_eta3[ra]], out[s,i_rho1[ra]], out[s,i_rho2[ra]], out[s,i_rho3[ra]], out[s,i_xi1[ra]], out[s,i_xi2[ra]], out[s,i_xi3[ra]], immRate[s,ra], n1[s,ra], n2[s,ra], n3[s,ra], n4[s,ra], n5[s,ra], n6[s,ra], n7[s,ra], n8[s,ra], n9[s,ra], n10[s,ra], n11[s,ra])
  # 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)
}


rcont <- cont / rowSums(cont)

# Save the results
save(cont, rcont, sens, elas, file='ContributionFigure.Rdata')



# Contribution through demographic rates
mean(apply(abs(cont[,1:16]), 1, sum) / apply(abs(cont), 1, sum))
quantile(apply(abs(cont[,1:16]), 1, sum) / apply(abs(cont), 1, sum), c(0.025, 0.975))

# Contribution through population structure
mean(apply(abs(cont[,17:27]), 1, sum) / apply(abs(cont), 1, sum))
quantile(apply(abs(cont[,17:27]), 1, sum) / apply(abs(cont), 1, sum), c(0.025, 0.975))


# Sum the contributions
mean(apply(cont, 1, sum))

# Compute the 'annual environmental growth rate' without and with demographic stochasticity (cf. Knape et al. 2023)
ndet <- lam.real.ES <- lam.real <- matrix(NA, nrow=nrow(out), ncol=46)
for (t in 1:46){
  ndet[,t] <- out[,i_N2F[t]] * out[,i_xi1[t]] * out[,i_eta1[t]] * out[,i_rho1[t]] * out[,i_s1F[t]] + out[,i_N4F[t]] * out[,i_xi2[t]] * out[,i_eta2[t]] * out[,i_rho2[t]] * out[,i_s1F[t]] + out[,i_N5F[t]] * out[,i_xi3[t]] * out[,i_eta3[t]] * out[,i_rho3[t]] * out[,i_s1F[t]] + 
  out[,i_N2F[t]] * (1-out[,i_xi1[t]]) * out[,i_eta1[t]] * out[,i_rho1[t]] * out[,i_s1M[t]] + out[,i_N4F[t]] * (1-out[,i_xi2[t]]) * out[,i_eta2[t]] * out[,i_rho2[t]] * out[,i_s1M[t]] + out[,i_N5F[t]] * (1-out[,i_xi3[t]]) * out[,i_eta3[t]] * out[,i_rho3[t]] * out[,i_s1M[t]] + 
  (out[,i_N1F[t]] + out[,i_N2F[t]]) * out[,i_s2F[t]] + (out[,i_N3F[t]] + out[,i_N4F[t]] + out[,i_N5F[t]]) * out[,i_s3F[t]] + (out[,i_N1M[t]] + out[,i_N2M[t]]) * out[,i_s2M[t]] + (out[,i_N3M[t]] + out[,i_N4M[t]] + out[,i_N5M[t]]) * (out[,i_s3M[t]] + immRate[,t]) 
  lam.real.ES[,t] <- ndet[,t] / Ntot[,t]      # without dem. stoch
  lam.real[,t] <- Ntot[,t+1] / Ntot[,t]       # with dem. stoch
}


# Taylor approximation:
taylor.diag <- apply(cont, 1, sum) / apply(lam.real.ES, 1, var)
rel.error <- 1-taylor.diag
mean(rel.error)
quantile(rel.error, c(0.025, 0.975))


# Compute the 'demographic growth rate'
lam.real.D <- lam.real - lam.real.ES

# Compute the relative contribution of environmental variation to variation in the realized growth rate:
rel.env.var <- (var(lam.real.ES) + cov(lam.real.ES, lam.real.D)) / var(lam.real)

cov.ES.D <- numeric()
for (i in 1:nrow(out)){
  cov.ES.D[i] <- cov(lam.real.ES[i,], lam.real.D[i,])
}

rel.env.var <- (apply(lam.real.ES, 1, var) + cov.ES.D) / apply(lam.real, 1, var)

# Contribution env. stoch:
mean(rel.env.var)
quantile(rel.env.var, c(0.025, 0.975))



###########################################
#
# Assessment of density-dependence at the population level
#
###########################################

load('Model1.Rdata')
res_ipm <- res_ipm1
out <- do.call(rbind, res_ipm$samples)


BreederF <- out[,i_N2F] + out[,i_N4F] + out[,i_N5F] 
FloaterF <- out[,i_N1F] + out[,i_N3F]

N <- BreederF + FloaterF             

T <- ncol(N)
n <- nrow(N)

r <- sigma.r <- alpha.star <- beta.star <- sigma.star <- alpha <- beta <- sigma <- numeric(n)
N.star <- matrix(NA, nrow=n, ncol=T)
for (i in 1:n){    # loop over all MCMC draw
  if(round(i %% 500 ) == 0) {cat(paste('** Processing draw', i, '**\n'))}
  # 1. Estimate the growth parameters under the exponential model (i.e.,    without DD)
  z <- lm(log(N[i,-1] / N[i,-T]) ~ 1)
  r[i] <- coef(z)
  sigma.r[i] <- summary(z)$sigma
  # 2. Simulate population trajectories under the exponential growth model
  N.star[i,1] <- N[i,1]
  for (t in 1:(T-1)){
    N.star[i,t+1] <- N.star[i,t] * exp(r[i] + rnorm(1, 0, sigma.r[i]))
  } #t
  # 3. Fit the DD model to simulated population sizes from step 2 
  z <- lm(log(N.star[i,-1]) ~ log(N.star[i,-T]))
  alpha.star[i] <- coef(z)[1]
  beta.star[i] <- coef(z)[2]-1
  sigma.star[i] <- summary(z)$sigma
  # 4. Fit DD model to estimated population sizes (from the IPM)
  z <- lm(log(N[i,-1]) ~ log(N[i,-T]))
  alpha [i] <- coef(z)[1]
  beta[i] <- coef(z)[2]-1
  sigma[i] <- summary(z)$sigma
} #i

# Test for DD based on the total female population
mean(beta < beta.star)

# Save result file (for producing figure 6)
save(beta, beta.star, file='DDFigure.Rdata')