Commit eecb65cd authored by matbuoro's avatar matbuoro
Browse files

ajout txt

parent c299779d
######################################################################################
## Cm_B[t,a]: Annual number of marked fish captured at Beauchamps per sea age category; 1: 1SW, 2: MSW
## Cum_B[t,a]: Annual number of unmarked fish captured at Beauchamps per sea age category; 1: 1SW, 2: MSW
######################################################################################
Cm_B[,1] Cm_B[,2] Cum_B[,1] Cum_B[,2]
6 0 2 1
20 1 4 0
9 1 5 1
6 3 3 1
0 1 1 0
NA NA NA NA
2 0 1 0
7 0 3 0
5 2 0 0
NA NA NA NA
1 0 1 1
2 1 15 0
3 3 3 2
1 0 0 0
8 2 4 1
1 0 1 5
NA NA NA NA
NA NA NA NA
13 5 12 4
0 1 0 0
11 2 7 1
37 2 24 3
13 3 6 6
1 0 5 0
5 0 4 1
12 4 1 0
10 2 3 3
15 1 0 0
7 1 2 0
5 0 3 3
3 2 3 3
12 2 0 2
7 1 4 0 #2015/2016
#END
#####################################################################################################
## C_Eu[t,a]: Annual number of fish captured at Eu per sea age category; 1:1SW, 2:MSW
## Cm_Eu[t,a]: Annual number of fish captured and marked at Eu per sea age category; 1:1SW, 2:MSW
#####################################################################################################
C_Eu[,1] C_Eu[,2] Cm_Eu[,1] Cm_Eu[,2]
43 26 43 25
78 19 78 19
119 19 119 19
98 25 97 25
53 12 53 12
121 20 121 20
50 33 50 33
118 19 118 18
123 35 122 35
33 16 33 16
26 5 26 5
23 1 23 1
18 16 18 16
24 10 24 10
177 10 177 10
12 13 12 13
7 4 7 4
8 1 8 1
56 13 56 13
13 10 13 10
45 9 45 9
159 30 158 30
73 66 73 66
31 13 31 13
78 3 78 3
71 34 69 34
105 24 105 24
123 24 123 24
97 24 97 24
93 42 93 42
64 28 64 28
106 33 106 33
116 15 116 15 #2016
#END
######################################################################################################################################
## Mean flow (l/s): - 15 june - 31 august for 1SW ([,1])
## - 15 april - 30 june for MSW ([,2])
## Covariate is standardized within WinBUGS
######################################################################################################################################
Q[,1] Q[,2]
5858.292117 6922.272878
6449.612232 7809.483502
5686.187501 7141.983616
7053.504676 8180.862324
7276.216404 9886.532623
6097.252582 7687.061265
4431.323747 5444.970317
4833.538045 5482.957108
4834.534743 5182.577258
5338.2195 6115.544548
7650.689918 10344.17178
6545.820166 9489.651519
3918.415905 4360.736871
4119.325152 4633.25528
4839.929364 5276.819081
6074.813107 8034.131962
7921.282051 9988.311688
9638.205128 13543.76623
7821.153846 10150.64935
6276.153846 7725.714286
5354.487179 5401.688312
4998.205128 5231.948052
5112.948718 5460
6130 6144.545455
7714.358974 9162.597403
5772.820513 6832.207792
5724.871795 7065.974026
4707.051282 5039.350649
6508.333333 7241.168831
7005 8451.428571
7907.307692 8912.857143
6358.846154 7587.532468
7754.615385 8476.883117 #2016
#END
#####################################################################################################
## Y: length of the time series (number of years)
## Q2pic[t]: Annual mean flow associated to the second peak of migration (1 octobre - 30 november)
## Flow data (l/s) are standardized within the model.
#####################################################################################################
list(Y=33,
Q2pic = c(
6466.290,
5326.312,
5401.795,
8410.899,
5943.663,
4929.979,
4399.319,
4509.769,
6260.526,
6037.201,
6168.087,
5038.144,
4073.442,
4339.881,
6523.104,
5206.278,
9596.557,
8615.902,
7455.902,
5324.754,
4621.803,
4770.492,
4761.803,
5714.590,
7377.869,
5544.590,
5962.131,
4414.918,
7025.738,
8210.000,
6952.295,
6001.967,
6428.361
)
)
list(
# METTRE A JOUR
lambda_tot = c(
101.8,137.8,207.4,185.3,129.9,
207.0,84.98,198.1,165.3,136.1,
71.26,135.3,65.63,66.28,279.3,
87.58,48.08,131.9,128.1,84.58,
161.4,357.2,134.2,151.9,180.1,
117.4,184.8,146.0,198.6,179.2,
180, 180),
# METTRE A JOUR /!\ TAILLE MATRICE
logit_pi_B = structure(.Data = c(
-2.496,-3.38,-1.372,-2.32,-2.6,
-2.831,-2.956,-2.549,-3.518,-3.105,
NA, NA,-2.793,-3.107,-3.062,
-2.998,-2.493,-3.884, NA, NA,
-3.226,-2.017,-2.034,-1.316,-2.486,
-1.573,-3.139,-3.074,-2.49,-2.94,
-2.965,-1.942, NA, NA, NA,
NA,-1.371,-1.124,-3.252,-3.075,
-2.032,-2.273,-1.121,-1.827,-1.782,
-2.796,-2.925,-3.714,-2.054,-1.113,
-1.612,-2.486,-2.649,-2.862,-2.107,
-2.132,-2.473,-2.559,-2.798,-2.454,
-2.132,-2.473,-2.559,-2.798),
.Dim = c(32,2)),
# METTRE A JOUR /!\ TAILLE MATRICE
logit_pi_Eu = structure(.Data = c(
0.343,1.078,1.487,0.906,0.9426,
1.136,0.2599,0.8924,-0.0344,1.519,
1.249,1.272,1.719,0.7151,1.138,
0.4976,1.947,2.237,-0.7839,0.2312,
-0.4747,0.7478,-1.035,0.3989,-1.214,
1.034,0.02703,1.256,1.021,0.9162,
-0.8177,0.1278,-0.742,-1.017,-1.053,
-2.073,0.2711,-0.2967,-1.547,0.4636,
-0.6245,1.732,0.5033,0.681,0.3818,
0.8959,-1.375,0.4715,0.4871,-0.9009,
1.564,1.294,1.121,0.4898,2.645,
1.779,1.479,1.005,0.09753,0.7002,
1.779,1.479,1.005,0.09753),
.Dim = c(32,2)),
# METTRE A JOUR /!\ TAILLE MATRICE
n = structure(.Data = c(
74.0,33.0,104.0,23.0,165.0,
27.0,168.0,41.0,97.0,13.0,
163.0,29.0,55.0,45.0,161.0,
25.0,139.0,35.0,120.0,37.0,
59.0,14.0,115.0,2.0,58.0,
21.0,49.0,13.0,240.0,14.0,
47.0,31.0,35.0,21.0,109.0,
14.0,112.0,24.0,65.0,22.0,
124.0,12.0,249.0,50.0,116.0,
102.0,129.0,19.0,129.0,10.0,
85.0,40.0,157.0,40.0,129.0,
28.0,120.0,33.0,147.0,69.0,
120.0,33.0,147.0,69.0),
.Dim = c(32,2)),
# NO UPDATE
lambda_tot0 = 104.2,
logit_flow_Eu = c(
-0.4262,-0.7079),
lflow_fall_Eu = c(-0.5,-0.5),
logit_int_Eu = c(
0.5082,0.8304),
mupi_B = c(
0.0827,0.1032),
pi_Eu00 = c(
0.5871,0.5992),
pi_Eu01 = c(
0.2006,0.7994),
rate_lambda = 0.04238,
s = c(
15.67,3.914),
shape_lambda = 5.995,
sigmapi_B = c(
0.644,0.6969),
sigmapi_Eu = c(
1.019,0.8865)
)
################################################################################
### Model of CMR data to estimate spawners population size ###
### of Salmo salar in Bresle river. ###
### Sabrina Servanty & Etienne Prévost ###
### January 2015 ###
################################################################################
model {
############################################################################################
## Difference with Sébastien's report (Delmotte et al. 2010):
## - set probability of capture at Eu in 2000 and 2001 to be smaller than the other years (partial trapping)
## - didn't consider recapture probability in 1989 and 1993 because trap in Beauchamps was not working those years (and so in 2000 and 2001 as in the report)
## - adding a flow effect in pi_Eu. Considering a different temporal window depending on sea age. Flow data are standardized (ln(Q[t]) - mean(ln(Q)))/sd(ln(Q)) within the model. Residuals are standardized and followed.
########### - 15 june - 31 august for 1SW
########### - 15 april - 30 june for MSW
## - adding another effect of flow corresponding to the second peak of migration.Flow data are standardized within the model.
########### Same temporal window for 1SW and MSW: 1 octobre - 30 november
## - calculating calculating R² (% of variation explained by the two covariates in the probability of capture at Eu)
###############################################################################################
############################################################################################
## Used indices:
## t: year; 1 to Y - from 1994 to Y ##
## a: sea age;
## 1-1SW (Grisle),
## 2-MSW (salmon)
#####################################################################################################
#####################################################################################################
## DATA
## --------------
## Y: length of the time series (number of years)
## C_Eu[t,a]: Annual number of fish captured at Eu per sea age category (1SW/MSW)
## Cm_Eu[t,a]: Annual number of fish captured and marked at Eu per sea age category
## Cm_B[t,a]: Annual number of marked fish captured at Beauchamps per sea age category
## Cum_B[t,a]: Annual number of unmarked fish captured at Beauchamps per sea age category
## Q[t,a]: Annual mean flow associated to sea age category
## Q2pic[t]: Annual mean flow associated to the second peak of migration
######################################################################################
#################################################################################
# PROBABILITIES (p)
# ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
## HYPER-PARAMETERS
## ------
#################################################################################
### Mean and standard deviation of probabilities depending on sea age.
# Probabilities to be captured at Eu.
for (a in 1:2) {
pi_Eu00[a] ~ dbeta(1,1) # year 2000 is considered to be different (partial trapping)
pi_Eu01[a] ~ dbeta(1,1) # year 2001 is considered to be different (partial trapping)
logit_int_Eu[a] ~ dunif(-10,10) #intercept
logit_flow_Eu[a] ~ dunif(-10,10) #slope for flow data (1SW at Eu: 15 june - 31 august, MSW: 15 april - 30 june)
sigmapi_Eu[a] ~ dunif(0,20)
lflow_fall_Eu[a] ~ dunif(-10,10) # slope for flow data, second peak of migration (same period for both 1SW and MSW)
# Probability to be captured at Beauchamps (after reproduction)
mupi_B[a] ~ dbeta(1,1)
sigmapi_B[a] ~ dunif(0,20)
}
################################################################################
## PROBABILITY DISTRIBUTIONS
## -------------------------
## pi_Eu[t,a]:annual probability to be captured at Eu given sea age
## pi_S[t,a]: annual survival until reproduction given sea age
## pi_B[t,a]: annual probability to be captured at Beauchamps given sea age (after reproduction)
####################################################################################
### Probabilities to be captured at Eu (time and sea age dependent)
for (t in 1:Y) {
lQ2pic[t] <- log(Q2pic[t]) # ln transformation of autumn covariate
stlQ2pic[t] <- (lQ2pic[t] - mean(lQ2pic[]))/sd(lQ2pic[]) # standardized covariate
} #end of loop over years
for (a in 1:2) {
varpi_Eu[a] <- (sigmapi_Eu[a])*(sigmapi_Eu[a])
precpi_Eu[a] <- 1/(varpi_Eu[a]) # precision
for (t in 1:Y) { #pi_Eu exchangeable from 1984 to 1999
logQ[t,a] <- log(Q[t,a]) # ln transformation of covariate
stlogQ[t,a] <- (logQ[t,a] - mean(logQ[,a]))/sd(logQ[,a]) # standardized covariate
logit_mupi_Eu[t,a] <- logit_int_Eu[a] + logit_flow_Eu[a] * stlogQ[t,a] + lflow_fall_Eu[a] * stlQ2pic[t]
logit_pi_Eu[t,a] ~ dnorm(logit_mupi_Eu[t,a],precpi_Eu[a])
pi_Eu[t,a] <- exp(logit_pi_Eu[t,a])/(1+exp(logit_pi_Eu[t,a])) # back-transformation on the probability scale
epsilon_Eu[t,a] <- (logit_pi_Eu[t,a] - logit_mupi_Eu[t,a]) / sigmapi_Eu[a] # standardized residuals
eps_Eu[t,a] <- logit_pi_Eu[t,a] - logit_mupi_Eu[t,a] # residuals not standardized
} ## End of loop over years
#Calculating R² = 1 -(E(variance of residuals (/!\ not standardized!) / E(variance of capture probabilities)))
# See Gelman & Pardoe 2006
sdeps_Eu[a] <- sd(eps_Eu[,a])
vareps_Eu[a] <- sdeps_Eu[a] * sdeps_Eu[a]
sdlpi_Eu[a] <- sd(logit_pi_Eu[,a])
varlpi_Eu[a] <- sdlpi_Eu[a] * sdlpi_Eu[a]
R2[a] <- 1 - (mean(vareps_Eu[a])/mean(varlpi_Eu[a]))
}# end of loop over sea age
test[1] <- step(logit_flow_Eu[1]) # is logit_flow >=0 for 1SW?
test[2] <- step(logit_flow_Eu[2]) # is logit_flow >=0 for MSW?
test[3] <- step(lflow_fall_Eu[1]) # is logit_flow_fall >=0 for 1SW?
test[4] <- step(lflow_fall_Eu[2]) # is logit_flow_fall >=0 for 1SW?
diff_flow <- logit_flow_Eu[1] - logit_flow_Eu[2]
test[5] <- step(diff_flow) #is difference in slope >=0 for 1SW? (1SW>MSW)
diff_fall <- lflow_fall_Eu[1] - lflow_fall_Eu[2]
test[6] <- step(diff_fall) # is difference in slope in fall >=0 for 1SW? (1SW>MSW)
### Probabilities to be captured at Beauchamps after reproduction (time and sea age dependent)
for (a in 1:2) {
logit_mupi_B[a] <- log(mupi_B[a]/(1-mupi_B[a])) # logit transformation
varpi_B[a] <- (sigmapi_B[a])*(sigmapi_B[a])
precpi_B[a] <- 1/(varpi_B[a]) # precision
for (t in 1:5) { ## pi_B: Exchangeable from 1984 to 1988
logit_pi_B[t,a] ~ dnorm(logit_mupi_B[a],precpi_B[a])
pi_B[t,a] <- exp(logit_pi_B[t,a])/(1+exp(logit_pi_B[t,a])) # back-transformation on the probability scale
} ## End of loop over years
pi_B[6,a] <- 0 # Trap not working in 1989
for (t in 7:9) { ## pi_B: Exchangeable from 1990 to 1992
logit_pi_B[t,a] ~ dnorm(logit_mupi_B[a],precpi_B[a])
pi_B[t,a] <- exp(logit_pi_B[t,a])/(1+exp(logit_pi_B[t,a])) # back-transformation on the probability scale
} ## End of loop over years
pi_B[10,a] <- 0 # Trap not working in 1993
for (t in 11:16) { ## pi_B: Exchangeable from 1994 to 1999
logit_pi_B[t,a] ~ dnorm(logit_mupi_B[a],precpi_B[a])
pi_B[t,a] <- exp(logit_pi_B[t,a])/(1+exp(logit_pi_B[t,a])) # back-transformation on the probability scale
} ## End of loop over years
pi_B[17,a] <- 0 # Trap not working in 2000
pi_B[18,a] <- 0 # Trap not working in 2001
for (t in 19:Y) { ## pi_B: Exchangeable from 2002 to now on
logit_pi_B[t,a] ~ dnorm(logit_mupi_B[a],precpi_B[a])
pi_B[t,a] <- exp(logit_pi_B[t,a])/(1+exp(logit_pi_B[t,a])) # back-transformation on the probability scale
} ## End of loop over years
} # end of loop over sea age
######################################################################################################################################
######################################################################################################################################
## DATA
## --------------
## C_Eu[t,a]: Annual number of fish captured at Eu per sea age category (1SW/MSW)
## Cm_Eu[t,a]: Annual number of fish captured and marked at Eu per sea age category
## Cm_B[t,a]: Annual number of marked fish captured at Beauchamps per sea age category
## Cum_B[t,a]: Annual number of unmarked fish captured at Beauchamps per sea age category
######################################################################################
## NUMBERS PER GROUP
## |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
## n_tot[t]: annual total number of fish
## n_1SW[t]: annual Number of 1SW
## n_MSW[t]: annual Number of MSW
## We included a hierarchy on breeding category.
####################################################
## PRIORS
## ------
####################################################################################
## Annual returns
for (t in 1:Y) {
n_tot[t] <- sum(n[t,]) # total number
n_1SW[t] <- n[t,1] # Number of males and females 1SW
n_MSW[t] <- n[t,2] # Number of males and females MSW
} ## End of loop over years
# Shape and rate parameter for gamma distribution for negative binomial (see Gelman, 2d edition, p446)
shape_lambda ~ dgamma(0.001,0.001)
rate_lambda ~ dgamma(0.001,0.001)
#To be able to initiate lambda_tot ensuite
lambda_tot0 ~ dgamma(shape_lambda,rate_lambda)
Plambda0[1:2] ~ ddirich(s[])
# Hyperprior for lambda_n
for (t in 1:Y) {
lambda_tot[t] ~ dgamma(shape_lambda,rate_lambda)
Plambda[t,1:2] ~ ddirich(s[])
} # end of loop over years
for (a in 1:2) {
s[a] ~ dexp(0.25)T(0.02,)
} # end of loop over sea age
for (t in 1:Y-1) {
lambda_n[t,1] <- lambda_tot[t] * Plambda[t,1]
lambda_n[t+1,2] <- lambda_tot[t] * Plambda[t,2]
} # end of loop over years
### lambda_n is calculated by smolt cohort (proportion and mean population size of MSW is one year yearlier than 1SW)
## First year for MSW
lambda_n[1,2] <- lambda_tot0 * Plambda0[2]
## Last year for 1SW
lambda_n[Y,1] <- lambda_tot[Y] * Plambda[Y,1]
for (a in 1:2) {
for (t in 1:Y) {
n[t,a] ~ dpois(lambda_n[t,a]) # annual number of fish per sea age
} #end of loop over years
} # end of loop over sea age
#################################
## CAPTURE AT EU
## -------
## n_um[t]: Annual number of fish not captured (not marked) at Eu
##################################
## Likelihood on trapped fish number
for (a in 1:2) {
for (t in 1:16) { # from 1984 to 1999
C_Eu[t,a] ~ dbin(pi_Eu[t,a],n[t,a])
## Fish number avoiding the trap, not marked fish (n_um).
n_um[t,a] <- n[t,a] - C_Eu[t,a]
} ## End of loop over years
#Year 2000
p_Eu00_tot[a] <- pi_Eu[17,a] * pi_Eu00[a] # probability in 2000 to be observed is set to be smaller than later
C_Eu[17,a] ~ dbin(p_Eu00_tot[a],n[17,a]) # 1994 is considered to be different.
n_um[17,a] <- n[17,a] - C_Eu[17,a] ## Fish number avoiding the trap, not marked fish (n_um) in 2000
#Year 2001
p_Eu01_tot[a] <- pi_Eu[18,a] * pi_Eu01[a] # probability in 2000 to be observed is set to be smaller than later
C_Eu[18,a] ~ dbin(p_Eu01_tot[a],n[18,a]) # 1994 is considered to be different.
n_um[18,a] <- n[18,a] - C_Eu[18,a] ## Fish number avoiding the trap, not marked fish (n_um) in 2001
for (t in 19:Y) { # from 2002 to now on
C_Eu[t,a] ~ dbin(pi_Eu[t,a],n[t,a])
## Fish number avoiding the trap, not marked fish (n_um).
n_um[t,a] <- n[t,a] - C_Eu[t,a]
} ## End of loop over years
} # end of loop over sea age
#####################################
## RECAPTURE OF FISH AT BEAUCHAMPS AFTER REPRODUCTION
## -------------------------------------
##
for (a in 1:2) {
for (t in 1:Y) {
## Marked fish
Cm_B[t,a] ~ dbin(pi_B[t,a],Cm_Eu[t,a])
## Unmarked fish
Cum_B[t,a] ~ dbin(pi_B[t,a],n_um[t,a])
} # end of loop over years
} # end of loop over sea age
} # end of the model
\ No newline at end of file
=============================
DIAGNOSTICS
=============================
---------------------------
Heidelberger and Welch's convergence diagnostic
heidel.diag is a run length control diagnostic based on a criterion of relative accuracy for the estimate of the mean. The default setting corresponds to a relative accuracy of two significant digits.
heidel.diag also implements a convergence diagnostic, and removes up to half the chain in order to ensure that the means are estimated from a chain that has converged.
Stationarity start p-value
test iteration
shape_lambda passed 1 0.976
rate_lambda passed 1 0.952
lambda_tot0 passed 1 0.624
Halfwidth Mean Halfwidth
test
shape_lambda passed 3.604 0.051319
rate_lambda passed 0.024 0.000343
lambda_tot0 passed 167.980 1.800868
---------------------------
Geweke's convergence diagnostic
Geweke (1992) proposed a convergence diagnostic for Markov chains based on a test for equality of the means of the first and last part of a Markov chain (by default the first 10% and the last 50%).
If the samples are drawn from the stationary distribution of the chain, the two means are equal and Geweke's statistic has an asymptotically standard normal distribution.
The test statistic is a standard Z-score: the difference between the two sample means divided by its estimated standard error. The standard error is estimated from the spectral density at zero and so takes into account any autocorrelation.
The Z-score is calculated under the assumption that the two parts of the chain are asymptotically independent, which requires that the sum of frac1 and frac2 be strictly less than 1.
Fraction in 1st window = 0.1
Fraction in 2nd window = 0.5
shape_lambda rate_lambda lambda_tot0
0.1124 -0.2537 1.4104
---------------------------
Raftery and Lewis's diagnostic
Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95
Burn-in Total Lower bound Dependence
(M) (N) (Nmin) factor (I)
shape_lambda 20 26665 3746 7.12
rate_lambda 20 24844 3746 6.63
lambda_tot0 12 16041 3746 4.28
=============================
DIAGNOSTICS
=============================
---------------------------
Heidelberger and Welch's convergence diagnostic
heidel.diag is a run length control diagnostic based on a criterion of relative accuracy for the estimate of the mean. The default setting corresponds to a relative accuracy of two significant digits.
heidel.diag also implements a convergence diagnostic, and removes up to half the chain in order to ensure that the means are estimated from a chain that has converged.
Stationarity start p-value
test iteration
shape_lambda passed 1 0.0754
rate_lambda passed 1 0.0988
lambda_tot0 passed 1 0.9441
Halfwidth Mean Halfwidth
test
shape_lambda passed 3.6837 0.128060
rate_lambda passed 0.0243 0.000869
lambda_tot0 passed 164.1569 4.714110