Commit 68da8925 authored by matbuoro's avatar matbuoro
Browse files

Update 2016

parent 7a3b17ab
......@@ -44,7 +44,8 @@ if(!file.exists(paste("inits/init-",site,"-",stade,year,".txt",sep=""))){
#load(paste('inits/inits_',stade,'.Rdata',sep="")) # chargement des inits
#if(site == "Bresle" && stade == "adult") {inits <- list(read.bugsdata(paste("inits/init-",site,"-",stade,year,".txt",sep="")))}
#if(site == "Nivelle") {inits <- list(read.bugsdata(paste("inits/init-",site,"-",stade,year,".txt",sep="")))}
inits <- list(read.bugsdata(paste("inits/init-",site,"-",stade,year,".txt",sep="")))
inits.tmp <- read.bugsdata(paste("inits/init-",site,"-",stade,year,".txt",sep=""))
inits <- rep(list(inits.tmp),2)
#------------------------MODEL----------------------------------##
model <- paste("model/model_",stade,"-",site,".R",sep="") # path of the model
......@@ -56,11 +57,11 @@ filename <- file.path(work.dir, model)
#---------------------------ANALYSIS-----------------------------##
nChains = length(inits) # Number of chains to run.
nChains = 2 #length(inits) # Number of chains to run.
adaptSteps = 1000 # Number of steps to "tune" the samplers.
nburnin=5000 # Number of steps to "burn-in" the samplers.
nstore=25000 # Total number of steps in chains to save.
nthin=2 # Number of steps to "thin" (1=keep every step).
nstore=10000 # Total number of steps in chains to save.
nthin=5 # Number of steps to "thin" (1=keep every step).
#nPerChain = ceiling( ( numSavedSteps * thinSteps ) / nChains ) # Steps per chain.
### Start of the run ###
......@@ -139,7 +140,7 @@ cat("=============================\n")
if (nChains > 1) {
cat("Convergence: gelman-Rubin R test\n")
gelman.diag(fit.mcmc[,which(varnames(fit.mcmc)%in%parameterstotest)])
gelman.diag(fit.mcmc[,which(varnames(fit.mcmc)%in%parameterstotest)],multivariate=TRUE)
}
......@@ -177,6 +178,7 @@ pdf(paste('results/Results_',stade,"_",year,'.pdf',sep=""))
traplot(fit.mcmc[,which(varnames(fit.mcmc)%in%parameterstotest)])
# caterplot(fit.mcmc,parameters[i])
#}
gelman.plot(fit.mcmc[,which(varnames(fit.mcmc)%in%parameterstotest)])
dev.off()
......
list(shape_lambda=5.00000E+00, rate_lambda=1.00000E-02, lambda_tot0=2.89400E+02, a_1.1SW=4.20000E+00, a_1.1SW=8.00000E-02, a_MSW=c(4.00000E-01, 9.05900E-01, 2.51600E+00, 4.55000E-01, 3.03700E-01), eta_1=1.35000E+00, eta_2=5.50000E+00, k_1=8.70000E-01, k_2=1.46000E+00, logit_p_11_1=c(-1.42000E+00, -1.42000E+00, -1.42000E+00, -1.42000E+00, -1.42000E+00, -1.42000E+00, -1.42000E+00, -1.42000E+00), logit_pi_EF=c(-1.68000E+00, -1.68000E+00, -1.68000E+00, -1.68000E+00, -1.68000E+00, -1.68000E+00, -1.68000E+00, -1.68000E+00), mup_11_1=2.10000E-01, mup_11_2=1.80000E-01, mup_21=4.20000E-01, mupi_EF=1.60000E-01, mupi_U=c(8.68900E-01, 8.14600E-01), rho=7.90000E-01, sigmap_11_1=2.10000E-01, sigmap_11_2=9.00000E-02, sigmap_12=c(5.50000E-01, 4.10000E-01, 2.90000E-01, 4.00000E-01), sigmap_21=4.30000E-01, sigmapi_EF=5.20000E-01, sigmapi_U=c(7.34100E-01, 5.64100E-01), alpha_1=c(2.96000E+00, 1.06000E+00, 1.07000E+00, 6.30000E-01, 1.09000E+00, 3.10000E-01, 3.90000E-01, 5.80000E-01, 2.10000E-01, 7.90000E-01, 6.80000E-01, 5.60000E-01, 5.20000E-01, 7.10000E-01, 7.20000E-01, 4.30000E-01, 7.70000E-01, 4.80000E-01, 1.47000E+00, 1.50000E+00, 1.39000E+00, 1.26000E+00, 1.37000E+00, 1.24000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00), alpha_2=c( NA, NA, NA, NA, NA, NA, 8.10000E-01, 9.70000E-01, 1.18000E+00, 2.08000E+00, 1.77000E+00, 1.22000E+00, 1.31000E+00, 1.62000E+00, 1.50000E+00, 1.41000E+00, 1.00000E+00, 9.30000E-01, 2.76000E+00, 2.13000E+00, 2.44000E+00, 1.67000E+00, 2.28000E+00, 2.22000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00), e_21= structure(.Data= c( NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 2.00000E+01, 2.30000E+01, 1.00000E+00, 3.00000E+00, 1.90000E+01, 1.10000E+01, 1.00000E+00, 6.00000E+00, 3.10000E+01, 2.80000E+01, 1.00000E+00, 7.00000E+00, 7.60000E+01, 6.70000E+01, 1.00000E+00, 6.00000E+00, 3.90000E+01, 3.20000E+01, 1.00000E+00, 6.00000E+00, 3.20000E+01, 3.00000E+01, 2.00000E+00, 3.00000E+00, 1.90000E+01, 3.80000E+01, 1.00000E+00, 1.20000E+01, 1.00000E+01, 2.20000E+01, 1.00000E+00, 3.00000E+00, 2.70000E+01, 2.60000E+01, 1.00000E+00, 4.00000E+00, 2.60000E+01, 2.60000E+01, 3.00000E+00, 4.00000E+00, 1.90000E+01, 1.40000E+01, NA, 3.00000E+00, 4.20000E+01, 3.00000E+01, 1.00000E+00, 8.00000E+00, 5.50000E+01, 6.00000E+01, 1.00000E+00, 8.00000E+00, 1.00000E+00, 5.00000E+00, 5.00000E+00, 1.00000E+01, 7.00000E+00, 5.00000E+00, NA, 2.00000E+00, 1.10000E+01, 1.10000E+01, NA, 3.00000E+00, 5.00000E+00, 3.00000E+00, 2.00000E+00, 5.00000E+00, 1.00000E+01, 7.00000E+00, 2.00000E+00, 3.00000E+00, 6.00000E+00, 5.00000E+00, NA, 7.00000E+00, 1.50000E+01, 6.00000E+00, NA, 6.00000E+00, 3.00000E+01, 2.50000E+01, 5.00000E+00, 5.00000E+00, 5.00000E+00, 5.00000E+00, 7.00000E+00, 8.00000E+00, 3.00000E+00, 9.00000E+00, 1.00000E+00, 3.00000E+00, 3.00000E+00, 9.00000E+00, 1.00000E+00, 3.00000E+00, 3.00000E+00, 9.00000E+00, 1.00000E+00, 3.00000E+00, 3.00000E+00, 2.00000E+00, 1.00000E+00, 3.00000E+00, 6.00000E+00, 4.00000E+00, 0.00000E+00, 5.00000E+00), .Dim=c(33, 4)), logit_p_11_2=c( NA, NA, NA, NA, NA, NA, NA, NA, -1.98000E+00, -1.98000E+00, -1.98000E+00, -1.98000E+00, -1.98000E+00, -1.98000E+00, -1.98000E+00, -1.98000E+00, -1.98000E+00, -1.98000E+00, -1.98000E+00, -1.98000E+00, -1.98000E+00, -1.98000E+00, -1.98000E+00, -1.98000E+00, -1.98000E+00, -1.98000E+00, -1.98000E+00, -1.98000E+00, -1.98000E+00, -1.98000E+00, -1.98000E+00, -1.98000E+00, -1.98000E+00), logit_p_21=c( NA, NA, NA, NA, NA, NA, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01, 2.80000E-01), logit_pi_U= structure(.Data= c(1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01, 1.77000E+00, 9.80000E-01), .Dim=c(33, 2)), logit_p_n12= structure(.Data= c( NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 7.14500E-01, 5.69500E-01, 3.75000E-01, 3.00900E-01, 8.80700E-01, 5.64500E-01, 3.83500E-01, 4.14300E-01, 6.82800E-01, 4.38800E-01, 7.57900E-01, 3.55800E-01, 9.00100E-01, 6.12300E-01, 4.70100E-01, 2.25500E-01, 6.79700E-01, 8.03500E-01, 5.83700E-01, 5.93600E-01, 6.25600E-01, 7.10600E-01, 6.62000E-01, 4.64900E-01, 6.67400E-01, 7.37100E-01, 6.38600E-01, 5.84000E-01, 8.12100E-01, 7.68100E-01, 6.27600E-01, 4.32600E-01, 6.75300E-01, 6.20000E-01, 5.34900E-01, 3.69400E-01, 8.18400E-01, 6.87500E-01, 4.17100E-01, 5.99000E-01, 5.36100E-01, 6.36200E-01, 4.27700E-01, 4.70300E-01, 6.01100E-01, 6.91900E-01, 7.27400E-01, 4.54800E-01, 4.88600E-01, 4.17400E-01, 3.24600E-01, 3.72300E-01, 7.50400E-01, 6.15400E-01, 5.45800E-01, 4.55200E-01, 5.78100E-01, 5.54400E-01, 6.49700E-01, 4.26600E-01, 7.84700E-01, 6.11400E-01, 6.46700E-01, 4.57700E-01, 5.89100E-01, 4.39400E-01, 3.23100E-01, 5.60600E-01, 7.60500E-01, 6.75400E-01, 5.46800E-01, 5.67300E-01, 7.69300E-01, 6.37300E-01, 7.64300E-01, 5.38200E-01, 7.23300E-01, 5.08200E-01, 8.75400E-01, 4.83400E-01, 6.22700E-01, 6.20400E-01, 6.58300E-01, 3.51300E-01, 5.04900E-01, 5.46700E-01, 4.61600E-01, 6.16200E-01, 5.04900E-01, 5.46700E-01, 4.61600E-01, 6.16200E-01, 5.04900E-01, 5.46700E-01, 4.61600E-01, 6.16200E-01, 5.04900E-01, 5.46700E-01, 4.61600E-01, 6.16200E-01), .Dim=c(33, 4)), n= structure(.Data= c(7.10000E+01, 1.01000E+02, 8.00000E+00, 1.50000E+01, 4.20000E+01, 3.50000E+01, 9.00000E+00, 3.80000E+01, 1.54000E+02, 1.58000E+02, 1.10000E+01, 3.00000E+01, 7.80000E+01, 1.34000E+02, 1.00000E+01, 4.10000E+01, 8.40000E+01, 5.80000E+01, 4.00000E+00, 4.10000E+01, 1.04000E+02, 1.04000E+02, 1.80000E+01, 4.90000E+01, 1.41000E+02, 1.50000E+02, 1.10000E+01, 2.90000E+01, 6.80000E+01, 7.50000E+01, 1.40000E+01, 4.90000E+01, 1.04000E+02, 1.14000E+02, 1.00000E+01, 4.60000E+01, 2.65000E+02, 3.24000E+02, 8.00000E+00, 5.50000E+01, 1.81000E+02, 2.17000E+02, 6.00000E+00, 5.00000E+01, 9.10000E+01, 1.04000E+02, 1.70000E+01, 4.30000E+01, 7.80000E+01, 1.10000E+02, 7.00000E+00, 5.20000E+01, 4.50000E+01, 8.00000E+01, 3.00000E+00, 1.20000E+01, 9.00000E+01, 8.20000E+01, 5.00000E+00, 1.40000E+01, 8.10000E+01, 7.80000E+01, 1.90000E+01, 2.80000E+01, 7.40000E+01, 5.90000E+01, 7.00000E+00, 2.40000E+01, 9.50000E+01, 8.90000E+01, 5.00000E+00, 3.20000E+01, 2.12000E+02, 2.00000E+02, 5.00000E+00, 3.60000E+01, 1.30000E+01, 1.80000E+01, 1.30000E+01, 4.90000E+01, 4.30000E+01, 4.60000E+01, 4.00000E+00, 1.70000E+01, 3.90000E+01, 4.80000E+01, 2.00000E+00, 1.30000E+01, 2.50000E+01, 2.10000E+01, 7.00000E+00, 2.70000E+01, 3.40000E+01, 3.10000E+01, 6.00000E+00, 1.50000E+01, 2.50000E+01, 3.00000E+01, 5.00000E+00, 2.50000E+01, 3.50000E+01, 2.00000E+01, 2.00000E+00, 2.00000E+01, 7.00000E+01, 6.00000E+01, 1.00000E+01, 1.50000E+01, 1.50000E+01, 2.50000E+01, 1.50000E+01, 3.00000E+01, 7.50000E+01, 1.67000E+02, 9.60000E+01, 8.30000E+01, 9.50000E+01, 1.87000E+02, 1.03000E+02, 9.30000E+01, 1.04000E+02, 1.14000E+02, 1.00000E+01, 4.60000E+01, 2.50000E+01, 1.50000E+01, 1.00000E+01, 2.50000E+01, 3.20000E+01, 9.00000E+00, 4.00000E+00, 2.40000E+01), .Dim=c(33, 4)), n_11= structure(.Data= c(9.00000E+00, 1.80000E+01, 0.00000E+00, 2.00000E+00, 6.00000E+00, 5.00000E+00, 1.00000E+00, 6.00000E+00, 2.80000E+01, 2.90000E+01, 1.00000E+00, 5.00000E+00, 1.10000E+01, 1.90000E+01, 1.00000E+00, 5.00000E+00, 9.00000E+00, 7.00000E+00, 0.00000E+00, 4.00000E+00, 2.10000E+01, 2.20000E+01, 2.00000E+00, 8.00000E+00, 1.80000E+01, 2.60000E+01, 1.00000E+00, 4.00000E+00, 1.20000E+01, 1.40000E+01, 2.00000E+00, 6.00000E+00, 1.60000E+01, 1.80000E+01, 1.00000E+00, 7.00000E+00, 3.40000E+01, 4.80000E+01, 0.00000E+00, 6.00000E+00, 2.50000E+01, 3.10000E+01, 0.00000E+00, 6.00000E+00, 1.30000E+01, 1.60000E+01, 1.00000E+00, 6.00000E+00, 1.10000E+01, 1.60000E+01, 1.00000E+00, 6.00000E+00, 6.00000E+00, 1.10000E+01, 0.00000E+00, 1.00000E+00, 1.30000E+01, 1.20000E+01, 0.00000E+00, 1.00000E+00, 1.20000E+01, 1.10000E+01, 1.00000E+00, 3.00000E+00, 1.20000E+01, 9.00000E+00, 0.00000E+00, 3.00000E+00, 2.00000E+00, 1.40000E+01, 0.00000E+00, 3.00000E+00, 3.50000E+01, 3.20000E+01, 0.00000E+00, 4.00000E+00, 1.00000E+00, 2.00000E+00, 1.00000E+00, 7.00000E+00, 6.00000E+00, 7.00000E+00, 0.00000E+00, 2.00000E+00, 5.00000E+00, 7.00000E+00, 0.00000E+00, 1.00000E+00, 3.00000E+00, 2.00000E+00, 0.00000E+00, 3.00000E+00, 4.00000E+00, 4.00000E+00, 0.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00), .Dim=c(33, 4)), p_1.1SW=c(1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 9.56700E-01, 1.00000E+00, 1.00000E+00, 1.00000E+00, 8.39700E-01, 9.36300E-01, 1.00000E+00, 1.00000E+00, 9.99700E-01, 9.99600E-01, 1.00000E+00, 9.98100E-01, 9.48700E-01, 1.00000E+00, 1.00000E+00, 9.93200E-01, 7.13100E-01, 1.00000E+00, 6.24200E-01, 9.88600E-01, 1.00000E+00, 9.49600E-01, 1.00000E+00, 1.00000E+00, 1.00000E+00, 9.80000E-01, 9.50000E-01, 9.50000E-01, 8.00000E-01), no_ech_1.1SW=c(1.57000E+02, 1.40000E+01, 7.80000E+01, 6.30000E+01, 6.30000E+01, 4.90000E+01, 8.80000E+01, 3.10000E+01, 3.60000E+01, 3.11000E+02, 2.21000E+02, 4.50000E+01, 3.30000E+01, 2.40000E+01, 3.40000E+01, 4.30000E+01, 2.60000E+01, 3.50000E+01, 1.54000E+02, 7.00000E+00, 1.30000E+01, 1.80000E+01, 7.00000E+00, 1.30000E+01, 8.00000E+00, 1.30000E+01, 5.00000E+00, 3.00000E+00, 2.00000E+00, 5.00000E+00, 5.00000E+00, 2.00000E+00, 5.00000E+00), no_ech_MSW= structure(.Data= c(2.00000E+00, 4.00000E+00, 2.00000E+00, 2.00000E+00, NA, 0.00000E+00, 2.00000E+00, 7.00000E+00, 1.00000E+00, NA, 2.00000E+00, 0.00000E+00, 7.00000E+00, 1.00000E+00, NA, 0.00000E+00, 1.00000E+00, 4.00000E+00, 5.00000E+00, NA, 0.00000E+00, 6.00000E+00, 8.00000E+00, 5.00000E+00, NA, 2.00000E+00, 9.00000E+00, 1.30000E+01, 0.00000E+00, NA, 0.00000E+00, 1.00000E+00, 6.00000E+00, 1.00000E+00, NA, 5.00000E+00, 0.00000E+00, 1.80000E+01, 1.00000E+00, NA, 0.00000E+00, 0.00000E+00, 1.20000E+01, 0.00000E+00, NA, 1.60000E+01, 1.30000E+01, 9.00000E+00, 4.00000E+00, NA, 1.00000E+00, 1.20000E+01, 1.50000E+01, 0.00000E+00, NA, 3.00000E+00, 0.00000E+00, 1.90000E+01, 3.00000E+00, NA, 0.00000E+00, 3.00000E+00, 2.00000E+00, 9.00000E+00, NA, 1.00000E+00, 2.00000E+00, 1.00000E+00, 0.00000E+00, NA, 0.00000E+00, 5.00000E+00, 1.00000E+00, 0.00000E+00, NA, 5.00000E+00, 8.00000E+00, 9.00000E+00, 3.00000E+00, NA, 0.00000E+00, 4.00000E+00, 2.00000E+00, 1.00000E+00, NA, 8.00000E+00, 0.00000E+00, 1.00000E+01, 1.00000E+00, NA, 0.00000E+00, 5.00000E+00, 7.00000E+00, 0.00000E+00, NA, 0.00000E+00, 0.00000E+00, 1.10000E+01, 0.00000E+00, NA, 1.00000E+00, 0.00000E+00, 4.00000E+00, 0.00000E+00, NA, 3.00000E+00, 0.00000E+00, 1.00000E+00, 1.00000E+00, NA, 4.00000E+00, 0.00000E+00, 5.00000E+00, 2.00000E+00, NA, 0.00000E+00, 3.00000E+00, 3.00000E+00, 2.00000E+00, NA, 0.00000E+00, 3.00000E+00, 4.00000E+00, 0.00000E+00, NA, 1.00000E+00, 0.00000E+00, 3.00000E+00, 0.00000E+00, NA, 2.00000E+00, 0.00000E+00, 2.00000E+00, 0.00000E+00, NA, 0.00000E+00, 0.00000E+00, 2.00000E+00, 5.00000E+00, NA, 2.00000E+00, 0.00000E+00, 2.00000E+00, 0.00000E+00, NA, 4.00000E+00, 0.00000E+00, 5.00000E+00, 2.00000E+00, NA, 8.00000E+00, 0.00000E+00, 1.00000E+01, 1.00000E+00, NA, 3.00000E+00, 1.00000E+00, 0.00000E+00, 0.00000E+00, NA, 0.00000E+00, 1.00000E+00, 5.00000E+00, 0.00000E+00, NA), .Dim=c(33, 5)))
This diff is collapsed.
modelCheck('/home/basp-meco88/Documents/RESEARCH/PROJECTS/ORE/Abundance/Nivelle/adult/bugs/model_adult-Nivelle.R.txt')
modelData('/home/basp-meco88/Documents/RESEARCH/PROJECTS/ORE/Abundance/Nivelle/adult/bugs/data.txt')
modelCompile(1)
modelCompile(2)
modelSetRN(1)
modelInits('/home/basp-meco88/Documents/RESEARCH/PROJECTS/ORE/Abundance/Nivelle/adult/bugs/inits1.txt',1)
modelInits('/home/basp-meco88/Documents/RESEARCH/PROJECTS/ORE/Abundance/Nivelle/adult/bugs/inits2.txt',2)
modelGenInits()
modelUpdate(500,1,500)
modelUpdate(5000,5,5000)
samplesSet(mup_11_1)
samplesSet(sigmap_11_1)
samplesSet(mup_11_2)
......@@ -139,7 +140,7 @@ summarySet(c_3SW)
summarySet(c_tot)
summarySet(P_1SW)
summarySet(P_MSW)
modelUpdate(1000,1,1000)
modelUpdate(10000,5,10000)
samplesCoda('*', '/home/basp-meco88/Documents/RESEARCH/PROJECTS/ORE/Abundance/Nivelle/adult/bugs//')
summaryStats('*')
modelQuit('y')
=============================
DIAGNOSTICS
=============================
Convergence: gelman-Rubin R test
Potential scale reduction factors:
Point est. Upper C.I.
mup_11_1 1.00 1.01
sigmap_11_1 1.00 1.01
mup_11_2 1.00 1.00
sigmap_11_2 1.02 1.06
mupi_EF 1.00 1.00
sigmapi_EF 1.00 1.00
mup_21 1.00 1.00
sigmap_21 1.00 1.01
k_1 1.00 1.00
k_2 1.00 1.01
eta_1 1.00 1.00
eta_2 1.00 1.00
rho 1.00 1.00
shape_lambda 1.00 1.00
rate_lambda 1.00 1.00
lambda_tot0 1.00 1.01
a_1.1SW 1.00 1.00
a_2.1SW 1.00 1.00
Multivariate psrf
1.01
---------------------------
Heidelberger and Welch's convergence diagnostic
......@@ -9,48 +35,94 @@ heidel.diag is a run length control diagnostic based on a criterion of relative
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
mup_11_1 passed 1 4.04e-01
sigmap_11_1 failed NA 2.25e-06
mup_11_2 passed 401 1.18e-01
sigmap_11_2 passed 1 9.81e-01
mupi_EF passed 101 1.96e-01
sigmapi_EF passed 1 5.78e-01
mup_21 passed 101 3.24e-01
sigmap_21 passed 1 8.06e-01
k_1 passed 1 5.27e-01
k_2 passed 1 4.58e-01
eta_1 passed 1 4.09e-01
eta_2 passed 1 1.45e-01
rho passed 1 5.06e-01
shape_lambda passed 1 6.97e-01
rate_lambda passed 1 5.58e-01
lambda_tot0 passed 1 6.10e-01
a_1.1SW passed 1 3.33e-01
a_2.1SW passed 1 3.16e-01
[[1]]
Stationarity start p-value
test iteration
mup_11_1 passed 1 0.2470
sigmap_11_1 passed 1 0.7451
mup_11_2 passed 1 0.6874
sigmap_11_2 passed 1 0.1751
mupi_EF passed 1 0.5407
sigmapi_EF passed 1 0.9894
mup_21 passed 1 0.5731
sigmap_21 passed 1 0.6131
k_1 passed 1 0.4446
k_2 passed 1 0.3467
eta_1 passed 1 0.3852
eta_2 passed 1 0.9147
rho passed 1 0.3876
shape_lambda passed 1 0.8234
rate_lambda passed 1 0.8821
lambda_tot0 passed 1 0.6623
a_1.1SW passed 1 0.1238
a_2.1SW passed 1 0.0793
Halfwidth Mean Halfwidth
test
mup_11_1 passed 0.241 0.002201
sigmap_11_1 passed 0.486 0.020994
mup_11_2 passed 0.123 0.001182
sigmap_11_2 failed 0.289 0.029270
mupi_EF passed 0.248 0.001220
sigmapi_EF passed 0.754 0.008462
mup_21 passed 0.652 0.000880
sigmap_21 passed 0.567 0.006365
k_1 passed 1.169 0.020959
k_2 passed 2.425 0.039284
eta_1 passed 3.408 0.030719
eta_2 passed 5.197 0.040380
rho passed 0.944 0.002132
shape_lambda passed 2.842 0.021423
rate_lambda passed 0.015 0.000125
lambda_tot0 passed 140.771 1.680774
a_1.1SW passed 10.305 0.084024
a_2.1SW passed 3.029 0.023497
[[2]]
Stationarity start p-value
test iteration
mup_11_1 passed 1 0.7690
sigmap_11_1 passed 1 0.5792
mup_11_2 passed 1 0.4289
sigmap_11_2 passed 1 0.2448
mupi_EF passed 1 0.0806
sigmapi_EF passed 1 0.7410
mup_21 passed 1 0.3183
sigmap_21 passed 1 0.3247
k_1 passed 1 0.1192
k_2 passed 1 0.1103
eta_1 passed 1 0.6892
eta_2 passed 1 0.2011
rho passed 1 0.1269
shape_lambda passed 1 0.1239
rate_lambda passed 1 0.1399
lambda_tot0 passed 1 0.2453
a_1.1SW passed 1 0.3396
a_2.1SW passed 1 0.1362
Halfwidth Mean Halfwidth
test
mup_11_1 passed 0.2437 0.013806
sigmap_11_1 <NA> NA NA
mup_11_2 passed 0.1283 0.004165
sigmap_11_2 failed 0.3165 0.083794
mupi_EF passed 0.2514 0.005447
sigmapi_EF passed 0.7237 0.043936
mup_21 passed 0.6517 0.005495
sigmap_21 passed 0.5410 0.041786
k_1 passed 1.1453 0.097843
k_2 passed 2.3712 0.151751
eta_1 passed 3.3987 0.118404
eta_2 passed 5.3360 0.181535
rho passed 0.9475 0.011612
shape_lambda passed 2.7168 0.174404
rate_lambda passed 0.0143 0.000985
lambda_tot0 passed 138.7444 8.752523
a_1.1SW passed 10.1821 0.526254
a_2.1SW passed 2.9996 0.150857
mup_11_1 passed 0.2435 0.002453
sigmap_11_1 passed 0.4806 0.022835
mup_11_2 passed 0.1231 0.000961
sigmap_11_2 passed 0.3167 0.017793
mupi_EF passed 0.2472 0.001151
sigmapi_EF passed 0.7548 0.008312
mup_21 passed 0.6517 0.000908
sigmap_21 passed 0.5603 0.006236
k_1 passed 1.1637 0.018626
k_2 passed 2.4293 0.037790
eta_1 passed 3.4121 0.027121
eta_2 passed 5.1754 0.042218
rho passed 0.9450 0.002336
shape_lambda passed 2.8507 0.022661
rate_lambda passed 0.0151 0.000137
lambda_tot0 passed 139.6302 1.621970
a_1.1SW passed 10.2410 0.081807
a_2.1SW passed 3.0105 0.022534
---------------------------
Geweke's convergence diagnostic
......@@ -61,23 +133,88 @@ The test statistic is a standard Z-score: the difference between the two sample
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.
[[1]]
Fraction in 1st window = 0.1
Fraction in 2nd window = 0.5
mup_11_1 sigmap_11_1 mup_11_2 sigmap_11_2 mupi_EF sigmapi_EF
-1.259 -0.558 -0.129 0.366 -0.816 0.348
mup_21 sigmap_21 k_1 k_2 eta_1 eta_2
0.182 0.163 -0.445 -0.464 0.442 1.062
rho shape_lambda rate_lambda lambda_tot0 a_1.1SW a_2.1SW
0.465 -1.091 -0.999 -0.580 -1.196 -1.534
[[2]]
Fraction in 1st window = 0.1
Fraction in 2nd window = 0.5
mup_11_1 sigmap_11_1 mup_11_2 sigmap_11_2 mupi_EF sigmapi_EF
-0.3224 -0.4088 -0.8603 -1.3132 1.5172 0.5646
1.539 0.175 1.876 -1.492 0.265 -1.001
mup_21 sigmap_21 k_1 k_2 eta_1 eta_2
-2.2853 -1.0174 -1.6749 -0.6190 1.0087 -2.4604
-0.166 1.130 -1.569 -1.256 -0.301 -0.260
rho shape_lambda rate_lambda lambda_tot0 a_1.1SW a_2.1SW
1.2729 -0.4802 -0.1552 0.0936 0.8232 0.6989
0.812 -1.653 -2.003 0.860 -0.814 -0.690
---------------------------
Raftery and Lewis's diagnostic
[[1]]
Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95
Burn-in Total Lower bound Dependence
(M) (N) (Nmin) factor (I)
mup_11_1 40 46360 3746 12.40
sigmap_11_1 200 212400 3746 56.70
mup_11_2 40 44790 3746 12.00
sigmap_11_2 2730 3035435 3746 810.00
mupi_EF 10 19235 3746 5.13
sigmapi_EF 20 24625 3746 6.57
mup_21 15 20770 3746 5.54
sigmap_21 80 84760 3746 22.60
k_1 90 103350 3746 27.60
k_2 90 100760 3746 26.90
eta_1 10 19325 3746 5.16
eta_2 10 18705 3746 4.99
rho 120 136350 3746 36.40
shape_lambda 20 24160 3746 6.45
rate_lambda 15 22050 3746 5.89
lambda_tot0 15 21870 3746 5.84
a_1.1SW 15 22050 3746 5.89
a_2.1SW 15 22665 3746 6.05
[[2]]
Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95
Burn-in Total Lower bound Dependence
(M) (N) (Nmin) factor (I)
mup_11_1 30 49520 3746 13.20
sigmap_11_1 340 370580 3746 98.90
mup_11_2 40 55430 3746 14.80
sigmap_11_2 225 259995 3746 69.40
mupi_EF 10 18680 3746 4.99
sigmapi_EF 25 27700 3746 7.39
mup_21 15 20370 3746 5.44
sigmap_21 70 73730 3746 19.70
k_1 210 239085 3746 63.80
k_2 125 140950 3746 37.60
eta_1 10 19615 3746 5.24
eta_2 15 20305 3746 5.42
rho 150 158410 3746 42.30
shape_lambda 20 25675 3746 6.85
rate_lambda 20 24770 3746 6.61
lambda_tot0 15 21510 3746 5.74
a_1.1SW 15 21160 3746 5.65
a_2.1SW 15 21225 3746 5.67
You need a sample size of at least 3746 with these values of q, r and s
......@@ -44,7 +44,8 @@ if(!file.exists(paste("inits/init-",site,"-",stade,year,".txt",sep=""))){
#load(paste('inits/inits_',stade,'.Rdata',sep="")) # chargement des inits
#if(site == "Bresle" && stade == "adult") {inits <- list(read.bugsdata(paste("inits/init-",site,"-",stade,year,".txt",sep="")))}
#if(site == "Nivelle") {inits <- list(read.bugsdata(paste("inits/init-",site,"-",stade,year,".txt",sep="")))}
inits <- list(read.bugsdata(paste("inits/init-",site,"-",stade,year,".txt",sep="")))
inits.tmp <- read.bugsdata(paste("inits/init-",site,"-",stade,year,".txt",sep=""))
inits <- rep(list(inits.tmp),2)
#------------------------MODEL----------------------------------##
model <- paste("model/model_",stade,"-",site,".R",sep="") # path of the model
......@@ -56,11 +57,11 @@ filename <- file.path(work.dir, model)
#---------------------------ANALYSIS-----------------------------##
nChains = length(inits) # Number of chains to run.
nChains = 2 #length(inits) # Number of chains to run.
adaptSteps = 1000 # Number of steps to "tune" the samplers.
nburnin=5000 # Number of steps to "burn-in" the samplers.
nstore=25000 # Total number of steps in chains to save.
nthin=2 # Number of steps to "thin" (1=keep every step).
nstore=10000 # Total number of steps in chains to save.
nthin=5 # Number of steps to "thin" (1=keep every step).
#nPerChain = ceiling( ( numSavedSteps * thinSteps ) / nChains ) # Steps per chain.
### Start of the run ###
......@@ -139,7 +140,7 @@ cat("=============================\n")
if (nChains > 1) {
cat("Convergence: gelman-Rubin R test\n")
gelman.diag(fit.mcmc[,which(varnames(fit.mcmc)%in%parameterstotest)])
gelman.diag(fit.mcmc[,which(varnames(fit.mcmc)%in%parameterstotest)],multivariate=TRUE)
}
......@@ -177,6 +178,7 @@ pdf(paste('results/Results_',stade,"_",year,'.pdf',sep=""))
traplot(fit.mcmc[,which(varnames(fit.mcmc)%in%parameterstotest)])
# caterplot(fit.mcmc,parameters[i])
#}
gelman.plot(fit.mcmc[,which(varnames(fit.mcmc)%in%parameterstotest)])
dev.off()
......
****** Sorry something went wrong in procedure Mark in module Kernel ******
This diff is collapsed.
This diff is collapsed.
modelCheck('/home/basp-meco88/Documents/RESEARCH/PROJECTS/ORE/Abundance/Nivelle/tacon/bugs/model_tacon-Nivelle.R.txt')
modelData('/home/basp-meco88/Documents/RESEARCH/PROJECTS/ORE/Abundance/Nivelle/tacon/bugs/data.txt')
modelCompile(1)
modelCompile(2)
modelSetRN(1)
modelInits('/home/basp-meco88/Documents/RESEARCH/PROJECTS/ORE/Abundance/Nivelle/tacon/bugs/inits1.txt',1)
modelInits('/home/basp-meco88/Documents/RESEARCH/PROJECTS/ORE/Abundance/Nivelle/tacon/bugs/inits2.txt',2)
modelGenInits()
modelUpdate(1000,1,1000)
modelUpdate(5000,5,5000)
samplesSet(mu_p_srem)
samplesSet(sd_logit_p_srem)
samplesSet(epsilon_p)
......@@ -67,7 +68,7 @@ summarySet(jres_ns_runs)
summarySet(dj)
summarySet(dj_nat)
summarySet(dj_comp)
modelUpdate(2000,1,2000)
modelUpdate(10000,5,10000)
samplesCoda('*', '/home/basp-meco88/Documents/RESEARCH/PROJECTS/ORE/Abundance/Nivelle/tacon/bugs//')
summaryStats('*')
modelQuit('y')
> cat("=============================\n")
=============================
> cat("DIAGNOSTICS\n")
DIAGNOSTICS
> cat("=============================\n")
=============================
> if (nChains > 1) {
+ cat("Convergence: gelman-Rubin R test\n")
+ gelman.diag(fit.mcmc[,which(varnames(fit.mcmc)%in%parameterstotest)])
+ }
> cat("\n---------------------------\n")
Convergence: gelman-Rubin R test
Potential scale reduction factors:
Point est. Upper C.I.
mu_p_srem 1.00 1.01
sd_logit_p_srem 1.00 1.00
epsilon_p 1.00 1.02
pi_dj 1.00 1.00
zeta_alpha_dj 1.00 1.00
mu_dj_nat 1.02 1.10
k_cpue 1.00 1.02
eta_cpue 1.02 1.03
rho_s 1.01 1.03
sd_s_rec 1.00 1.00
p_cpue 1.00 1.01
Multivariate psrf
1.02
---------------------------
> cat("Heidelberger and Welch's convergence diagnostic\n")
Heidelberger and Welch's convergence diagnostic
> cat("
+ heidel.diag is a run length control diagnostic based on a criterion of relative accuracy for the estimate of the mean. The default setting c ..." ... [TRUNCATED]
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.
> heidel.diag(fit.mcmc[,which(varnames(fit.mcmc)%in%parameterstotest)], eps=0.1, pvalue=0.05)
[[1]]
Stationarity start p-value
test iteration
mu_p_srem passed 1 0.698
sd_logit_p_srem passed 1 0.201
epsilon_p passed 1 0.540
pi_dj passed 1 0.157
zeta_alpha_dj passed 1 0.133
mu_dj_nat passed 1 0.366
k_cpue passed 1 0.488
eta_cpue passed 1 0.728
rho_s passed 1 0.163
sd_s_rec passed 1 0.258
p_cpue passed 1 0.946
Halfwidth Mean Halfwidth
test
mu_p_srem passed 0.516 0.00542
sd_logit_p_srem passed 1.091 0.00951
epsilon_p passed 0.561 0.00855
pi_dj passed 0.726 0.00418
zeta_alpha_dj passed 4.022 0.06083
mu_dj_nat passed 0.128 0.00198
k_cpue passed 124.367 0.90915
eta_cpue passed 2.259 0.18766
rho_s passed 0.162 0.00595
sd_s_rec passed 0.444 0.00377
p_cpue passed 0.631 0.01055
[[2]]
Stationarity start p-value
test iteration
mu_p_srem passed 601 0.0669
sd_logit_p_srem passed 601 0.2366
epsilon_p passed 801 0.4381
pi_dj passed 601 0.1008
zeta_alpha_dj passed 1 0.1208
mu_dj_nat passed 1 0.0641
k_cpue passed 801 0.1195
eta_cpue passed 1 0.5970
rho_s passed 1 0.3307
sd_s_rec passed 1 0.3200
p_cpue passed 1 0.5303
mu_p_srem passed 1 0.732
sd_logit_p_srem passed 1 0.783
epsilon_p passed 1 0.527
pi_dj passed 1 0.456
zeta_alpha_dj passed 1 0.903
mu_dj_nat passed 1 0.444
k_cpue passed 1 0.410
eta_cpue passed 2001 0.118
rho_s passed 1 0.368
sd_s_rec passed 1 0.527
p_cpue passed 1 0.354
Halfwidth Mean Halfwidth
test
mu_p_srem passed 0.544 0.01652
sd_logit_p_srem passed 1.024 0.03291
epsilon_p passed 0.504 0.01790
pi_dj passed 0.731 0.01765
zeta_alpha_dj passed 3.906 0.10157
mu_dj_nat passed 0.129 0.00664
k_cpue passed 130.697 3.53667
eta_cpue failed 2.127 0.79120
rho_s failed 0.177 0.01827
sd_s_rec passed 0.442 0.01248
p_cpue passed 0.628 0.04318
> cat("\n---------------------------\n")
mu_p_srem passed 0.516 0.00531
sd_logit_p_srem passed 1.096 0.00882
epsilon_p passed 0.560 0.00833
pi_dj passed 0.724 0.00451
zeta_alpha_dj passed 4.100 0.06321
mu_dj_nat passed 0.130 0.00194
k_cpue passed 124.085 0.92403
eta_cpue passed 2.274 0.19046
rho_s passed 0.152 0.00602
sd_s_rec passed 0.444 0.00347
p_cpue passed 0.637 0.01037
---------------------------
> cat("Geweke's convergence diagnostic\n")
---------------------------
Geweke's convergence diagnostic
> cat("
+ 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 ..." ... [TRUNCATED]
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.
[[1]]
Fraction in 1st window = 0.1
Fraction in 2nd window = 0.5
> geweke.diag(fit.mcmc[,which(varnames(fit.mcmc)%in%parameterstotest)], frac1 = 0.1, frac2 = 0.5)
mu_p_srem sd_logit_p_srem epsilon_p pi_dj zeta_alpha_dj
-0.8963 -1.0340 1.6530 -0.8663 -0.2139
mu_dj_nat k_cpue eta_cpue rho_s sd_s_rec
0.9978 -1.1527 -0.0156 -0.7487 1.2660
p_cpue
0.5001
[[2]]
Fraction in 1st window = 0.1
Fraction in 2nd window = 0.5
mu_p_srem sd_logit_p_srem epsilon_p pi_dj zeta_alpha_dj mu_dj_nat k_cpue eta_cpue
-6.641 7.805 6.139 1.559 -5.074 0.192 -5.779 -1.319
rho_s sd_s_rec p_cpue
-0.594 0.359 -1.285
mu_p_srem sd_logit_p_srem epsilon_p pi_dj zeta_alpha_dj
0.823 -0.791 -0.658 -1.549 -0.150
mu_dj_nat k_cpue eta_cpue rho_s sd_s_rec
-0.814 0.444 -1.350 0.307 1.182
p_cpue
-1.539
> cat("\n---------------------------\n")
---------------------------
> cat("Raftery and Lewis's diagnostic\n")
Raftery and Lewis's diagnostic
[[1]]
> raftery.diag(fit.mcmc[,which(varnames(fit.mcmc)%in%parameterstotest)], q=0.025, r=0.005, s=0.95, converge.eps=0.001)
Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95
Burn-in Total Lower bound Dependence
(M) (N) (Nmin) factor (I)
mu_p_srem 120 152320 3746 40.70
sd_logit_p_srem 80 106440 3746 28.40
epsilon_p 150 164200 3746 43.80
pi_dj 25 29910 3746 7.98
zeta_alpha_dj 15 21510 3746 5.74
mu_dj_nat 225 232775 3746 62.10
k_cpue 120 143550 3746 38.30
eta_cpue 60 75220 3746 20.10
rho_s 80 100200 3746 26.70
sd_s_rec 50 55950 3746 14.90
p_cpue 60 75220 3746 20.10
[[2]]
Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95
Burn-in Total Lower bou