### ---------------------  APPENDIX 2
### JAGS model code used to model freshwater fish population dynamics for two different size-classes
###
### Variables' names :
### - ln.Xad = matrix of log-abundance of >0+ individuals with time series (i.e. sites) in column and years in rows
### - ln.X0 = matrix of log-abundance of 0+ individuals with time series (i.e. sites) in column and years in rows
### - time = vector of the time series length
### - var_temp = matrix of temperature variability with time series (i.e. sites) in column and years in rows
### - mu_temp = matrix of average temperature with time series (i.e. sites) in column and years in rows
### - samp = matrix of sampling area with time series (i.e. sites) in column and years in rows
### - Nts =  total number of time series
### - Y = vector of latitude
### - Alt = vector of altitude (m)
### - Dist = vector of distance to the geographic center (km)
###
### Parameters' names:
### - lambda_ad (and lambda_juv) = estimated log-abundance of >0+ (and 0+) individuals 
### - tau.proc.ad (and tau.proc.juv) = precision (inverse of variance) of the normal distributions -> process error variance around estimated log-abundances
### - inter.ad (and inter.juv) = site-specific intercepts
### - mu.inter.ad (mu.inter.juv) = global intercepts 
### - tau.inter.ad (and tau.inter.juv) = precision (inverse of variance) of site intercepts
### - beta.ad = coefficients tailored to the dynamic of >0+ individuals (1=linear effect of average T°C, 2=linear effect of T°C variability, 
###                                                                      3=quadratic effect of average T°C, 4=quadratic effect of T°C variability, 
###                                                                      5=effect of the density of 0+ individuals at t-1[transition probability between size classes], 
###                                                                      6=effect of the density of >0+ individuals at t-1[density dependence])
### - beta.juv = coefficients tailored to the dynamic of 0+ individuals (1=linear effect of average T°C, 2=linear effect of T°C variability, 
###                                                                      3=quadratic effect of average T°C, 4=quadratic effect of T°C variability, 
###                                                                      5=effect of the density of >0+ individuals at t[recruitment rate])

model{
    
    ### LIKELIHOOD
	
	#-- Population dynamic provess
    for(ts in 1:Nts){
      
      for(t in 2:time[ts]){
	  
        #ADULTS
        ln.Xad[t,ts] ~ dnorm(lambda_ad[t,ts], tau.proc.ad)	
		
        lambda_ad[t,ts] <- inter.ad[ts] + mu.inter.ad + ln.Xad[t-1,ts] +
								beta.ad[1, ts]*mu_temp[t,ts] + beta.ad[2,ts]*var_temp[t,ts] + 
								beta.ad[3, ts]*pow(mu_temp[t,ts],2) + beta.ad[4,ts]*pow(var_temp[t,ts],2) + 
								beta.ad[5, ts]*(ln.X0[t-1,ts]/log(samp[t-1,ts])) + # Recruitment to the >0+ class
								beta.ad[6, ts]*(ln.Xad[t-1,ts]/log(samp[t-1,ts])) + # Gompertz type of DD		
								log(samp[t,ts]/samp[t-1,ts])
    
        #JUVENILES
        ln.X0[t,ts] ~ dnorm(lambda_juv[t,ts], tau.proc.juv)
		
		lambda_juv[t,ts] <- inter.juv[ts] + mu.inter.juv +
								 beta.juv[1, ts]*mu_temp[t,ts] + beta.juv[2, ts]*var_temp[t,ts] + 
								 beta.juv[3, ts]*pow(mu_temp[t,ts],2) + beta.juv[4, ts]*pow(var_temp[t,ts],2) +
								 beta.juv[5, ts]*(ln.Xad[t,ts]/log(samp[t,ts])) + log(samp[t,ts])

      }

      inter.ad[ts] ~ dnorm(0, tau.inter.ad)
      inter.juv[ts] ~ dnorm(0, tau.inter.juv)
	  for(z1 in 1:6){beta.ad[z1, ts] ~ dnorm(mean.beta.ad[z1], tau.beta.ad[z1])}
	  for(z1 in 1:5){beta.juv[z1, ts] ~ dnorm(mean.beta.juv[z1], tau.beta.juv[z1])}	  
		
    }

    ### PRIORS    
    tau.proc.ad <- pow(sd.proc.ad,-2)
    sd.proc.ad ~ dt(0,1,1)T(0,)

    tau.proc.juv <- pow(sd.proc.juv,-2)
    sd.proc.juv ~ dt(0,1,1)T(0,)

    mu.inter.ad ~ dnorm(0, 0.1)
    tau.inter.ad <- pow(sd.inter.ad,-2)
    sd.inter.ad ~ dt(0,1,1)T(0,)

    mu.inter.juv ~ dnorm(0, 0.1)
    tau.inter.juv <- pow(sd.inter.juv,-2)
    sd.inter.juv ~ dt(0,1,1)T(0,)

    for(z1 in 1:6){
		mean.beta.ad[z1] ~ dnorm(0, 0.1) 
		tau.beta.ad[z1] <- pow(sd.beta.ad[z1], -2)
		sd.beta.ad[z1] ~ dt(0,1,1)T(0,)
	}
	for(z1 in 1:5){
		mean.beta.juv[z1] ~ dnorm(0, 0.1) 
		tau.beta.juv[z1] <- pow(sd.beta.juv[z1], -2)
		sd.beta.juv[z1] ~ dt(0,1,1)T(0,)
	}
	
    ### Posterior predictive checks on standardized residuals for both adults and juveniles
    for(ts in 1:Nts){
	
      for(t in 2:time[ts]){
	  
        # Discrepancy for observed data
        res.ad.1[t,ts] <- (ln.Xad[t,ts]-lambda_ad[t,ts])/(sd.proc.ad)
		res.juv.1[t,ts] <- (ln.X0[t,ts]-lambda_juv[t,ts])/(sd.proc.juv)

        # Discrepancy for replicate data
        ln.X0.new[t,ts] ~ dnorm(lambda_juv[t,ts], tau.proc.juv)
        ln.Xad.new[t,ts] ~ dnorm(lambda_ad[t,ts], tau.proc.ad)
        res.ad.2[t,ts] <- (ln.Xad.new[t,ts]-lambda_ad[t,ts])/(sd.proc.ad)
		res.juv.2[t,ts] <- (ln.X0.new[t,ts]-lambda_juv[t,ts])/(sd.proc.juv)
      }
	  
      tmp.fit.ad[ts] <- sum(pow(res.ad.1[2:time[ts],ts],2))
      tmp.fit.new.ad[ts] <- sum(pow(res.ad.2[2:time[ts],ts],2))
	  
	  tmp.fit.juv[ts] <- sum(pow(res.juv.1[2:time[ts],ts],2))
      tmp.fit.new.juv[ts] <- sum(pow(res.juv.2[2:time[ts],ts],2))
    }
	
    fit.ad <- sum(tmp.fit.ad[])
    fit.new.ad <- sum(tmp.fit.new.ad[])
	test.ad <- step(fit.new.ad-fit.ad)
    bpvalue.ad <- mean(test.ad)
	
	fit.juv <- sum(tmp.fit.juv[])
    fit.new.juv <- sum(tmp.fit.new.juv[])
	test.juv <- step(fit.new.juv-fit.juv)
    bpvalue.juv <- mean(test.juv)
}