# R-Code to supplementary material 3 # Korner-Nievergelt, F., Liechti, F., Thorup, K. (2013): # A bird distribution model for ring recovery data: Where do the European Robins go? # submitted to Methods in Ecology and Evolution #------------------------------------------------------------------------------- # simulation study to assess minimal sample size needed # note: here we use two different models: V1 and V2 # V1 allows that all sets of birds (ij) released at the same region in the same month behave independently with respect to migration # V2 constrains the birds released during the three winter months and during June and July respectively to behave similarily # in the paper, we only present model V2 # the model does not run in WinBugs! #------------------------------------------------------------------------------- # note: please report bugs and comments to fraenzi.korner@vogelwarte.ch #------------------------------------------------------------------------------- setwd("....") bugsworkingdir <- getwd() # load libraries library(arm) library(R2jags) #------------------- ---------------------------------------------------------- #------------------------------------------------------------------------------- # Model definition sink(file.path(bugsworkingdir, "distributionmodelV2.txt")) cat(" model{ ### Likelihood for(i in 1:(npop*nrelocc)){ recmatlong[i,1:33]~ dmulti(probmatlong[i,1:33], nreleased[i]) } # survival matrix A <- pow(s,11)*(1-s)/(1-pow(s,12)) for(j in 1:nrelocc){ # diagonal survivalmatcomp[j, j] <- A } # upper diagonal for(j in 1:(nrelocc-1)){ for(lc in (j+1):nrelocc){ survivalmatcomp[j, lc] <- pow(s,(lc-j-1))*(1-s) + pow(s,(lc-j))*A } } # lower diagonal for(j in 2:nrelocc){ for(lc in 1:(j-1)){ survivalmatcomp[j, lc] <- pow(s,(11-j+lc))*(1-s) + pow(s,(12-j+lc))*A } } # merge recovery seasons (winter and summer) for(j in 1:nrelocc){ survivalmat[j,1] <- survivalmatcomp[j,1] + survivalmatcomp[j,2] + survivalmatcomp[j,12] for(l in 2:4){ survivalmat[j,l] <- survivalmatcomp[j,l+1] } survivalmat[j,5] <- sum(survivalmatcomp[j,6:8]) for(l in 6:8){ survivalmat[j,l] <- survivalmatcomp[j,l+3] } } # constraint for distributions for(i in 1:npop){ for(j in 1:nrelocc){ # proportion of birds in A m0[i,j,1,1] <- 0 # winter m0[i,j,1,2] <- 0 # march m0[i,j,1,3] <- m0a[i,relseason[j],1,3] # april m0[i,j,1,4] <- m0a[i,relseason[j],1,4] # mai m0[i,j,1,5] <- m0a[i,relseason[j],1,5] # juni-aug m0[i,j,1,6] <- 0 # sept m0[i,j,1,7] <- 0 # oct m0[i,j,1,8] <- 0 # nov for(l in 1:nrecocc){ m0[i,j,2,l] <- m0a[i,relseason[j],2,l] # proportion of birds in B } for(k in 3:ndest){ # in C and D m0[i,j,k,1] <- m0a[i,relseason[j],k,1] # winter m0[i,j,k,2] <- m0a[i,relseason[j],k,2] # march m0[i,j,k,3] <- m0a[i,relseason[j],k,3] # april m0[i,j,k,4] <- m0a[i,relseason[j],k,4] # mai m0[i,j,k,5] <- 0 # juni-aug m0[i,j,k,6] <- m0a[i,relseason[j],k,6] # sept m0[i,j,k,7] <- m0a[i,relseason[j],k,7] # oct m0[i,j,k,8] <- m0a[i,relseason[j],k,8] # nov } #k } #j } #i # proportions are equal for those birds ringed during the winter and summer months for(i in 1:npop){ for(j in 1:8){ # proportion of birds in A m0a[i,j,1,1] <- 0 # winter m0a[i,j,1,2] <- 0 # march m0a[i,j,1,3] ~ dunif(0,1) # april m0a[i,j,1,4] ~ dunif(0,1) # mai m0a[i,j,1,5] ~ dunif(0,1) # juni-aug m0a[i,j,1,6] <- 0 # sept m0a[i,j,1,7] <- 0 # oct m0a[i,j,1,8] <- 0 # nov for(l in 1:nrecocc){ m0a[i,j,2,l] ~ dunif(0,1) # proportion of birds in B } for(k in 3:ndest){ # in C and D m0a[i,j,k,1] ~ dunif(0,1) # winter m0a[i,j,k,2] ~ dunif(0,1) # march m0a[i,j,k,3] ~ dunif(0,1) # april m0a[i,j,k,4] ~ dunif(0,1) # mai m0a[i,j,k,5] <- 0 # juni-aug m0a[i,j,k,6] ~ dunif(0,1) # sept m0a[i,j,k,7] ~ dunif(0,1) # oct m0a[i,j,k,8] ~ dunif(0,1) # nov } #k } #j } #i # scale to sum of m over destination areas = 1 for(i in 1:npop){ for(j in 1:nrelocc){ for(l in 1:nrecocc){ summ[i,j,l] <- sum(m0[i,j,1:ndest,l]) for(k in 1:ndest){ m[i,j,k,l] <- m0[i,j,k,l]/summ[i,j,l] } } } } # fill up probability matrix for(i in 1:npop){ for(j in 1:nrelocc){ for(k in 1:ndest){ for(l in 1:nrecocc){ probmat[i,j,k,l] <- survivalmat[j,l] * r[k] * m[i,j,k,l] } } } } # rearrange probabilities into long format for(i in 1:npop){ for(j in 1:nrelocc){ for(k in 1:ndest){ probmatlong[(i-1)*nrelocc+j, ((k-1)*nrecocc+1):(k*nrecocc)] <- probmat[i,j,k,1:nrecocc] } probmatlong[(i-1)*nrelocc+j, nrecocc*ndest+1] <- 1-sum(probmatlong[(i-1)*nrelocc+j, 1:(nrecocc*ndest)]) } } ### priors s ~ dunif(0,1) for(k in 1:ndest){ r[k] ~ dunif(0,1) } } ",fill=TRUE) sink() sink(file.path(bugsworkingdir, "distributionmodelV1.txt")) cat(" model{ ### Likelihood for(i in 1:(npop*nrelocc)){ recmatlong[i,1:33]~ dmulti(probmatlong[i,1:33], nreleased[i]) } # survival matrix A <- pow(s,11)*(1-s)/(1-pow(s,12)) for(j in 1:nrelocc){ # diagonal survivalmatcomp[j, j] <- A } # upper diagonal for(j in 1:(nrelocc-1)){ for(lc in (j+1):nrelocc){ survivalmatcomp[j, lc] <- pow(s,(lc-j-1))*(1-s) + pow(s,(lc-j))*A } } # lower diagonal for(j in 2:nrelocc){ for(lc in 1:(j-1)){ survivalmatcomp[j, lc] <- pow(s,(11-j+lc))*(1-s) + pow(s,(12-j+lc))*A } } # merge recovery seasons (winter and summer) for(j in 1:nrelocc){ survivalmat[j,1] <- survivalmatcomp[j,1] + survivalmatcomp[j,2] + survivalmatcomp[j,12] for(l in 2:4){ survivalmat[j,l] <- survivalmatcomp[j,l+1] } survivalmat[j,5] <- sum(survivalmatcomp[j,6:8]) for(l in 6:8){ survivalmat[j,l] <- survivalmatcomp[j,l+3] } } # constraint for distributions for(i in 1:npop){ for(j in 1:nrelocc){ # proportion of birds in A m0[i,j,1,1] <- 0 # winter m0[i,j,1,2] <- 0 # march m0[i,j,1,3] ~ dunif(0,1) # april m0[i,j,1,4] ~ dunif(0,1) # mai m0[i,j,1,5] ~ dunif(0,1) # juni-aug m0[i,j,1,6] <- 0 # sept m0[i,j,1,7] <- 0 # oct m0[i,j,1,8] <- 0 # nov for(l in 1:nrecocc){ m0[i,j,2,l] ~ dunif(0,1) # proportion of birds in B } for(k in 3:ndest){ # in C and D m0[i,j,k,1] ~ dunif(0,1) # winter m0[i,j,k,2] ~ dunif(0,1) # march m0[i,j,k,3] ~ dunif(0,1) # april m0[i,j,k,4] ~ dunif(0,1) # mai m0[i,j,k,5] <- 0 # juni-aug m0[i,j,k,6] ~ dunif(0,1) # sept m0[i,j,k,7] ~ dunif(0,1) # oct m0[i,j,k,8] ~ dunif(0,1) # nov } } } # scale to sum of m over destination areas = 1 for(i in 1:npop){ for(j in 1:nrelocc){ for(l in 1:nrecocc){ summ[i,j,l] <- sum(m0[i,j,1:ndest,l]) for(k in 1:ndest){ m[i,j,k,l] <- m0[i,j,k,l]/summ[i,j,l] } } } } # fill up probability matrix for(i in 1:npop){ for(j in 1:nrelocc){ for(k in 1:ndest){ for(l in 1:nrecocc){ probmat[i,j,k,l] <- survivalmat[j,l] * r[k] * m[i,j,k,l] } } } } # rearrange probabilities into long format for(i in 1:npop){ for(j in 1:nrelocc){ for(k in 1:ndest){ probmatlong[(i-1)*nrelocc+j, ((k-1)*nrecocc+1):(k*nrecocc)] <- probmat[i,j,k,1:nrecocc] } probmatlong[(i-1)*nrelocc+j, nrecocc*ndest+1] <- 1-sum(probmatlong[(i-1)*nrelocc+j, 1:(nrecocc*ndest)]) } } ### priors s ~ dunif(0,1) for(k in 1:ndest){ r[k] ~ dunif(0,1) } } ",fill=TRUE) sink() #---------------------------------------------------------------------------------- # Define parameters to be monitored parameters <- c("m","s","r") initfunV2 <- function(){ ncells <- npop*8*ndest*nrecocc m0a <- array(runif(ncells, 0,1), dim=c(npop, 8, ndest, nrecocc)) for(i in 1:npop){ for(j in 1:8){ # proportion of birds in A m0a[i,j,1,1] <- NA # winter m0a[i,j,1,2] <- NA # march m0a[i,j,1,6] <- NA # sept m0a[i,j,1,7] <- NA # oct m0a[i,j,1,8] <- NA # nov for(k in 3:ndest){ # in C and D m0a[i,j,k,5] <- NA # juni-aug } } } list(r=runif(ndest, 0.01, 0.1), s=runif(1, 0,1), m0a=m0a) } initfunV1 <- function(){ ncells <- npop*nrelocc*ndest*nrecocc m0 <- array(runif(ncells, 0,1), dim=c(npop, nrelocc, ndest, nrecocc)) for(i in 1:npop){ for(j in 1:nrelocc){ # proportion of birds in A m0[i,j,1,1] <- NA # winter m0[i,j,1,2] <- NA # march m0[i,j,1,6] <- NA # sept m0[i,j,1,7] <- NA # oct m0[i,j,1,8] <- NA # nov for(k in 3:ndest){ # in C and D m0[i,j,k,5] <- NA # juni-aug } } } list(r=runif(ndest, 0.01, 0.1), s=runif(1, 0,1), m0=m0) } #-------------------------------------------------------------------------------- # MCMC settings niter <- 20000 nthin <- 2 nburn <- 5000 nchains <- 2 # set fixed true model parameters s and r s <- 0.93 # survival matrix survivalmatcomp <- matrix(ncol=12, nrow=12) survivalmat <- matrix(ncol=8, nrow=12) A <- s^11*(1-s)/(1-s^12) for(j in 1:nrelocc){ # diagonal survivalmatcomp[j, j] <- A } # upper diagonal for(j in 1:(nrelocc-1)){ for(lc in (j+1):nrelocc){ survivalmatcomp[j, lc] <- s^(lc-j-1)*(1-s) + s^(lc-j)*A } } # lower diagonal for(j in 2:nrelocc){ for(lc in 1:(j-1)){ survivalmatcomp[j, lc] <- s^(11-j+lc)*(1-s) + s^(12-j+lc)*A } } # merge recovery seasons (winter and summer) for(j in 1:nrelocc){ survivalmat[j,1] <- survivalmatcomp[j,1] + survivalmatcomp[j,2] + survivalmatcomp[j,12] for(l in 2:4){ survivalmat[j,l] <- survivalmatcomp[j,l+1] } survivalmat[j,5] <- sum(survivalmatcomp[j,6:8]) for(l in 6:8){ survivalmat[j,l] <- survivalmatcomp[j,l+3] } } r <- c(0.0006, 0.001, 0.0015, 0.002) npop <- 2 # number of release regions nrelocc <- 12 # number of release occasions ndest <- 4 # number of destination areas (recovery regions) nrecocc <- 8 # number of recovery occasions #------------------------------------------------ # start simulations nsim <- 20 result <- expand.grid(n=c(50000, 10000, 1000), n0=c(0,2,5)) # prepare vectors for the estimates estsV1 <- array(dim=c(nsim, nrow(result))) # dim1=number of simulations, dim2 = number of settings estsV2 <- array(dim=c(nsim, nrow(result))) lowerV1.s <- array(dim=c(nsim, nrow(result))) upperV1.s <- array(dim=c(nsim, nrow(result))) lowerV2.s <- array(dim=c(nsim, nrow(result))) upperV2.s <- array(dim=c(nsim, nrow(result))) estrV1 <- array(dim=c(nsim, nrow(result), 4)) # dim1=number of simulations, dim2 = number of settings, dim3= number of parameters estrV2 <- array(dim=c(nsim, nrow(result), 4)) lower.rV1 <- array(dim=c(nsim, nrow(result), 4)) lower.rV2 <- array(dim=c(nsim, nrow(result), 4)) upper.rV1 <- array(dim=c(nsim, nrow(result), 4)) upper.rV2 <- array(dim=c(nsim, nrow(result), 4)) estmV1 <- array(dim=c(nsim, nrow(result), npop*nrelocc*ndest*nrecocc)) # dim1=number of simulations, dim2 = number of settings, dim3= number of parameters estmV2 <- array(dim=c(nsim, nrow(result), npop*nrelocc*ndest*nrecocc)) lower.mV1 <- array(dim=c(nsim, nrow(result), npop*nrelocc*ndest*nrecocc)) # dim1=number of simulations, dim2 = number of settings, dim3= number of parameters lower.mV2 <- array(dim=c(nsim, nrow(result), npop*nrelocc*ndest*nrecocc)) upper.mV1 <- array(dim=c(nsim, nrow(result), npop*nrelocc*ndest*nrecocc)) # dim1=number of simulations, dim2 = number of settings, dim3= number of parameters upper.mV2 <- array(dim=c(nsim, nrow(result), npop*nrelocc*ndest*nrecocc)) truemV1 <- array(dim=c(nsim, nrow(result), npop*nrelocc*ndest*nrecocc)) truemV2 <- array(dim=c(nsim, nrow(result), npop*nrelocc*ndest*nrecocc)) nrecoveriesV1 <- array(0, dim=c(nsim, nrow(result))) # prepare vector with the number of recoveries nrecoveriesV2 <- array(0, dim=c(nsim, nrow(result))) # prepare vector with the number of recoveries set.seed(45165) for(rowres in 4:nrow(result)){ nreleased <- rep(result$n[rowres], npop*nrelocc) if(result$n0[rowres]>0) nreleased[sample(1:length(nreleased),result$n0[rowres], replace=FALSE)] <- 3 nrsims <- 1 for(simulnr in 1:nsim){ if(nrsims>50) next nrsims <- 1 while(nrecoveriesV2[simulnr, rowres]<5){ m0 <- array(dim=c(npop, nrelocc, ndest, nrecocc)) m <- array(dim=c(npop, nrelocc, ndest, nrecocc)) mV2 <- array(dim=c(npop, nrelocc, ndest, nrecocc)) summV1 <- array(dim=c(npop, nrelocc, nrecocc)) summV2 <- array(dim=c(npop, nrelocc, nrecocc)) probmat <- array(dim=c(npop, nrelocc, ndest, nrecocc)) probmatlong <- matrix(nrow=npop*nrelocc, ncol=ndest*nrecocc+1) recmatlong <- matrix(nrow=npop*nrelocc, ncol=ndest*nrecocc+1) probmatV2 <- array(dim=c(npop, nrelocc, ndest, nrecocc)) probmatlongV2 <- matrix(nrow=npop*nrelocc, ncol=ndest*nrecocc+1) recmatlongV2 <- matrix(nrow=npop*nrelocc, ncol=ndest*nrecocc+1) # simulate distribution for(i in 1:npop){ for(j in 1:nrelocc){ # proportion of birds in A m0[i,j,1,1] <- 0 # winter m0[i,j,1,2] <- 0 # march m0[i,j,1,3] <-runif(1, 0,1) # april m0[i,j,1,4] <-runif(1, 0,1) # mai m0[i,j,1,5] <-runif(1, 0,1) # juni-aug m0[i,j,1,6] <- 0 # sept m0[i,j,1,7] <- 0 # oct m0[i,j,1,8] <- 0 # nov for(l in 1:nrecocc){ m0[i,j,2,l] <-runif(1, 0,1) # proportion of birds in B } for(k in 3:ndest){ # in C and D m0[i,j,k,1] <-runif(1, 0,1) # winter m0[i,j,k,2] <-runif(1, 0,1) # march m0[i,j,k,3] <-runif(1, 0,1) # april m0[i,j,k,4] <- 0 # mai m0[i,j,k,5] <- 0 # juni-aug m0[i,j,k,6] <-runif(1, 0,1) # sept m0[i,j,k,7] <-runif(1, 0,1) # oct m0[i,j,k,8] <-runif(1, 0,1) # nov } } } m0V1 <- m0 m0V2 <- m0 for(i in 1:npop){ for(l in 1:nrecocc){ for(k in 1:ndest){ m0V2[i,1,k,l] <- m0V2[i,2,k,l] m0V2[i,12,k,l] <- m0V2[i,2,k,l] m0V2[i,6,k,l] <- m0V2[i,7,k,l] m0V2[i,8,k,l] <- m0V2[i,7,k,l] } } } # scale to sum of m over destination areas = 1 for(i in 1:npop){ for(j in 1:nrelocc){ for(l in 1:nrecocc){ summV1[i,j,l] <- sum(m0[i,j,1:ndest,l]) summV2[i,j,l] <- sum(m0V2[i,j,1:ndest,l]) for(k in 1:ndest){ m[i,j,k,l] <- m0[i,j,k,l]/summV1[i,j,l] mV2[i,j,k,l] <- m0V2[i,j,k,l]/summV2[i,j,l] } } } } # fill up probability matrix for(i in 1:npop){ for(j in 1:nrelocc){ for(k in 1:ndest){ for(l in 1:nrecocc){ probmat[i,j,k,l] <- survivalmat[j,l] * r[k] * m[i,j,k,l] probmatV2[i,j,k,l] <- survivalmat[j,l] * r[k] * mV2[i,j,k,l] } } } } # rearrange probabilities into long format for(i in 1:npop){ for(j in 1:nrelocc){ for(k in 1:ndest){ probmatlong[(i-1)*nrelocc+j, ((k-1)*nrecocc+1):(k*nrecocc)] <- probmat[i,j,k,1:nrecocc] probmatlongV2[(i-1)*nrelocc+j, ((k-1)*nrecocc+1):(k*nrecocc)] <- probmatV2[i,j,k,1:nrecocc] } probmatlong[(i-1)*nrelocc+j, nrecocc*ndest+1] <- 1-sum(probmatlong[(i-1)*nrelocc+j, 1:(nrecocc*ndest)]) probmatlongV2[(i-1)*nrelocc+j, nrecocc*ndest+1] <- 1-sum(probmatlongV2[(i-1)*nrelocc+j, 1:(nrecocc*ndest)]) } } for(i in 1:(npop*nrelocc)){ recmatlong[i,1:33] <- rmultinom(1, prob=probmatlong[i,1:33],size=nreleased[i]) recmatlongV2[i,1:33] <- rmultinom(1, prob=probmatlongV2[i,1:33],size=nreleased[i]) } nrecoveriesV1[simulnr, rowres] <- sum(recmatlong[,-33]) nrecoveriesV2[simulnr, rowres] <- sum(recmatlongV2[,-33]) nrsims <- nrsims +1 if(nrsims>50) break } # close while-loop for simulation if(nrsims>50) cat("rownr", rowres, "cannot be computed due to too low sample size\n") if(nrsims<50){ dataxsimV1 <- list(recmatlong=recmatlong, nreleased=nreleased, npop=npop, nrelocc=nrelocc, ndest=ndest, nrecocc=nrecocc) dataxsimV2 <- list(recmatlong=recmatlongV2, nreleased=nreleased, npop=npop, nrelocc=nrelocc, ndest=ndest, nrecocc=nrecocc, relseason=relseason) modV1 <- jags(dataxsimV1, inits=initfunV1, parameters, "distributionmodelV1.txt", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, working.directory=bugsworkingdir) modV2 <- jags(dataxsimV2, inits=initfunV2, parameters, "distributionmodelV2.txt", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, working.directory=bugsworkingdir) # look at the bias bugsoutV1 <- modV1$BUGSoutput bugsoutV2 <- modV2$BUGSoutput estsV1[simulnr, rowres] <- bugsoutV1$mean$s estsV2[simulnr, rowres] <- bugsoutV2$mean$s lowerV1.s[simulnr, rowres] <- bugsoutV1$summary["s", "2.5%"] upperV1.s[simulnr, rowres] <- bugsoutV1$summary["s", "97.5%"] lowerV2.s[simulnr, rowres] <- bugsoutV2$summary["s", "2.5%"] upperV2.s[simulnr, rowres] <- bugsoutV2$summary["s", "97.5%"] estrV1[simulnr,rowres,] <- bugsoutV1$mean$r estrV2[simulnr,rowres,] <- bugsoutV2$mean$r lower.rV1[simulnr,rowres,] <- bugsoutV1$summary[paste("r[", 1:ndest, "]", sep=""), "2.5%"] lower.rV2[simulnr,rowres,] <- bugsoutV2$summary[paste("r[", 1:ndest, "]", sep=""), "2.5%"] upper.rV1[simulnr,rowres,] <- bugsoutV1$summary[paste("r[", 1:ndest, "]", sep=""), "97.5%"] upper.rV2[simulnr,rowres,] <- bugsoutV2$summary[paste("r[", 1:ndest, "]", sep=""), "97.5%"] estmV1[simulnr,rowres,] <- bugsoutV1$summary[2:769,"mean"] estmV2[simulnr,rowres,] <- bugsoutV2$summary[2:769,"mean"] lower.mV1[simulnr,rowres,] <- bugsoutV1$summary[2:769,"2.5%"] lower.mV2[simulnr,rowres,] <- bugsoutV2$summary[2:769,"2.5%"] upper.mV1[simulnr,rowres,] <- bugsoutV1$summary[2:769,"97.5%"] upper.mV2[simulnr,rowres,] <- bugsoutV2$summary[2:769,"97.5%"] truemV1[simulnr, rowres,] <- as.numeric(m) truemV2[simulnr, rowres,] <- as.numeric(mV2) save.image("simresults.RData") } # if-command: jags is only done when enough recoveries } # simulnr } # rowres # ----------------------------------------------------------------------------- # look at the simulation results # open the saved R-Console "simresults.RData" # ----------------------------------------------------------------------------- # mean number of recoveries per settning result$nrecV1 <- apply(nrecoveriesV1, 2, mean) result$nrecV2 <- apply(nrecoveriesV2, 2, mean) result$bias.sV1 <- apply(estsV1-s, 2, mean) result$se.sV1 <- apply(estsV1-s, 2, sd) incl.sV1 <- lowerV1.ss result$incl.sV1 <- apply(incl.sV1, 2, mean) result$bias.sV2 <- apply(estsV2-s, 2, mean) result$se.sV2 <- apply(estsV2-s, 2, sd) incl.sV2 <- lowerV2.ss result$incl.sV2 <- apply(incl.sV2, 2, mean) truer <- array(rep(r, each=nsim*nrow(result)), dim=dim(estrV1)) truermat <- matrix(ncol=dim(truer)[2], nrow=dim(truer)[1]*dim(truer)[3]) for(i in 1:dim(truer)[3]) truermat[((i-1)*dim(truer)[1]+1):(i*dim(truer)[1]),] <- truer[,,i] estrV1mat <- matrix(ncol=dim(truer)[2], nrow=dim(truer)[1]*dim(truer)[3]) for(i in 1:dim(truer)[3]) estrV1mat[((i-1)*dim(truer)[1]+1):(i*dim(truer)[1]),] <- estrV1[,,i] estrV2mat <- matrix(ncol=dim(truer)[2], nrow=dim(truer)[1]*dim(truer)[3]) for(i in 1:dim(truer)[3]) estrV2mat[((i-1)*dim(truer)[1]+1):(i*dim(truer)[1]),] <- estrV2[,,i] lower.rV1mat <- matrix(ncol=dim(truer)[2], nrow=dim(truer)[1]*dim(truer)[3]) for(i in 1:dim(truer)[3]) lower.rV1mat[((i-1)*dim(truer)[1]+1):(i*dim(truer)[1]),] <- lower.rV1[,,i] upper.rV1mat <- matrix(ncol=dim(truer)[2], nrow=dim(truer)[1]*dim(truer)[3]) for(i in 1:dim(truer)[3]) upper.rV1mat[((i-1)*dim(truer)[1]+1):(i*dim(truer)[1]),] <- upper.rV1[,,i] lower.rV2mat <- matrix(ncol=dim(truer)[2], nrow=dim(truer)[1]*dim(truer)[3]) for(i in 1:dim(truer)[3]) lower.rV2mat[((i-1)*dim(truer)[1]+1):(i*dim(truer)[1]),] <- lower.rV2[,,i] upper.rV2mat <- matrix(ncol=dim(truer)[2], nrow=dim(truer)[1]*dim(truer)[3]) for(i in 1:dim(truer)[3]) upper.rV2mat[((i-1)*dim(truer)[1]+1):(i*dim(truer)[1]),] <- upper.rV2[,,i] incl.rV1 <- lower.rV1mattruermat incl.rV2 <- lower.rV2mattruermat result$bias.rV1 <- apply(estrV1mat-truermat, 2, mean) result$se.rV1 <- apply(estrV1mat-truermat, 2, sd) result$incl.rV1 <- apply(incl.rV1, 2, mean) result$bias.rV2 <- apply(estrV2mat-truermat, 2, mean) result$se.rV2 <- apply(estrV2mat-truermat, 2, sd) result$incl.rV2 <- apply(incl.rV2, 2, mean) estmV1mat <- matrix(ncol=dim(estmV1)[2], nrow=dim(estmV1)[1]*dim(estmV1)[3]) for(i in 1:dim(estmV1)[3]) estmV1mat[((i-1)*dim(estmV1)[1]+1):(i*dim(estmV1)[1]),] <- estmV1[,,i] estmV2mat <- matrix(ncol=dim(estmV1)[2], nrow=dim(estmV1)[1]*dim(estmV1)[3]) for(i in 1:dim(estmV1)[3]) estmV2mat[((i-1)*dim(estmV1)[1]+1):(i*dim(estmV1)[1]),] <- estmV2[,,i] truemV1mat <- matrix(ncol=dim(estmV1)[2], nrow=dim(estmV1)[1]*dim(estmV1)[3]) for(i in 1:dim(estmV1)[3]) truemV1mat[((i-1)*dim(estmV1)[1]+1):(i*dim(estmV1)[1]),] <- truemV1[,,i] truemV2mat <- matrix(ncol=dim(estmV1)[2], nrow=dim(estmV1)[1]*dim(estmV1)[3]) for(i in 1:dim(estmV1)[3]) truemV2mat[((i-1)*dim(estmV1)[1]+1):(i*dim(estmV1)[1]),] <- truemV2[,,i] lower.mV1mat <- matrix(ncol=dim(estmV1)[2], nrow=dim(estmV1)[1]*dim(estmV1)[3]) for(i in 1:dim(estmV1)[3]) lower.mV1mat[((i-1)*dim(estmV1)[1]+1):(i*dim(estmV1)[1]),] <- lower.mV1[,,i] upper.mV1mat <- matrix(ncol=dim(estmV1)[2], nrow=dim(estmV1)[1]*dim(estmV1)[3]) for(i in 1:dim(estmV1)[3]) upper.mV1mat[((i-1)*dim(estmV1)[1]+1):(i*dim(estmV1)[1]),] <- upper.mV1[,,i] lower.mV2mat <- matrix(ncol=dim(estmV1)[2], nrow=dim(estmV1)[1]*dim(estmV1)[3]) for(i in 1:dim(estmV1)[3]) lower.mV2mat[((i-1)*dim(estmV1)[1]+1):(i*dim(estmV1)[1]),] <- lower.mV2[,,i] upper.mV2mat <- matrix(ncol=dim(estmV1)[2], nrow=dim(estmV1)[1]*dim(estmV1)[3]) for(i in 1:dim(estmV1)[3]) upper.mV2mat[((i-1)*dim(estmV1)[1]+1):(i*dim(estmV1)[1]),] <- upper.mV2[,,i] incl.mV1 <- lower.mV1mat<=truemV1mat&upper.mV1mat>=truemV1mat incl.mV2 <- lower.mV2mat<=truemV2mat&upper.mV2mat>=truemV2mat incl.mV1[truemV1mat==0] <- NA # the values fixed to 0 are not estimated! incl.mV2[truemV2mat==0] <- NA result$bias.mV1 <- apply(estmV1mat- truemV1mat, 2, mean) result$se.mV1 <- apply(estmV1mat- truemV1mat, 2, sd) result$incl.mV1 <- apply(incl.mV1, 2, function(x)sum(x, na.rm=TRUE)/sum(!is.na(x))) result$bias.mV2 <- apply(estmV2mat- truemV2mat, 2, mean) result$se.mV2 <- apply(estmV2mat- truemV2mat, 2, sd) result$incl.mV2 <- apply(incl.mV2, 2, function(x)sum(x, na.rm=TRUE)/sum(!is.na(x))) names(result) round(result[,c("n", "n0", "nrecV1", "nrecV2", "bias.sV1", "se.sV1", "incl.sV1", "bias.sV2", "se.sV2", "incl.sV2")], 3) round(result[,c("n", "n0", "bias.rV1", "se.rV1", "incl.rV1", "bias.rV2", "se.rV2", "incl.rV2")], 6) round(result[,c("n", "n0", "bias.mV1", "se.mV1","incl.mV1", "bias.mV2", "se.mV2", "incl.mV2")], 6) summary(result) round(result[,c("n", "n0","nrecV2", "bias.sV2","se.sV2","incl.sV2", "bias.rV2", "se.rV2", "incl.rV2", "bias.mV2", "se.mV2", "incl.mV2")], 6) # survival probability s # rows = number of simulations, columns = simulation settning xvalues <- 1:nsim n <- unique(result$n) n0 <- unique(result$n0) range(lowerV1.s) par(mfrow=c(3,3), mar=c(2,4,3,1), oma=c(0,0,3,0)) for(i in n){ for(j in n0){ index <- result$n==i & result$n0==j plot(xvalues, estsV1[, index], ylim=c(0.68,1), ylab="Survival", xlab="", main=paste("n=", i, ", n0=", j)) segments(xvalues, lowerV1.s[, index], xvalues, upperV1.s[,index]) abline(h=s, lty=3) } } mtext("Version 1", outer=TRUE, side=3) par(mfrow=c(3,3), mar=c(2,4,3,1), oma=c(0,0,3,0)) for(i in n){ for(j in n0){ index <- result$n==i & result$n0==j plot(xvalues, estsV2[, index], ylim=c(0.68,1), ylab="Survival", xlab="", main=paste("n=", i, ", n0=", j)) segments(xvalues, lowerV2.s[, index], xvalues, upperV2.s[,index]) abline(h=s, lty=3) } } mtext("Version 2", outer=TRUE, side=3) par(mfrow=c(1,1)) plot(nrecoveriesV1, estsV1, ylim=c(0.68, 1), ylab="survival estimate", xlab="number of recoveries", main="version 1", cex.lab=1.4, cex.axis=1.1, las=1) segments(nrecoveriesV1, lowerV1.s,nrecoveriesV1, upperV1.s) abline(h=s, lty=3) par(mfrow=c(1,1)) plot(nrecoveriesV2, estsV2, ylim=c(0.68, 1), ylab="survival estimate", xlab="number of recoveries", main="version 2", cex.lab=1.4, cex.axis=1.1, las=1) segments(nrecoveriesV2, lowerV2.s,nrecoveriesV2, upperV2.s) abline(h=s, lty=3) # recovery probability par(mfrow=c(3,3), mar=c(2,4,3,1), oma=c(0,0,3,0)) for(i in n){ for(j in n0){ index <- result$n==i & result$n0==j plot(truer[, index,], estrV1[, index,], ylim=c(0,0.005), ylab="Estimated recovery probability", xlab="True recovery probability", main=paste("n=", i, ", n0=", j)) segments(truer[, index,], lower.rV1[, index,], truer[, index,], upper.rV1[,index,]) abline(0,1, lty=3) } } mtext("Version 1", outer=TRUE, side=3) par(mfrow=c(3,3), mar=c(2,4,3,1), oma=c(0,0,3,0)) for(i in n){ for(j in n0){ index <- result$n==i & result$n0==j plot(truer[, index,], estrV2[, index,], ylim=c(0,0.005), ylab="Recovery probability", xlab="", main=paste("n=", i, ", n0=", j)) segments(truer[, index,], lower.rV2[, index,], truer[, index,], upper.rV2[,index,]) abline(0,1, lty=3) } } mtext("Version 2", outer=TRUE, side=3) # bird distribution par(mfrow=c(3,3), mar=c(4,4,3,1), oma=c(0,0,3,0)) for(i in n){ for(j in n0){ index <- result$n==i & result$n0==j corvalue <- round(cor(as.numeric(truemV1[, index,]), as.numeric(estmV1[, index,])), 2) plot(truemV1[, index,], estmV1[, index,], ylim=c(0,1), ylab="Estimated bird distribution", xlab="True bird distribution", main=paste("n=", i, ", n0=",j, ", r=", corvalue, sep=""), col=rgb(0,0,0,0.2)) #segments(truemV1[, index,], lower.mV1[, index,], truemV1[, index,], upper.mV1[,index,], col=rgb(0,0,0,0.2)) abline(0,1, lty=3, col="orange") } } mtext("Version 1", outer=TRUE, side=3) plot(truemV1[1, 1, 200:300], estmV1[1, 1, 200:300], ylim=c(0,1), ylab="Estimated bird distribution", xlab="True bird distribution", main=paste("n=", i, ", n0=",j, ", r=", corvalue, sep="")) segments(truemV1[1, 1, 200:300], lower.mV1[1, 1, 200:300], truemV1[1, 1, 200:300], upper.mV1[1, 1, 200:300]) abline(0,1, lty=3, col="orange") cor(truemV1, estmV1) #---------------------------------------------- # figure for appendix 2 n <-c(50000, 10000, 1000) n0 <- c(0,2,5) par(mfrow=c(3,3), mar=c(4,4,3,1), oma=c(0,0,3,0)) for(i in n){ for(j in n0){ index <- result$n==i & result$n0==j corvalue <- round(cor(as.numeric(truemV2[, index,]), as.numeric(estmV2[, index,])), 2) plot(truemV2[, index,], estmV2[, index,], ylim=c(0,1), ylab="Estimated bird distribution", xlab="True bird distribution", main=paste("n=", i, ", n0=",j, ", r=", corvalue, sep=""), col=rgb(0,0,0,0.2)) #segments(truemV2[, index,], lower.mV2[, index,], truemV2[, index,], upper.mV2[,index,], col=rgb(0,0,0,0.2)) abline(0,1, lty=3, col="orange") } } # end of fig for supplement #------------------------------------------------------------------------------------------------------------------ cor(truemV2, estmV2) biasV1 <- estmV1-truemV1 plot(nrecoveriesV1, biasV1[,,1], col=rgb(0,0,0,0.2)) for(i in 2:768) points(nrecoveriesV1, biasV1[,,i], col=rgb(0,0,0,0.2))