#------------------------------------------------------------------------------- # R-Code for the simulation studies to assess effects of violations of model assumptions # for the model given in # Korner-Nievergelt, F., Liechti, F., Hahn, S. (2012) Migratory connectivity derived from sparse ring reencounter data with unknown numbers of ringed birds. Journal of Ornithology # DOI 10.1007/s10336-011-0793-z #------------------------------------------------------------------------------- # 6 different simulation studies to assess effects of violations of model assumptions #----------------------------------------------------------------------------------- # 1) perfect data # 2) heterogeneous p # 3) heterogeneous phi # 4) positive correlation between r and latitute of breeding area # 5) negative correlation between r and latitute of breeding area # 6) heterogeneous r #--------------------------------------------------------------------------------- # load functions and libraries pfadname<-"...." setwd(pfadname) source("..../simul.ch.r") source("..../marray.recapture.r") source("..../shapeparameter.r") upper <- function(x) quantile(x, prob=0.975) lower <- function(x) quantile(x, prob=0.025) library(BRugs) library(R2WinBUGS) ################################################################################ # model model<-function(){ # priors für N-modell S~dunif(0,1) # return probability (phi) for(j in 1:npop){ p[j]~dunif(0,1) # recapture probability in breeding area } # multinomial likelihood füfor CJS-Modell for(t in 1:(J-1)){ MArecaptCE[t,1:J]~dmulti(prCE[t,], rCE[t]) MArecaptEE[t,1:J]~dmulti(prEE[t,], rEE[t]) MArecaptSE[t,1:J]~dmulti(prSE[t,], rSE[t]) MArecaptNEE[t,1:J]~dmulti(prNEE[t,], rNEE[t]) } # number of released birds for(t in 1:(J-1)){ rCE[t] <- sum(MArecaptCE[t,]) rEE[t] <- sum(MArecaptEE[t,]) rSE[t] <- sum(MArecaptSE[t,]) rNEE[t] <- sum(MArecaptNEE[t,]) } # cell probabilities for(j in 1:npop){ q[j] <- 1-p[j] } for(t in 1:(J-1)){ # main diagonal prCE[t,t] <- S*p[1] prEE[t,t] <- S*p[2] prSE[t,t] <- S*p[3] prNEE[t,t]<- S*p[4] for(k in (t+1):(J-1)){ # above main diagonal prCE[t,k] <- pow(S, (k-t+1))*pow(q[1], (k-t))*p[1] prEE[t,k] <- pow(S, (k-t+1))*pow(q[2], (k-t))*p[2] prSE[t,k] <- pow(S, (k-t+1))*pow(q[3], (k-t))*p[3] prNEE[t,k] <- pow(S, (k-t+1))*pow(q[4], (k-t))*p[4] } } for(t in 2:(J-1)){ # belwo main diagonal for(k in 1:(t-1)){ prCE[t,k] <- 0 prEE[t,k] <- 0 prSE[t,k] <- 0 prNEE[t,k] <- 0 } } # last column: probability of non recapture for(t in 1:(J-1)){ prCE[t,J] <- 1-sum(prCE[t,1:(J-1)]) prEE[t,J] <- 1-sum(prEE[t,1:(J-1)]) prSE[t,J] <- 1-sum(prSE[t,1:(J-1)]) prNEE[t,J] <- 1-sum(prNEE[t,1:(J-1)]) } # number of ringed in the population for(j in 1:npop){ N[j] <- n[j]/(1-(1-S)/(1-S*(1-p[j]))) } # migration rate model for(j in 1:npop){ for(a in 1:3){ r.korr[j,a]<-r[a]*m[j,a] muW[j,a] <- r.korr[j,a] * N[j] wirec[j,a]~dpois(muW[j,a]) } } # priors für migration rate model for(a in 1:3){ r[a]~dunif(0,1) # re-encounter probabilities } for(j in 1:npop){ m[j,1]~dunif(0,1) mquer[j] <- 1-m[j,1] m[j,2]~dunif(0, mquer[j] ) m[j,3] <- max(1-sum(m[j,1:2]), 0) } } writeModel(model, "model.txt") #------------------------------------------------------------------------------- # simulate perfect data # set the true parameter values N <- rep(50000, 4) p <- c(0.15, 0.2, 0.1, 0.1) S <- 0.5 m <- read.table("mestnoninfprior100920.txt", header=TRUE) mrate <- matrix(m$mean, ncol=3, nrow=4, byrow=TRUE) r <- rep(c(0.001, 0.01, 0.003), 4) J <- 8 # n occations nsim <- 20 biasmmean <- matrix(ncol=12, nrow=nsim) biasSmean <- numeric(nsim) biaspmean <- matrix(ncol=4, nrow=nsim) biasrmean <- matrix(ncol=3, nrow=nsim) biasNmean <- matrix(ncol=4, nrow=nsim) biasmmedian <- matrix(ncol=12, nrow=nsim) biasSmedian <- numeric(nsim) biaspmedian <- matrix(ncol=4, nrow=nsim) biasrmedian <- matrix(ncol=3, nrow=nsim) biasNmedian <- matrix(ncol=4, nrow=nsim) for(sim in 1:nsim){ # simulate breeding data pmin1wf <- 1-(1-S)/(1-S*(1-p)) n <- rbinom(4, N, pmin1wf) chPop1 <- simul.ch(S, p[1], c(n[1], rep(0, J-1)), J) chPop2 <- simul.ch(S, p[2], c(n[2], rep(0, J-1)), J) chPop3 <- simul.ch(S, p[3], c(n[3], rep(0, J-1)), J) chPop4 <- simul.ch(S, p[4], c(n[4], rep(0, J-1)), J) MAPop1 <- marray.recapture(chPop1)$marray MAPop2 <- marray.recapture(chPop2)$marray MAPop3 <- marray.recapture(chPop3)$marray MAPop4 <- marray.recapture(chPop4)$marray # simulate non-breeding data wf <- rbinom(3*4, size=rep(N, 3), prob=m$mean*r) wirec <- matrix(wf, ncol=3, nrow=4, byrow=TRUE) datax<-list(MArecaptCE=MAPop1, MArecaptEE=MAPop2, MArecaptSE=MAPop3, MArecaptNEE=MAPop4,n=n, J=J, wirec=wirec, npop=npop) bugsData(datax, fileName = "data.txt") minit<-wirec/apply(wirec,1,sum) minit <- matrix(c(0.61000, 0.16570, 0.22440, 0.02289, 0.30590, 0.67120, 0.01112, 0.64550, 0.34340, 0.04396, 0.08658, 0.86950), ncol=3, nrow=4, byrow=TRUE) minit[,3]<-1-apply(minit[,c(1:2)], 1, sum) minit <- minit[,c(2,3,1)] minit[,3] <- NA init1 <- list(p=runif(npop, 0.3,0.7), m=minit, r=runif(3, 0,0.4)) init2 <- list(p=runif(npop, 0.3,0.7), m=minit, r=runif(3, 0,0.4)) inits<-list(init1, init2) bugsInits(inits, numChains = 2, c("init1.txt", "init2.txt"), digits = 5) mod <- bugs(datax, inits=list(init1, init2), parameters.to.save=c("p", "S", "N", "m", "r") , model.file="model.txt", n.chains=2, n.iter=250000, n.burnin=200000, n.sims = 1000, debug=FALSE, DIC=TRUE, digits=5, codaPkg=FALSE, bugs.directory="c:/Programme/WinBUGS14/", program="WinBUGS", working.directory=getwd()) biaspmean[sim, ] <- mod$summary[1:4,1]-p biaspmedian[sim, ] <- mod$summary[1:4,5]-p biasSmean[sim] <- mod$summary[5,1]-S biasSmedian[sim] <- mod$summary[5,5]-S biasNmean[sim, ] <- mod$summary[6:9,1]-N biasNmedian[sim, ] <- mod$summary[6:9,5]-N biasmmean[sim, ] <- mod$summary[10:21,1]-m$mean biasmmedian[sim, ] <- mod$summary[10:21,5]-m$mean biasrmean[sim, ] <- mod$summary[22:24,1]-r[1:3] biasrmedian[sim, ] <- mod$summary[22:24,5]-r[1:3] # save results bias <- list(biaspmean=biaspmean, biaspmedian=biaspmedian, biasSmean=biasSmean, biasSmedian=biasSmedian, biasNmean=biasNmean, biasNmedian=biasNmedian, biasmmean=biasmmean,biasmmedian=biasmmedian, biasrmean=biasrmean, biasrmedian=biasrmedian) dput(bias, "bias_perfectdata.txt") } # close sim #------------------------------------------------------------------------------ # simulate data with heterogeneity in recapture probability p between individuals # set the true parameter values pa <- 2; pb <- 9 biasmmean <- matrix(ncol=12, nrow=nsim) biasSmean <- numeric(nsim) biaspmean <- matrix(ncol=4, nrow=nsim) biasrmean <- matrix(ncol=3, nrow=nsim) biasNmean <- matrix(ncol=4, nrow=nsim) biasmmedian <- matrix(ncol=12, nrow=nsim) biasSmedian <- numeric(nsim) biaspmedian <- matrix(ncol=4, nrow=nsim) biasrmedian <- matrix(ncol=3, nrow=nsim) biasNmedian <- matrix(ncol=4, nrow=nsim) for(sim in 1:nsim){ # simulate breeding data pm <- pa/(pa+pb) # mean of beta distribution pmin1wf <- 1-(1-S)/(1-S*(1-p)) # Wahrscheinlichkeit für mind 1 Wiederfang n <- rbinom(4, N, pmin1wf) p1 <- rbeta(n[1], pa, pb) p2 <- rbeta(n[2], pa, pb) p3 <- rbeta(n[3], pa, pb) p4 <- rbeta(n[4], pa, pb) chPop1 <- simul.ch(S, p1, c(n[1], rep(0, J-1)), J, type="S.p(i)") chPop2 <- simul.ch(S, p2, c(n[2], rep(0, J-1)), J, type="S.p(i)") chPop3 <- simul.ch(S, p3, c(n[3], rep(0, J-1)), J, type="S.p(i)") chPop4 <- simul.ch(S, p4, c(n[4], rep(0, J-1)), J, type="S.p(i)") MAPop1 <- marray.recapture(chPop1)$marray MAPop2 <- marray.recapture(chPop2)$marray MAPop3 <- marray.recapture(chPop3)$marray MAPop4 <- marray.recapture(chPop4)$marray # simulate non-breeding data wf <- rbinom(3*4, size=rep(N, 3), prob=m$mean*r) wirec <- matrix(wf, ncol=3, nrow=4, byrow=TRUE) datax<-list(MArecaptCE=MAPop1, MArecaptEE=MAPop2, MArecaptSE=MAPop3, MArecaptNEE=MAPop4,n=n, J=J, wirec=wirec, npop=npop) bugsData(datax, fileName = "data.txt") minit <- matrix(c(0.61000, 0.16570, 0.22440, 0.02289, 0.30590, 0.67120, 0.01112, 0.64550, 0.34340, 0.04396, 0.08658, 0.86950), ncol=3, nrow=4, byrow=TRUE) minit[,3]<-1-apply(minit[,c(1:2)], 1, sum) minit <- minit[,c(2,3,1)] minit[,3] <- NA init1 <- list(p=runif(npop, 0.3,0.7), m=minit, r=runif(3, 0,0.4)) init2 <- list(p=runif(npop, 0.3,0.7), m=minit, r=runif(3, 0,0.4)) inits<-list(init1, init2) bugsInits(inits, numChains = 2, c("init1.txt", "init2.txt"), digits = 5) mod <- bugs(datax, inits=list(init1, init2), parameters.to.save=c("p", "S", "N", "m", "r") , model.file="model.txt", n.chains=2, n.iter=250000, n.burnin=200000, n.sims = 1000, debug=FALSE, DIC=TRUE, digits=5, codaPkg=FALSE, bugs.directory="c:/Programme/WinBUGS14/", program="WinBUGS", working.directory=getwd()) biaspmean[sim, ] <- mod$summary[1:4,1]-pm biaspmedian[sim, ] <- mod$summary[1:4,5]-pm biasSmean[sim] <- mod$summary[5,1]-S biasSmedian[sim] <- mod$summary[5,5]-S biasNmean[sim, ] <- mod$summary[6:9,1]-N biasNmedian[sim, ] <- mod$summary[6:9,5]-N biasmmean[sim, ] <- mod$summary[10:21,1]-m$mean biasmmedian[sim, ] <- mod$summary[10:21,5]-m$mean biasrmean[sim, ] <- mod$summary[22:24,1]-r[1:3] biasrmedian[sim, ] <- mod$summary[22:24,5]-r[1:3] # save results bias <- list(biaspmean=biaspmean, biaspmedian=biaspmedian, biasSmean=biasSmean, biasSmedian=biasSmedian, biasNmean=biasNmean, biasNmedian=biasNmedian, biasmmean=biasmmean,biasmmedian=biasmmedian, biasrmean=biasrmean, biasrmedian=biasrmedian) dput(bias, "bias_p_i.txt") } # close sim #------------------------------------------------------------------------------ # simulate data with heterogeneity in survival probability S between individuals # set the true parameter values Sa <- 3; Sb <- 3 biasmmean <- matrix(ncol=12, nrow=nsim) biasSmean <- numeric(nsim) biaspmean <- matrix(ncol=4, nrow=nsim) biasrmean <- matrix(ncol=3, nrow=nsim) biasNmean <- matrix(ncol=4, nrow=nsim) biasmmedian <- matrix(ncol=12, nrow=nsim) biasSmedian <- numeric(nsim) biaspmedian <- matrix(ncol=4, nrow=nsim) biasrmedian <- matrix(ncol=3, nrow=nsim) biasNmedian <- matrix(ncol=4, nrow=nsim) for(sim in 1:nsim){ # simulate breeding data Sm <- Sa/(Sa+Sb) # mean Survival probability pmin1wf <- 1-(1-Sm)/(1-Sm*(1-p)) # Wahrscheinlichkeit für mind 1 Wiederfang n <- rbinom(4, N, pmin1wf) S1 <- rbeta(n[1], Sa, Sb) S2 <- rbeta(n[2], Sa, Sb) S3 <- rbeta(n[3], Sa, Sb) S4 <- rbeta(n[4], Sa, Sb) chPop1 <- simul.ch(S1, p[1], c(n[1], rep(0, J-1)), J, type="S(i).p") chPop2 <- simul.ch(S2, p[2], c(n[2], rep(0, J-1)), J, type="S(i).p") chPop3 <- simul.ch(S3, p[3], c(n[3], rep(0, J-1)), J, type="S(i).p") chPop4 <- simul.ch(S4, p[4], c(n[4], rep(0, J-1)), J, type="S(i).p") MAPop1 <- marray.recapture(chPop1)$marray MAPop2 <- marray.recapture(chPop2)$marray MAPop3 <- marray.recapture(chPop3)$marray MAPop4 <- marray.recapture(chPop4)$marray # simulate non-breeding data wf <- rbinom(3*4, size=rep(N, 3), prob=m$mean*r) wirec <- matrix(wf, ncol=3, nrow=4, byrow=TRUE) datax<-list(MArecaptCE=MAPop1, MArecaptEE=MAPop2, MArecaptSE=MAPop3, MArecaptNEE=MAPop4,n=n, J=J, wirec=wirec, npop=npop) bugsData(datax, fileName = "data.txt") minit <- matrix(c(0.61000, 0.16570, 0.22440, 0.02289, 0.30590, 0.67120, 0.01112, 0.64550, 0.34340, 0.04396, 0.08658, 0.86950), ncol=3, nrow=4, byrow=TRUE) minit[,3]<-1-apply(minit[,c(1:2)], 1, sum) minit <- minit[,c(2,3,1)] minit[,3] <- NA init1 <- list(p=runif(npop, 0.3,0.7), m=minit, r=runif(3, 0,0.4)) init2 <- list(p=runif(npop, 0.3,0.7), m=minit, r=runif(3, 0,0.4)) inits<-list(init1, init2) bugsInits(inits, numChains = 2, c("init1.txt", "init2.txt"), digits = 5) mod <- bugs(datax, inits=list(init1, init2), parameters.to.save=c("p", "S", "N", "m", "r") , model.file="model.txt", n.chains=2, n.iter=250000, n.burnin=200000, n.sims = 1000, debug=FALSE, DIC=TRUE, digits=5, codaPkg=FALSE, bugs.directory="c:/Programme/WinBUGS14/", program="WinBUGS", working.directory=getwd()) biaspmean[sim, ] <- mod$summary[1:4,1]-p biaspmedian[sim, ] <- mod$summary[1:4,5]-p biasSmean[sim] <- mod$summary[5,1]-(Sa/(Sa+Sb)) biasSmedian[sim] <- mod$summary[5,5]-(Sa/(Sa+Sb)) biasNmean[sim, ] <- mod$summary[6:9,1]-N biasNmedian[sim, ] <- mod$summary[6:9,5]-N biasmmean[sim, ] <- mod$summary[10:21,1]-m$mean biasmmedian[sim, ] <- mod$summary[10:21,5]-m$mean biasrmean[sim, ] <- mod$summary[22:24,1]-r[1:3] biasrmedian[sim, ] <- mod$summary[22:24,5]-r[1:3] # save results bias <- list(biaspmean=biaspmean, biaspmedian=biaspmedian, biasSmean=biasSmean, biasSmedian=biasSmedian, biasNmean=biasNmean, biasNmedian=biasNmedian, biasmmean=biasmmean,biasmmedian=biasmmedian, biasrmean=biasrmean, biasrmedian=biasrmedian) dput(bias, "bias_S_i.txt") } # close sim #------------------------------------------------------------------------------ # simulate data with negative correlation of recovery probability with latitude of breeding area # set the true parameter values biasmmean <- matrix(ncol=12, nrow=nsim) biasSmean <- numeric(nsim) biaspmean <- matrix(ncol=4, nrow=nsim) biasrmean <- matrix(ncol=3, nrow=nsim) biasNmean <- matrix(ncol=4, nrow=nsim) biasmmedian <- matrix(ncol=12, nrow=nsim) biasSmedian <- numeric(nsim) biaspmedian <- matrix(ncol=4, nrow=nsim) biasrmedian <- matrix(ncol=3, nrow=nsim) biasNmedian <- matrix(ncol=4, nrow=nsim) for(sim in 1:nsim){ # simulate breeding data pmin1wf <- 1-(1-S)/(1-S*(1-p)) # Wahrscheinlichkeit für mind 1 Wiederfang n <- rbinom(4, N, pmin1wf) chPop1 <- simul.ch(S, p[1], c(n[1], rep(0, J-1)), J, type="S.p") chPop2 <- simul.ch(S, p[2], c(n[2], rep(0, J-1)), J, type="S.p") chPop3 <- simul.ch(S, p[3], c(n[3], rep(0, J-1)), J, type="S.p") chPop4 <- simul.ch(S, p[4], c(n[4], rep(0, J-1)), J, type="S.p") MAPop1 <- marray.recapture(chPop1)$marray MAPop2 <- marray.recapture(chPop2)$marray MAPop3 <- marray.recapture(chPop3)$marray MAPop4 <- marray.recapture(chPop4)$marray # simulate non-breeding data korrf <- c(1,1,1,1,1,1,0.5,0.5,0.5, 2,2,2) # the more north the smaller recovery probability wf <- rbinom(3*4, size=rep(N, 3), prob=m$mean*r*korrf) wirec <- matrix(wf, ncol=3, nrow=4, byrow=TRUE) datax<-list(MArecaptCE=MAPop1, MArecaptEE=MAPop2, MArecaptSE=MAPop3, MArecaptNEE=MAPop4,n=n, J=J, wirec=wirec, npop=npop) bugsData(datax, fileName = "data.txt") minit <- matrix(c(0.61000, 0.16570, 0.22440, 0.02289, 0.30590, 0.67120, 0.01112, 0.64550, 0.34340, 0.04396, 0.08658, 0.86950), ncol=3, nrow=4, byrow=TRUE) minit[,3]<-1-apply(minit[,c(1:2)], 1, sum) minit <- minit[,c(2,3,1)] minit[,3] <- NA init1 <- list(p=runif(npop, 0.3,0.7), m=minit, r=runif(3, 0,0.4)) init2 <- list(p=runif(npop, 0.3,0.7), m=minit, r=runif(3, 0,0.4)) inits<-list(init1, init2) bugsInits(inits, numChains = 2, c("init1.txt", "init2.txt"), digits = 5) mod <- bugs(datax, inits=list(init1, init2), parameters.to.save=c("p", "S", "N", "m", "r") , model.file="model.txt", n.chains=2, n.iter=250000, n.burnin=200000, n.sims = 1000, debug=FALSE, DIC=TRUE, digits=5, codaPkg=FALSE, bugs.directory="c:/Programme/WinBUGS14/", program="WinBUGS", working.directory=getwd()) biaspmean[sim, ] <- mod$summary[1:4,1]-p biaspmedian[sim, ] <- mod$summary[1:4,5]-p biasSmean[sim] <- mod$summary[5,1]-S biasSmedian[sim] <- mod$summary[5,5]-S biasNmean[sim, ] <- mod$summary[6:9,1]-N biasNmedian[sim, ] <- mod$summary[6:9,5]-N biasmmean[sim, ] <- mod$summary[10:21,1]-m$mean biasmmedian[sim, ] <- mod$summary[10:21,5]-m$mean biasrmean[sim, ] <- mod$summary[22:24,1]-r[1:3] biasrmedian[sim, ] <- mod$summary[22:24,5]-r[1:3] # save results bias <- list(biaspmean=biaspmean, biaspmedian=biaspmedian, biasSmean=biasSmean, biasSmedian=biasSmedian, biasNmean=biasNmean, biasNmedian=biasNmedian, biasmmean=biasmmean,biasmmedian=biasmmedian, biasrmean=biasrmean, biasrmedian=biasrmedian) dput(bias, "bias_negcor_rlat.txt") } # close sim #------------------------------------------------------------------------------ # simulate data with positive correlation of recovery probability with latitude of breeding area biasmmean <- matrix(ncol=12, nrow=nsim) biasSmean <- numeric(nsim) biaspmean <- matrix(ncol=4, nrow=nsim) biasrmean <- matrix(ncol=3, nrow=nsim) biasNmean <- matrix(ncol=4, nrow=nsim) biasmmedian <- matrix(ncol=12, nrow=nsim) biasSmedian <- numeric(nsim) biaspmedian <- matrix(ncol=4, nrow=nsim) biasrmedian <- matrix(ncol=3, nrow=nsim) biasNmedian <- matrix(ncol=4, nrow=nsim) for(sim in 1:nsim){ # simulate breeding data pmin1wf <- 1-(1-S)/(1-S*(1-p)) # Wahrscheinlichkeit für mind 1 Wiederfang n <- rbinom(4, N, pmin1wf) chPop1 <- simul.ch(S, p[1], c(n[1], rep(0, J-1)), J, type="S.p") chPop2 <- simul.ch(S, p[2], c(n[2], rep(0, J-1)), J, type="S.p") chPop3 <- simul.ch(S, p[3], c(n[3], rep(0, J-1)), J, type="S.p") chPop4 <- simul.ch(S, p[4], c(n[4], rep(0, J-1)), J, type="S.p") MAPop1 <- marray.recapture(chPop1)$marray MAPop2 <- marray.recapture(chPop2)$marray MAPop3 <- marray.recapture(chPop3)$marray MAPop4 <- marray.recapture(chPop4)$marray # simulate non-breeding data korrf <- c(1,1,1,1,1,1,2,2,2, 0.5,0.5,0.5) # the more north the higher recovery probability wf <- rbinom(3*4, size=rep(N, 3), prob=m$mean*r*korrf) wirec <- matrix(wf, ncol=3, nrow=4, byrow=TRUE) datax<-list(MArecaptCE=MAPop1, MArecaptEE=MAPop2, MArecaptSE=MAPop3, MArecaptNEE=MAPop4,n=n, J=J, wirec=wirec, npop=npop) bugsData(datax, fileName = "data.txt") minit <- matrix(c(0.61000, 0.16570, 0.22440, 0.02289, 0.30590, 0.67120, 0.01112, 0.64550, 0.34340, 0.04396, 0.08658, 0.86950), ncol=3, nrow=4, byrow=TRUE) minit[,3]<-1-apply(minit[,c(1:2)], 1, sum) minit <- minit[,c(2,3,1)] minit[,3] <- NA init1 <- list(p=runif(npop, 0.3,0.7), m=minit, r=runif(3, 0,0.4)) init2 <- list(p=runif(npop, 0.3,0.7), m=minit, r=runif(3, 0,0.4)) inits<-list(init1, init2) bugsInits(inits, numChains = 2, c("init1.txt", "init2.txt"), digits = 5) mod <- bugs(datax, inits=list(init1, init2), parameters.to.save=c("p", "S", "N", "m", "r") , model.file="model.txt", n.chains=2, n.iter=250000, n.burnin=200000, n.sims = 1000, debug=FALSE, DIC=TRUE, digits=5, codaPkg=FALSE, bugs.directory="c:/Programme/WinBUGS14/", program="WinBUGS", working.directory=getwd()) biaspmean[sim, ] <- mod$summary[1:4,1]-p biaspmedian[sim, ] <- mod$summary[1:4,5]-p biasSmean[sim] <- mod$summary[5,1]-S biasSmedian[sim] <- mod$summary[5,5]-S biasNmean[sim, ] <- mod$summary[6:9,1]-N biasNmedian[sim, ] <- mod$summary[6:9,5]-N biasmmean[sim, ] <- mod$summary[10:21,1]-m$mean biasmmedian[sim, ] <- mod$summary[10:21,5]-m$mean biasrmean[sim, ] <- mod$summary[22:24,1]-r[1:3] biasrmedian[sim, ] <- mod$summary[22:24,5]-r[1:3] # save results bias <- list(biaspmean=biaspmean, biaspmedian=biaspmedian, biasSmean=biasSmean, biasSmedian=biasSmedian, biasNmean=biasNmean, biasNmedian=biasNmedian, biasmmean=biasmmean,biasmmedian=biasmmedian, biasrmean=biasrmean, biasrmedian=biasrmedian) dput(bias, "bias_poscor_rlat.txt") } # close sim #------------------------------------------------------------------------------ # simulate data with random between-individual variance in reencounter probability # set the true parameter values N <- rep(50000, 4) p <- c(0.15, 0.2, 0.1, 0.1) S <- 0.5 m <- read.table("mestnoninfprior100920.txt", header=TRUE) mrate <- matrix(m$mean, ncol=3, nrow=4, byrow=TRUE) r <- rep(c(0.001, 0.01, 0.003), 4) # mean of reencounter probability ra <- shapeparameter(r, r/10, r*10)$a rb <- shapeparameter(r, r/10, r*10)$b J <- 8 # n occations x <- seq(0, 1, by=0.01) y <- dbeta(x, ra[1], rb[1]) plot(x,y, type="l") lines(x, dbeta(x, ra[2], rb[2]), col=2) lines(x, dbeta(x, ra[3], rb[3]), col=3) nsim <- 20 biasmmean <- matrix(ncol=12, nrow=nsim) biasSmean <- numeric(nsim) biaspmean <- matrix(ncol=4, nrow=nsim) biasrmean <- matrix(ncol=3, nrow=nsim) biasNmean <- matrix(ncol=4, nrow=nsim) biasmmedian <- matrix(ncol=12, nrow=nsim) biasSmedian <- numeric(nsim) biaspmedian <- matrix(ncol=4, nrow=nsim) biasrmedian <- matrix(ncol=3, nrow=nsim) biasNmedian <- matrix(ncol=4, nrow=nsim) for(sim in 1:nsim){ # simulate breeding data pmin1wf <- 1-(1-S)/(1-S*(1-p)) # Wahrscheinlichkeit für mind 1 Wiederfang n <- rbinom(4, N, pmin1wf) chPop1 <- simul.ch(S, p[1], c(n[1], rep(0, J-1)), J) chPop2 <- simul.ch(S, p[2], c(n[2], rep(0, J-1)), J) chPop3 <- simul.ch(S, p[3], c(n[3], rep(0, J-1)), J) chPop4 <- simul.ch(S, p[4], c(n[4], rep(0, J-1)), J) MAPop1 <- marray.recapture(chPop1)$marray MAPop2 <- marray.recapture(chPop2)$marray MAPop3 <- marray.recapture(chPop3)$marray MAPop4 <- marray.recapture(chPop4)$marray # simulate non-breeding data rindpop1 <- cbind(rbeta(n[1], ra[1], rb[1]), rbeta(n[1], ra[2], rb[2]), rbeta(n[1], ra[3], rb[3])) rindpop2 <- cbind(rbeta(n[2], ra[1], rb[1]), rbeta(n[2], ra[2], rb[2]), rbeta(n[2], ra[3], rb[3])) rindpop3 <- cbind(rbeta(n[3], ra[1], rb[1]), rbeta(n[3], ra[2], rb[2]), rbeta(n[3], ra[3], rb[3])) rindpop4 <- cbind(rbeta(n[4], ra[1], rb[1]), rbeta(n[4], ra[2], rb[2]), rbeta(n[4], ra[3], rb[3])) recpop1 <- c(sum(rbinom(n[1], size=1, prob=m$mean[1]*rindpop1[,1])),sum(rbinom(n[1], size=1, prob=m$mean[2]*rindpop1[,2])), sum(rbinom(n[1], size=1, prob=m$mean[3]*rindpop1[,3]))) recpop2 <- c(sum(rbinom(n[2], size=1, prob=m$mean[4]*rindpop2[,1])),sum(rbinom(n[2], size=1, prob=m$mean[5]*rindpop2[,2])), sum(rbinom(n[2], size=1, prob=m$mean[6]*rindpop2[,3]))) recpop3 <- c(sum(rbinom(n[3], size=1, prob=m$mean[7]*rindpop3[,1])),sum(rbinom(n[3], size=1, prob=m$mean[8]*rindpop3[,2])), sum(rbinom(n[3], size=1, prob=m$mean[9]*rindpop3[,3]))) recpop4 <- c(sum(rbinom(n[4], size=1, prob=m$mean[10]*rindpop4[,1])),sum(rbinom(n[4], size=1, prob=m$mean[11]*rindpop4[,2])), sum(rbinom(n[4], size=1, prob=m$mean[12]*rindpop4[,3]))) wf <- c(recpop1, recpop2, recpop3, recpop4) wirec <- matrix(wf, ncol=3, nrow=4, byrow=TRUE) datax<-list(MArecaptCE=MAPop1, MArecaptEE=MAPop2, MArecaptSE=MAPop3, MArecaptNEE=MAPop4,n=n, J=J, wirec=wirec, npop=npop) bugsData(datax, fileName = "data.txt") minit<-wirec/apply(wirec,1,sum) minit <- matrix(c(0.61000, 0.16570, 0.22440, 0.02289, 0.30590, 0.67120, 0.01112, 0.64550, 0.34340, 0.04396, 0.08658, 0.86950), ncol=3, nrow=4, byrow=TRUE) minit[,3]<-1-apply(minit[,c(1:2)], 1, sum) minit <- minit[,c(2,3,1)] minit[,3] <- NA init1 <- list(p=runif(npop, 0.3,0.7), m=minit, r=runif(3, 0,0.4)) init2 <- list(p=runif(npop, 0.3,0.7), m=minit, r=runif(3, 0,0.4)) inits<-list(init1, init2) bugsInits(inits, numChains = 2, c("init1.txt", "init2.txt"), digits = 5) mod <- bugs(datax, inits=list(init1, init2), parameters.to.save=c("p", "S", "N", "m", "r") , model.file="model.txt", n.chains=2, n.iter=250000, n.burnin=200000, n.sims = 1000, debug=FALSE, DIC=TRUE, digits=5, codaPkg=FALSE, bugs.directory="c:/Programme/WinBUGS14/", program="WinBUGS", working.directory=getwd()) biaspmean[sim, ] <- mod$summary[1:4,1]-p biaspmedian[sim, ] <- mod$summary[1:4,5]-p biasSmean[sim] <- mod$summary[5,1]-S biasSmedian[sim] <- mod$summary[5,5]-S biasNmean[sim, ] <- mod$summary[6:9,1]-N biasNmedian[sim, ] <- mod$summary[6:9,5]-N biasmmean[sim, ] <- mod$summary[10:21,1]-m$mean biasmmedian[sim, ] <- mod$summary[10:21,5]-m$mean biasrmean[sim, ] <- mod$summary[22:24,1]-r[1:3] biasrmedian[sim, ] <- mod$summary[22:24,5]-r[1:3] # save results bias <- list(biaspmean=biaspmean, biaspmedian=biaspmedian, biasSmean=biasSmean, biasSmedian=biasSmedian, biasNmean=biasNmean, biasNmedian=biasNmedian, biasmmean=biasmmean,biasmmedian=biasmmedian, biasrmean=biasrmean, biasrmedian=biasrmedian) dput(bias, "heterogene.r.txt") } # close sim #------------------------------------------------------------------------------ # plot results bias.perfect <- dget("bias_perfectdata.txt") bias.hetp <- dget("bias_p_i.txt") bias.hetpsi <- dget("bias_S_i.txt") bias.poscor <- dget("bias_poscor_rlat.txt") bias.negcor <- dget("bias_negcor_rlat.txt") bias.hetr <- dget("heterogene.r.txt") lower <- function(x) quantile(x, prob=0.025) upper <- function(x) quantile(x, prob=0.975) tsize <- 1.2 # size of title cex.points <- 1.2 cextextin <- 1.2 pdf("fig3.pdf", width=13, height=7) bmp("fig3.bmp", width = 2*480, height = 480) postscript("fig3.eps", width=13, height=7) par(mfrow=c(5, 6), mar=c(0,0,0,0), oma=c(0.2, 7,3,0.2)) # detection probability plot(1:4, seq(-0.1, 0.1, length=4), type="n", xlab="", ylab="", xaxt="n", las=1, cex.lab=1.2, main="", xlim=c(0.5, 4.5)) abline(h=0, lty=3) mtext(expression(p[B]), 2, las=1, line=3.5, cex=1.4) mtext("(1) perfect data", 3, line=0.7, cex=tsize, adj=0) ordind <- c(1,3,4,2) segments(1:4, apply(bias.perfect$biaspmean, 2, lower)[ordind], 1:4, apply(bias.perfect$biaspmean, 2, upper)[ordind], lwd=2) points(1:4, apply(bias.perfect$biaspmean, 2, mean)[ordind], pch=21, bg="white", cex=cex.points) text(1:4, rep(-0.095, 4), c("W ", "E", "S", "NE")[ordind], cex=cextextin) plot(1:4, seq(-0.1, 0.1, length=4), type="n", xlab="", ylab="", xaxt="n", yaxt="n", main="", xlim=c(0.5, 4.5)) abline(h=0, lty=3) mtext("(2) heterogeneous p", 3, line=0.7, cex=tsize, adj=0) segments(1:4, apply(bias.hetp$biaspmean, 2, lower)[ordind], 1:4, apply(bias.hetp$biaspmean, 2, upper)[ordind], lwd=2) points(1:4, apply(bias.hetp$biaspmean, 2, mean)[ordind], pch=21, bg="white", cex=cex.points) text(1:4, rep(-0.095, 4), c("W ", "E", "S", "NE")[ordind], cex=cextextin) plot(1:4, seq(-0.1, 0.1, length=4), type="n", xlab="", ylab="", xaxt="n", yaxt="n", main="", xlim=c(0.5, 4.5)) abline(h=0, lty=3) mtext(expression(paste("(3) heterogeneous ", psi)), 3, line=0.7, cex=tsize, adj=0) segments(1:4, apply(bias.hetpsi$biaspmean, 2, lower)[ordind], 1:4, apply(bias.hetpsi$biaspmean, 2, upper)[ordind], lwd=2) points(1:4, apply(bias.hetpsi$biaspmean, 2, mean)[ordind], pch=21, bg="white", cex=cex.points) text(1:4, rep(-0.095, 4), c("W ", "E", "S", "NE")[ordind], cex=cextextin) plot(1:4, seq(-0.1, 0.1, length=4), type="n", xlab="", ylab="", xaxt="n", yaxt="n", main="", xlim=c(0.5, 4.5)) abline(h=0, lty=3) mtext("(4) pos cor r-latitude", 3, line=0.7, cex=tsize, adj=0) segments(1:4, apply(bias.poscor$biaspmean, 2, lower)[ordind], 1:4, apply(bias.poscor$biaspmean, 2, upper)[ordind], lwd=2) text(1:4, rep(-0.095, 4), c("W ", "E", "S", "NE")[ordind], cex=cextextin) points(1:4, apply(bias.poscor$biaspmean, 2, mean)[ordind], pch=21, bg="white", cex=cex.points) plot(1:4, seq(-0.1, 0.1, length=4), type="n", xlab="", ylab="", xaxt="n", yaxt="n", main="", xlim=c(0.5, 4.5)) abline(h=0, lty=3) mtext("(5) neg cor r-latitude", 3, line=0.7, cex=tsize, adj=0) segments(1:4, apply(bias.negcor$biaspmean, 2, lower)[ordind], 1:4, apply(bias.negcor$biaspmean, 2, upper)[ordind], lwd=2) points(1:4, apply(bias.negcor$biaspmean, 2, mean)[ordind], pch=21, bg="white", cex=cex.points) text(1:4, rep(-0.095, 4), c("W ", "E", "S", "NE")[ordind], cex=cextextin) plot(1:4, seq(-0.1, 0.1, length=4), type="n", xlab="", ylab="", xaxt="n", yaxt="n", main="", xlim=c(0.5, 4.5)) abline(h=0, lty=3) mtext("(6) heterogeneous r", 3, line=0.7, cex=tsize, adj=0) segments(1:4, apply(bias.hetr$biaspmean, 2, lower)[ordind], 1:4, apply(bias.hetr$biaspmean, 2, upper)[ordind], lwd=2) points(1:4, apply(bias.hetr$biaspmean, 2, mean)[ordind], pch=21, bg="white", cex=cex.points) text(1:4, rep(-0.095, 4), c("W ", "E", "S", "NE")[ordind], cex=cextextin) # return probabilities plot(1:4, seq(-0.2, 0.2, length=4), type="n", xlab="", ylab="", xaxt="n", las=1, cex.lab=1.2, main="", xlim=c(0.5, 4.5)) abline(h=0, lty=3) mtext(expression(psi), 2, las=1, line=3.5, cex=1.4) segments(2.5, lower(bias.perfect$biasSmean), 2.5, upper(bias.perfect$biasSmean), lwd=2) points(2.5, median(bias.perfect$biasSmean), pch=21, bg="white", cex=cex.points) plot(1:4, seq(-0.2, 0.2, length=4), type="n", xlab="", ylab="", xaxt="n", yaxt="n", main="", xlim=c(0.5, 4.5)) abline(h=0, lty=3) segments(2.5, lower(bias.hetp$biasSmean), 2.5, upper(bias.hetp$biasSmean), lwd=2) points(2.5, median(bias.hetp$biasSmean), pch=21, bg="white", cex=cex.points) plot(1:4, seq(-0.2, 0.2, length=4), type="n", xlab="", ylab="", xaxt="n", yaxt="n", main="", xlim=c(0.5, 4.5)) abline(h=0, lty=3) segments(2.5, lower(bias.hetpsi$biasSmean), 2.5, upper(bias.hetpsi$biasSmean), lwd=2) points(2.5, median(bias.hetpsi$biasSmean), pch=21, bg="white", cex=cex.points) plot(1:4, seq(-0.2, 0.2, length=4), type="n", xlab="", ylab="", xaxt="n", yaxt="n", main="", xlim=c(0.5, 4.5)) abline(h=0, lty=3) segments(2.5, lower(bias.poscor$biasSmean), 2.5, upper(bias.poscor$biasSmean), lwd=2) points(2.5, median(bias.poscor$biasSmean), pch=21, bg="white", cex=cex.points) plot(1:4, seq(-0.2, 0.2, length=4), type="n", xlab="", ylab="", xaxt="n", yaxt="n", main="", xlim=c(0.5, 4.5)) abline(h=0, lty=3) segments(2.5, lower(bias.negcor$biasSmean), 2.5, upper(bias.negcor$biasSmean), lwd=2) points(2.5, median(bias.negcor$biasSmean), pch=21, bg="white", cex=cex.points) plot(1:4, seq(-0.2, 0.2, length=4), type="n", xlab="", ylab="", xaxt="n", yaxt="n", main="", xlim=c(0.5, 4.5)) abline(h=0, lty=3) segments(2.5, lower(bias.hetr$biasSmean), 2.5, upper(bias.hetr$biasSmean), lwd=2) points(2.5, median(bias.hetr$biasSmean), pch=21, bg="white", cex=cex.points) # migration rates ordind <- c(3, 1, 2, 9, 7,8,12, 10, 11, 6, 4, 5) segheight <- 0.4 textheight <- 0.46 plot(1:12, seq(-0.5, 0.5, length=12), type="n", xlab="", ylab="", xaxt="n", las=1, cex.lab=1.2, main="", xlim=c(0.5, 12.5)) abline(h=0, lty=3) mtext(expression(m[BD]), 2, las=1, line=3.5, cex=1.4) segments(1:12, apply(bias.perfect$biasmmean, 2, lower)[ordind], 1:12, apply(bias.perfect$biasmmean, 2, upper)[ordind], lwd=2, col=c(1,1,1, grey(0.5), grey(0.5), 1, 1, 1, 1, grey(0.5), grey(0.5), 1)[ordind]) points(1:12, apply(bias.perfect$biasmmean, 2, mean)[ordind], pch=21, bg="white", cex=cex.points) segments(c(1,4,7,10), segheight, c(1,4,7,10)+2, segheight) text(c(2, 5, 8, 11), textheight, c("W ", "E", "S", "NE"), cex=cextextin) #text(1:12, rep(-0.35, 12), c("W C-EU\nto CENTRAL", "W C-EU\nto EAST", "W C-EU\nto WEST", "E-EU\nto CENTRAL", "E-EU\nto EAST", "E-EU\nto WEST", # "S C-EU\nto CENTRAL", "S C-EU\nto EAST", "S C-EU\nto WEST", "NE C-EU\nto CENTRAL", "NE C-EU\nto EAST", "NE C-EU\nto WEST"), cex=0.8) text(1:12, -textheight, rep(c("w", "c", "e"), 4), cex=cextextin) plot(1:12, seq(-0.5, 0.5, length=12), type="n", xlab="", ylab="", xaxt="n",yaxt="n", main="", xlim=c(0.5, 12.5)) abline(h=0, lty=3) segments(1:12, apply(bias.hetp$biasmmean, 2, lower)[ordind], 1:12, apply(bias.hetp$biasmmean, 2, upper)[ordind], lwd=2, col=c(1,1,1, grey(0.5), grey(0.5), 1, 1, 1, 1, grey(0.5), grey(0.5), 1)[ordind]) points(1:12, apply(bias.hetp$biasmmean, 2, mean)[ordind], pch=21, bg="white", cex=cex.points) segments(c(1,4,7,10), segheight, c(1,4,7,10)+2, segheight) text(c(2, 5, 8, 11), textheight, c("W ", "E", "S", "NE"), cex=cextextin) text(1:12, -textheight, rep(c("w", "c", "e"), 4), cex=cextextin) plot(1:12, seq(-0.5, 0.5, length=12), type="n", xlab="", ylab="", xaxt="n",yaxt="n", main="", xlim=c(0.5, 12.5)) abline(h=0, lty=3) segments(1:12, apply(bias.hetpsi$biasmmean, 2, lower)[ordind], 1:12, apply(bias.hetpsi$biasmmean, 2, upper)[ordind], lwd=2, col=c(1,1,1, grey(0.5), grey(0.5), 1, 1, 1, 1, grey(0.5), grey(0.5), 1)[ordind]) points(1:12, apply(bias.hetpsi$biasmmean, 2, mean)[ordind], pch=21, bg="white", cex=cex.points) segments(c(1,4,7,10), segheight, c(1,4,7,10)+2, segheight) text(c(2, 5, 8, 11), textheight, c("W ", "E", "S", "NE"), cex=cextextin) text(1:12, -textheight, rep(c("w", "c", "e"), 4), cex=cextextin) plot(1:12, seq(-0.5, 0.5, length=12), type="n", xlab="", ylab="", xaxt="n",yaxt="n", main="", xlim=c(0.5, 12.5)) abline(h=0, lty=3) segments(1:12, apply(bias.poscor$biasmmean, 2, lower)[ordind], 1:12, apply(bias.poscor$biasmmean, 2, upper)[ordind], lwd=2, col=c(1,1,1, grey(0.5), grey(0.5), 1, 1, 1, 1, grey(0.5), grey(0.5), 1)[ordind]) points(1:12, apply(bias.poscor$biasmmean, 2, mean)[ordind], pch=21, bg="white", cex=cex.points) segments(c(1,4,7,10), segheight, c(1,4,7,10)+2, segheight) text(c(2, 5, 8, 11), textheight, c("W ", "E", "S", "NE"), cex=cextextin) text(1:12, -textheight, rep(c("w", "c", "e"), 4), cex=cextextin) plot(1:12, seq(-0.5, 0.5, length=12), type="n", xlab="", ylab="", xaxt="n",yaxt="n", main="", xlim=c(0.5, 12.5)) abline(h=0, lty=3) segments(1:12, apply(bias.negcor$biasmmean, 2, lower)[ordind], 1:12, apply(bias.negcor$biasmmean, 2, upper)[ordind], lwd=2, col=c(1,1,1, grey(0.5), grey(0.5), 1, 1, 1, 1, grey(0.5), grey(0.5), 1)[ordind]) points(1:12, apply(bias.negcor$biasmmean, 2, mean)[ordind], pch=21, bg="white", cex=cex.points) segments(c(1,4,7,10), segheight+c(0,0,-0.3, -0.3), c(1,4,7,10)+2, segheight+c(0,0,-0.3, -0.3)) text(c(2, 5, 8, 11), textheight+c(0,0,-0.3, -0.3), c("W ", "E", "S", "NE"), cex=cextextin) text(1:12, -textheight+c(rep(0, 8), 0.3, 0, 0, 0.16), rep(c("w", "c", "e"), 4), cex=cextextin) plot(1:12, seq(-0.5, 0.5, length=12), type="n", xlab="", ylab="", xaxt="n",yaxt="n", main="", xlim=c(0.5, 12.5)) abline(h=0, lty=3) segments(1:12, apply(bias.hetr$biasmmean, 2, lower)[ordind], 1:12, apply(bias.hetr$biasmmean, 2, upper)[ordind], lwd=2, col=c(1,1,1, grey(0.5), grey(0.5), 1, 1, 1, 1, grey(0.5), grey(0.5), 1)[ordind]) points(1:12, apply(bias.hetr$biasmmean, 2, mean)[ordind], pch=21, bg="white", cex=cex.points) segments(c(1,4,7,10), segheight, c(1,4,7,10)+2, segheight) text(c(2, 5, 8, 11), textheight, c("W ", "E", "S", "NE"), cex=cextextin) text(1:12, -textheight, rep(c("w", "c", "e"), 4), cex=cextextin) # recovery probability ordind <- c(3,1,2) plot(1:3, seq(-0.1, 0.1, length=3), type="n", xlab="", ylab="", xaxt="n", las=1, cex.lab=1.2, main="", xlim=c(0.5, 3.5)) abline(h=0, lty=3) mtext(expression(r[D]), 2, las=1, line=3.5, cex=1.4) segments(1:3, apply(bias.perfect$biasrmean, 2, lower)[ordind], 1:3, apply(bias.perfect$biasrmean, 2, upper)[ordind], lwd=2) points(1:3, apply(bias.perfect$biasrmean, 2, mean)[ordind], pch=21, bg="white", cex=cex.points) text(1:3, rep(-0.095, 3), c("c", "e", "w")[ordind], cex=cextextin) plot(1:3, seq(-0.1, 0.1, length=3), type="n", xlab="", ylab="", xaxt="n",yaxt="n", main="", xlim=c(0.5, 3.5)) abline(h=0, lty=3) segments(1:3, apply(bias.hetp$biasrmean, 2, lower)[ordind], 1:3, apply(bias.hetp$biasrmean, 2, upper)[ordind], lwd=2) points(1:3, apply(bias.hetp$biasrmean, 2, mean)[ordind], pch=21, bg="white", cex=cex.points) text(1:3, rep(-0.095, 3), c("c", "e", "w")[ordind], cex=cextextin) plot(1:3, seq(-0.1, 0.1, length=3), type="n", xlab="", ylab="", xaxt="n",yaxt="n", main="", xlim=c(0.5, 3.5)) abline(h=0, lty=3) segments(1:3, apply(bias.hetpsi$biasrmean, 2, lower)[ordind], 1:3, apply(bias.hetpsi$biasrmean, 2, upper)[ordind], lwd=2) points(1:3, apply(bias.hetpsi$biasrmean, 2, mean)[ordind], pch=21, bg="white", cex=cex.points) text(1:3, rep(-0.095, 3), c("c", "e", "w")[ordind], cex=cextextin) plot(1:3, seq(-0.1, 0.1, length=3), type="n", xlab="", ylab="", xaxt="n",yaxt="n", main="", xlim=c(0.5, 3.5)) abline(h=0, lty=3) segments(1:3, apply(bias.poscor$biasrmean, 2, lower)[ordind], 1:3, apply(bias.poscor$biasrmean, 2, upper)[ordind], lwd=2) points(1:3, apply(bias.poscor$biasrmean, 2, mean)[ordind], pch=21, bg="white", cex=cex.points) text(1:3, rep(-0.095, 3), c("c", "e", "w")[ordind], cex=cextextin) plot(1:3, seq(-0.1, 0.1, length=3), type="n", xlab="", ylab="", xaxt="n",yaxt="n", main="", xlim=c(0.5, 3.5)) abline(h=0, lty=3) segments(1:3, apply(bias.negcor$biasrmean, 2, lower)[ordind], 1:3, apply(bias.negcor$biasrmean, 2, upper)[ordind], lwd=2) points(1:3, apply(bias.negcor$biasrmean, 2, mean)[ordind], pch=21, bg="white", cex=cex.points) text(1:3, rep(-0.095, 3), c("c", "e", "w")[ordind], cex=cextextin) plot(1:3, seq(-0.1, 0.1, length=3), type="n", xlab="", ylab="", xaxt="n",yaxt="n", main="", xlim=c(0.5, 3.5)) abline(h=0, lty=3) segments(1:3, apply(bias.hetr$biasrmean, 2, lower)[ordind], 1:3, apply(bias.hetr$biasrmean, 2, upper)[ordind], lwd=2) points(1:3, apply(bias.hetr$biasrmean, 2, mean)[ordind], pch=21, bg="white", cex=cex.points) text(1:3, rep(-0.095, 3), c("c", "e", "w")[ordind], cex=cextextin) dev.off() #------------------------------------------------------------------------------- # number of ringed birds plot(1:4, seq(-0.5, 0.5, length=4), type="n", xlab="", ylab="", xaxt="n", las=1, cex.lab=1.2, main="", xlim=c(0.5, 4.5)) abline(h=0, lty=3) mtext(expression(N[B]), 2, las=1, line=3.5, cex=1.4) segments(1:4, apply(bias.perfect$biasNmean, 2, lower)/50000, 1:4, apply(bias.perfect$biasNmean, 2, upper)/50000, lwd=2) points(1:4, apply(bias.perfect$biasNmean, 2, mean)/50000, pch=21, bg="white") text(1:4, rep(-0.49, 4), c("W C-EU", "E-EU", "S C-EU", "NE C-EU")) plot(1:4, seq(-0.5, 0.5, length=4), type="n", xlab="", ylab="", xaxt="n", yaxt="n", main="", xlim=c(0.5, 4.5)) abline(h=0, lty=3) segments(1:4, apply(bias.hetp$biasNmean, 2, lower)/50000, 1:4, apply(bias.hetp$biasNmean, 2, upper)/50000, lwd=2) points(1:4, apply(bias.hetp$biasNmean, 2, mean)/50000, pch=21, bg="white") text(1:4, rep(-0.49, 4), c("W C-EU", "E-EU", "S C-EU", "NE C-EU")) plot(1:4, seq(-0.5, 0.5, length=4), type="n", xlab="", ylab="", xaxt="n", yaxt="n", main="", xlim=c(0.5, 4.5)) abline(h=0, lty=3) segments(1:4, apply(bias.hetpsi$biasNmean, 2, lower)/50000, 1:4, apply(bias.hetpsi$biasNmean, 2, upper)/50000, lwd=2) points(1:4, apply(bias.hetpsi$biasNmean, 2, mean)/50000, pch=21, bg="white") text(1:4, rep(-0.49, 4), c("W C-EU", "E-EU", "S C-EU", "NE C-EU")) plot(1:4, seq(-0.5, 0.5, length=4), type="n", xlab="", ylab="", xaxt="n", yaxt="n", main="", xlim=c(0.5, 4.5)) abline(h=0, lty=3) segments(1:4, apply(bias.poscor$biasNmean, 2, lower)/50000, 1:4, apply(bias.poscor$biasNmean, 2, upper)/50000, lwd=2) points(1:4, apply(bias.poscor$biasNmean, 2, mean)/50000, pch=21, bg="white") text(1:4, rep(-0.49, 4), c("W C-EU", "E-EU", "S C-EU", "NE C-EU")) plot(1:4, seq(-0.5, 0.5, length=4), type="n", xlab="", ylab="", xaxt="n", yaxt="n", main="", xlim=c(0.5, 4.5)) abline(h=0, lty=3) segments(1:4, apply(bias.negcor$biasNmean, 2, lower)/50000, 1:4, apply(bias.negcor$biasNmean, 2, upper)/50000, lwd=2) points(1:4, apply(bias.negcor$biasNmean, 2, mean)/50000, pch=21, bg="white") text(1:4, rep(-0.49, 4), c("W C-EU", "E-EU", "S C-EU", "NE C-EU")) plot(1:4, seq(-0.5, 0.5, length=4), type="n", xlab="", ylab="", xaxt="n", yaxt="n", main="", xlim=c(0.5, 4.5)) abline(h=0, lty=3) segments(1:4, apply(bias.hetr$biasNmean, 2, lower)/50000, 1:4, apply(bias.hetr$biasNmean, 2, upper)/50000, lwd=2) points(1:4, apply(bias.hetr$biasNmean, 2, mean)/50000, pch=21, bg="white") text(1:4, rep(-0.49, 4), c("W C-EU", "E-EU", "S C-EU", "NE C-EU"))