################################################### ### chunk number 1: ################################################### #To produce the .tex file use #Sweave("RY2008EcologyCommentSupplement.Rnw") #To produce a file with just the R commands use #Stangle("RY2008EcologyCommentSupplement.Rnw") ################################################### ### chunk number 2: opt ################################################### #this is used do that code does not go out of page in dvi and pdf options(width=70) ################################################### ### chunk number 3: fdt ################################################### library(xtable) library(R2WinBUGS) ################################################### ### chunk number 4: create.data ################################################### create.data<-function(seed=123,min=7,max=23,N=200,Times=5,p=0.25,delta=3,nzeroes=250,sigma=1){ #This function creates and packs appropriately the required data to feed into WinBUGS for replicating #the data analysis and simulation exercise from Royle & Young 2008. #Inputs: #seed a random seed, used so that the same results can be obtained again #min the lower x,y coordinates of the wider area S that defines the prior distribution for the animal's home range centers (default value 7 from RY2008) #max the upper x,y coordinates of the wider area S that defines the prior distribution for the animal's home range centers (default value 23 from RY2008) #N the number of home range centers within S (default value 200 from RY2008) #Times the number of sampling occasions (5 by default) #p the detection probability, conditional on an animal's location being within D (D is the sampling 10 by 10 inner square) #nzeroes the number of all zero's capture histories to be added to the data (for data augmentation) #sigma the movement parameter (a Gaussian standard deviation) set.seed(seed) #defining some required quatities Xl<-Yl<-min; Xu<-Yu<-max;sigma1<-sigma2<-sigma #generating the home range centers S1i<-runif(N,Xl,Xu); S2i<-runif(N,Yl,Yu); S<-data.frame(S1i,S2i) #holders for the animals positions over sampling ocasions U1x<-matrix(NA,nrow=N,ncol=Times);U2y<-matrix(NA,nrow=N,ncol=Times) #generating those positions for(j in 1:Times) { U1x[,j]<-rnorm(N,S1i,sigma1) U2y[,j]<-rnorm(N,S2i,sigma2) } #now get capture histories #check who's inside the plot D flag1<- U1x > Xl+(delta*sigma); flag2<- U1x < Xu-(delta*sigma); flag3<- U2y > Yl+(delta*sigma); flag4<- U2y < Yu-(delta*sigma) inplot<-flag1*flag2*flag3*flag4 #check who gets detected in each ocasion detected<-matrix(rbinom(n=N*Times, size=1, prob=p),nrow=N,ncol=Times) detected<-detected*inplot #number of individuals detected nind<-sum(apply(detected,1,max)) #the detected locations U1<-U1x*detected;U2<-U2y*detected #removing the capture histories for animals which were not observed at all capturehist<-detected[apply(U1,1,sum)>0,] #now keeping only locations for the animals observed at least once U1<-U1[apply(U1,1,sum)>0,];U2<-U2[apply(U2,1,sum)>0,] #and recoding the locations for which the animals were not detected as missing values U1<-ifelse(U1==0,NA,U1);U2<-ifelse(U1==0,NA,U2) #and the missing values for capture locations corresponding to the data augmented histories missing<-matrix(NA,nrow=nzeroes,ncol=Times) U1<-rbind(U1,missing);U2<-rbind(U2,missing) #adding the all zeroes capture histories zeros<-matrix(0,nrow=nzeroes,ncol=Times) Y<-rbind(capturehist,zeros) #the required data, as defined in Royle & Young supplementary material (R 'list' format) #data<-list(Xl=Xl,Xu=Xu,Yl=Yl,Yu=Yu,delta=delta,nind=nind,nzeroes=nzeroes,T=Times,Y=Y,U1=U1,U2=U2) data<-list(Xl=Xl,Xu=Xu,Yl=Yl,Yu=Yu,delta=delta,nind=nind,nzeroes=nzeroes,T=Times,Y=Y,U1=U1,U2=U2,sigma=sigma) return(data)} ################################################### ### chunk number 5: MCMC.inits ################################################### #creating inital values for WinBUGS MCMC.inits<-function(chains,nind,nzeroes,Xl,Xu,Yl,Yu){ #creating initial values inits<-vector("list", chains) for(j in 1:chains) { inits[[j]]<-list(psi=runif(1,0,1),p=runif(1,0,1),sigma1=runif(1,0,10),sigma2=runif(1,0,10),z=rbinom(nind+nzeroes,1,0.5),s1=runif(nind+nzeroes,Xl,Xu),s2=runif(nind+nzeroes,Yl,Yu))} return(inits)} ################################################### ### chunk number 6: mode ################################################### post.mode<-function (x, adjust = 1, ...) { find.mode <- function(x, adjust, ...) { dx <- density(x, adjust = adjust, ...) dx$x[which.max(dx$y)] } apply(as.matrix(x), 2, find.mode, adjust = adjust, ...) } ################################################### ### chunk number 7: run.sims ################################################### run.simRY.M2<-function(B,burn,iter,thin,chains,monitor=c("psi","p","sigma1","sigma2","N","ND"),bug.file,N,sigma,min=7,max=23,delta=3,bugs.directory="c:/x/WinBUGS14/",nzeroes,p=0.25,Times=5,name="OUT"){ #this function returns an array, which 3rd dimension refers to the simulation number, second to the parameters monitored, and first the quantities recorded from each parameter. #INPUTS: # B - number of simulations to run for a given scenario # burn,iter,thin,chains - the MCMC details, namely number of burn in iterations, iterations for inference, thining factor and number of chains # monitor - the parameters to monitor, currently and by default, "psi","p","sigma1","sigma2","N","ND" (they appear in this order in the created object, followed by the Deviance) # bug.file - the name of the model file for WinBUGS # N - The number of animals in area S (the wider area, or equivalently the way the prior for home ranges is defined) # sigma - The movement parameter # min,max,delta - the limits of the area S and the delta that defines the inner square as a 10 by 10 meter square # bugs.directory - where can R find WinBUGS # nzeroes - the augmented all 0's capture histories #Returns an array with dimensions [7, 7 or 9,B] #Within each level of the third dimention, i.e. for each simulation, columns represent the mean, sd, 2.5%, 25%, 50%, 75%, 97.5%, median for each variable #+ (3 values in the first 3 rows of the last column, namely: 1. the random seed used to create the data; 2. the number of animals captured; 3. the number of nzeroes used #the rows represent each one of the statistics being monitored, i.e. "psi","p","sigma1","sigma2","N","ND" (they appear in this order in the created object, followed by the Deviance) #set a seed for data creation routine semente<-round(runif(B,1,1000000)) #if there are not at least 2 chains, no convergence results are reported, so no R.hat nor n.eff stat<-ifelse(chains==1,8,10) #get stats names stat.names<-get.summary.names(chains) #create an object with correct dimension and appropriate names to store the summary stats ptms<-c("psi","p","sigma1","sigma2","N","ND","Deviance") storage<-array(data = NA, dim = c(7,stat+2,B),dimnames=list(parameter=ptms,statistic=stat.names,simulation=1:B)) #create an object with correct dimension and appropriate names to store the MCMC outpt #note that array's first dimention is the simulation number, second is the number of iterations, thirs the number of chains and fourth the number of parameters being monitored ci<-(iter-burn)/thin #number of iterations for inference MCMCout.names<-list(sim=paste("sim",1:B,sep=""),MCMCi=1:ci,chain=paste("C",1:chains,sep=""),pars=c(monitor,"Deviance")) MCMCout<-array(data = NA, dim = c(B,ci ,chains,length(monitor)+1), dimnames = MCMCout.names) #create an object with correct dimension and appropriate names to store the data DATAout<-vector("list", B) #for timing the procedure start.time<-Sys.time() #for each bootstrap iteration for(i in 1:B) { #print when it started cat("Bootstrap rep",i,"started at",date(),"\n") #create the simulated data ds<-create.data(semente[i],N=N,sigma=sigma,nzeroes=nzeroes,min=min,max=max,delta=delta,p=p,Times=Times) #creating initial values inits<-MCMC.inits(chains,nind=ds$nind,nzeroes=ds$nzeroes,Xl=ds$Xl,Xu=ds$Xu,Yl=ds$Yl,Yu=ds$Yu) #calling winbugs and run the analysis rs<- bugs(ds, inits, parameters.to.save=monitor, model.file=bug.file,n.chains = chains, n.iter = iter, n.burnin = burn,n.thin=thin,bugs.directory=bugs.directory) #storing MCMC results MCMCout[i,,,]<-rs$sims.array #storing the data DATAout[[i]]<-ds #storing summary results storage<-get.summary.results(rs,storage,stat,i,N,sigma,semente,ds) #save what's gone on so far save(storage,file="storage.Rdata") #save what's gone on so far #save(MCMCout,file="MCMCout.Rdata") #save what's gone on so far #save(dataOUT,file="dataOUT.Rdata") finish.time<-Sys.time() elapsed.time.per.rep<-as.numeric(difftime(finish.time,start.time,units="secs"))/i cat("Finished rep ",i,". Estimated completion time ",as.character(Sys.time()+elapsed.time.per.rep*(B-i)),"\n",sep="") } #saving data and MCMC output as adequately named objects assign(paste("DATAout",name, sep = ""),DATAout) save(list=paste("DATAout",name, sep = ""),file=paste("DATAout",name,".Rdata",sep="")) assign(paste("MCMCout",name, sep = ""),MCMCout) save(list=paste("MCMCout",name, sep = ""),file=paste("MCMCout",name,".Rdata",sep="")) return(storage)} ################################################### ### chunk number 8: more.bits.needed ################################################### get.summary.names<-function(chains){ #This function returns the names of the summary object #which are a function of the number of chains if(chains==1) { stat.names<-c("mean","sd","2.5%","25%","50%","75%","97.5%","median","seed","mode") } else { stat.names<-c("mean","sd","2.5%","25%","50%","75%","97.5%","median","Rhat","n.eff","seed","mode") } return(stat.names) } get.summary.results<-function(rs,storage,stat,i,N,sigma,semente,ds){ #gets summary information regarding current iteration #mean, sd and quantiles for each parameter storage[,1:(stat-1),i]<-rs$summary[(nrow(rs$summary)-6):nrow(rs$summary),] #median storage[,stat,i]<-as.numeric(rs$median) #the mode - requires bespoke function for(j in 1:7) { storage[j,stat+2,i]<-post.mode(as.numeric(rs$sims.array[,,j])) } #the random seed used storage[1,stat+1,i]<-semente[i] #the nind used storage[2,stat+1,i]<-ds$nind #the nzeroes effectivelly used storage[3,stat+1,i]<-ds$nzeroes #the N used storage[4,stat+1,i]<-N #the sigma storage[5,stat+1,i]<-sigma return(storage) } ################################################### ### chunk number 9: inputs ################################################### ################################################### ### Required DATA for sims ################################################### #the densities Ds<-c(0.234,0.352,0.469,0.586,0.781) #the sigmas sigmas<-c(1,2,4) #the required deltas are defined as 3 times the sigma (Royle mentions 2 times SE in his Table 1 notes ???) deltas<-3*sigmas #the inner sampling square side square.side<-10 #inner sampling square area chi<-square.side^2 #the larger conceptual S area Ss<-(square.side+2*deltas)^2 #these are how many animals need to be generated within a S of sizes Ss (the area size depends on sigma) NSs<-matrix(NA,nrow=3,ncol=5) for(i in 1:3) { NSs[i,]<-Ds*Ss[i]} rownames(NSs)<-paste("sigma =",sigmas) colnames(NSs)<-paste("Ds =",Ds) #these are the number of animals within the inner sampling square inChi<-matrix(NA,nrow=3,ncol=5) for(i in 1:3) inChi[i,]<-Ds*100 #these are the proposed number of extra zeroes to be added; this was chose to match the 250 used by RY2008 for the initial example extrazeroes<-NSs*1.25 #animals which should be inside S but outside the sampling frame inSnotinChi<-NSs-inChi #the mean probability of detection for individuals in the simulated populations #(inSnotinChi*0+inChi*0.25)/NSs #for simulations we need integer values roundNs<-round(NSs) ################################################### ### chunk number 10: sims eval=FALSE ################################################### ## ################################################### ## ### chunk number 34: ReplicaOfRYtable2 eval=FALSE ## ################################################### ## #required strings for making code easier to read ## #required strings for making code easier to read ## ptm<-c("psi","p","sigma1","sigma2","N","ND") ## BugsDir<-"C:/x/WinBUGS14" ## fs<-"WBMF.txt" ## #simulation parameters ## B<-100;thin<-1;chains<-1;iter<-16000;burn<-6000 ## #The Simulations - note this takes a long time to run ## #----------------------------------------------------------------- D=0.234 ## RYD0.234sigmat1<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=ptm,bug.file=fs,N=roundNs[1,1],sigma=1,bugs.directory=BugsDir,nzeroes=roundNs[1,1]*3,name="RYD0.234sigmat1") ## save(RYD0.234sigmat1, file = "RYD0.234sigmat1.Rdata") ## RYD0.234sigmat2<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=ptm,bug.file=fs,N=roundNs[2,1],sigma=2,min=4,max=26,delta=3,bugs.directory=BugsDir,nzeroes=roundNs[2,1]*3,name="RYD0.234sigmat2") ## save(RYD0.234sigmat2, file = "RYD0.234sigmat2.Rdata") ## RYD0.234sigmat4<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=ptm,bug.file=fs,N=roundNs[3,1],sigma=4,min=1,max=35,delta=3,bugs.directory=BugsDir,nzeroes=roundNs[3,1]*3,name="RYD0.234sigmat4") ## save(RYD0.234sigmat4, file = "RYD0.234sigmat4.Rdata") ## #----------------------------------------------------------------- D=0.352 ## RYD0.352sigmat1<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=ptm,bug.file=fs,N=roundNs[1,2],sigma=1,bugs.directory=BugsDir,nzeroes=roundNs[1,2]*3,name="RYD0.352sigmat1") ## save(RYD0.352sigmat1, file = "RYD0.352sigmat1.Rdata") ## RYD0.352sigmat2<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=ptm,bug.file=fs,N=roundNs[2,2],sigma=2,min=4,max=26,delta=3,bugs.directory=BugsDir,nzeroes=roundNs[2,2]*3,name="RYD0.352sigmat2") ## save(RYD0.352sigmat2, file = "RYD0.352sigmat2.Rdata") ## RYD0.352sigmat4<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=ptm,bug.file=fs,N=roundNs[3,2],sigma=4,min=1,max=35,delta=3,bugs.directory=BugsDir,nzeroes=roundNs[3,2]*3,name="RYD0.352sigmat4") ## save(RYD0.352sigmat4, file = "RYD0.352sigmat4.Rdata") ## #----------------------------------------------------------------- D=0.469 ## RYD0.469sigmat1<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=ptm,bug.file=fs,N=roundNs[1,3],sigma=1,bugs.directory=BugsDir,nzeroes=roundNs[1,3]*3,name="RYD0.469sigmat1") ## save(RYD0.469sigmat1, file = "RYD0.469sigmat1.Rdata") ## RYD0.469sigmat2<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=ptm,bug.file=fs,N=roundNs[2,3],sigma=2,min=4,max=26,delta=3,bugs.directory=BugsDir,nzeroes=roundNs[2,3]*3,name="RYD0.469sigmat2") ## save(RYD0.469sigmat2, file = "RYD0.469sigmat2.Rdata") ## RYD0.469sigmat4<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=ptm,bug.file=fs,N=roundNs[3,3],sigma=4,min=1,max=35,delta=3,bugs.directory=BugsDir,nzeroes=roundNs[3,3]*3,name="RYD0.469sigmat4") ## save(RYD0.469sigmat4, file = "RYD0.469sigmat4.Rdata") ## #----------------------------------------------------------------- D=0.586 ## RYD0.586sigmat1<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=ptm,bug.file=fs,N=roundNs[1,4],sigma=1,bugs.directory=BugsDir,nzeroes=roundNs[1,4]*3,name="RYD0.586sigmat1") ## save(RYD0.586sigmat1, file = "RYD0.586sigmat1.Rdata") ## RYD0.586sigmat2<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=ptm,bug.file=fs,N=roundNs[2,4],sigma=2,min=4,max=26,delta=3,bugs.directory=BugsDir,nzeroes=roundNs[2,4]*3,name="RYD0.586sigmat2")# ## save(RYD0.586sigmat2, file = "RYD0.586sigmat2.Rdata") ## RYD0.586sigmat4<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=ptm,bug.file=fs,N=roundNs[3,4],sigma=4,min=1,max=35,delta=3,bugs.directory=BugsDir,nzeroes=roundNs[3,4]*3,name="RYD0.586sigmat4") ## save(RYD0.586sigmat4, file = "RYD0.586sigmat4.Rdata") ## #----------------------------------------------------------------- D=0.781 ## RYD0.781sigmat1<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=ptm,bug.file=fs,N=roundNs[1,5],sigma=1,bugs.directory=BugsDir,nzeroes=roundNs[1,5]*3,name="RYD0.781sigmat1") ## save(RYD0.781sigmat1, file = "RYD0.781sigmat1.Rdata") ## RYD0.781sigmat2<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=ptm,bug.file=fs,N=roundNs[2,5],sigma=2,min=4,max=26,delta=3,bugs.directory=BugsDir,nzeroes=3*roundNs[2,5],name="RYD0.781sigmat2") ## save(RYD0.781sigmat2, file = "RYD0.781sigmat2.Rdata") ## RYD0.781sigmat4<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=ptm,bug.file=fs,N=roundNs[3,5],sigma=4,min=1,max=35,delta=3,bugs.directory=BugsDir,nzeroes=3*roundNs[3,5],name="RYD0.781sigmat4") ## save(RYD0.781sigmat4, file = "RYD0.781sigmat4.Rdata") ## ## #--------------------------------------------------------------------------------------------------------------------------------------------------------------- ################################################### ### chunk number 11: sigma4D0.781T10 eval=FALSE ################################################### ## B<-100;thin<-1;chains<-1;iter<-20000;burn<-5000 ## RYD0.781sigmat4.T10.da3N<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=c("psi","p","sigma1","sigma2","N","ND"),bug.file=fs, ## N=roundNs[3,5],sigma=4,min=1+20,max=35+20,delta=3,bugs.directory=BugsDir,nzeroes=3*roundNs[3,5],Times=10,name="RYD0.781sigmat4.T10.da3N") ## save(RYD0.781sigmat4.T10.da3N, file = "RYD0.781sigmat4.T10.da3N.Rdata") ################################################### ### chunk number 12: lizard.sim eval=FALSE ################################################### ## B<-100;thin<-1;chains<-1;iter<-16000;burn<-6000 ## lizard<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=c("psi","p","sigma1","sigma2","N","ND"), ## bug.file="WBMF.txt",N=155,sigma=0.15,min=2.4,max=6.6,delta=4,bugs.directory=BugsDir,nzeroes=200,p=0.122,Times=14,name="lizard") ## save(lizard, file = "lizard.Rdata") ################################################### ### chunk number 13: brokentillworks.sim eval=FALSE ################################################### ## B<-100;thin<-1;chains<-1;iter<-16000;burn<-6000 ## #ORIGINAL SCENARIO ## original<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=ptm,bug.file=fs,N=roundNs[1,1],sigma=1,bugs.directory=BugsDir,nzeroes=roundNs[1,1]*3,name="original") ## save(original, file = "original.Rdata") ## #B<-100;thin<-1;chains<-1;iter<-20000;burn<-5000 ## #-------------- Less sigma ## sigma0.5p0.25D0.234T5<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=ptm,bug.file=fs,N=roundNs[1,1],sigma=0.5,bugs.directory=BugsDir,nzeroes=roundNs[1,1]*3,name="sigma0.5p0.25D0.234T5") ## save(sigma0.5p0.25D0.234T5, file = "sigma0.5p0.25D0.234T5.Rdata") ## #-------------- MORE p ## sigma1p0.5D0.234T5<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=ptm,bug.file=fs,N=roundNs[1,1],sigma=1,bugs.directory=BugsDir,nzeroes=roundNs[1,1]*3,p=0.5,name="sigma1p0.5D0.234T5") ## save(sigma1p0.5D0.234T5, file = "sigma1p0.5D0.234T5.Rdata") ## #-------------- MORE Times ## sigma1p0.25D0.234T10<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=ptm,bug.file=fs,N=roundNs[1,1],sigma=1,bugs.directory=BugsDir,nzeroes=roundNs[1,1]*3,p=0.25,Times=10,name="sigma1p0.25D0.234T10") ## save(sigma1p0.5D0.234T10, file = "sigma1p0.25D0.234T10.Rdata") ## #-------------- Less sigma MORE p ## sigma0.5p0.5D0.234T5<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=ptm,bug.file=fs,N=roundNs[1,1],sigma=0.5,bugs.directory=BugsDir,nzeroes=roundNs[1,1]*3,p=0.5,name="sigma0.5p0.5D0.234T5") ## save(sigma0.5p0.5D0.234T5, file = "sigma0.5p0.5D0.234T5.Rdata") ## #-------------- MORE p MORE Times ## sigma1p0.5D0.234T10<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=ptm,bug.file=fs,N=roundNs[1,1],sigma=1,bugs.directory=BugsDir,nzeroes=roundNs[1,1]*3,p=0.5,Times=10,name="sigma1p0.5D0.234T10") ## save(sigma1p0.5D0.234T10, file = "sigma1p0.5D0.234T10.Rdata") ## #-------------- Less sigma MORE Times ## sigma0.5p0.25D0.234T10<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=ptm,bug.file=fs,N=roundNs[1,1],sigma=0.5,bugs.directory=BugsDir,nzeroes=roundNs[1,1]*3,p=0.25,Times=10,name="sigma0.5p0.25D0.234T10") ## save(sigma0.5p0.5D0.234T10, file = "sigma0.5p0.25D0.234T10.Rdata") ## #-------------- Less sigma MORE p MORE Times ## sigma0.5p0.5D0.234T10<-run.simRY.M2(B=B,burn=burn,iter=iter,thin=thin,chains=chains,monitor=ptm,bug.file=fs,N=roundNs[1,1],sigma=0.5,bugs.directory=BugsDir,nzeroes=roundNs[1,1]*3,p=0.5,Times=10,name="sigma0.5p0.5D0.234T10") ## save(sigma0.5p0.5D0.234T10, file = "sigma0.5p0.5D0.234T10.Rdata") ################################################### ### chunk number 14: codetable1RY2008 eval=FALSE ################################################### ## #how many times to repeat ## qtd<-1000 ## nam.r<-c("D0.234","D0.352","D0.469","D0.586","D0.781") ## nam.c<-c("Ds","Mean1","Mean2","Mean4") ## Ds<-c(0.234,0.352,0.469,0.586,0.781) ## #object to store results ## EDF<-matrix(NA,nrow=15,ncol=10,dimnames=list(c(11:15,6:10,1:5),c("sigma","Ds","0","1","2","3","4","5","En","Ep"))) ## sigmas<-c(1,2,4);mins<-c(7,4,1);maxs<-c(23,26,35) ## ct1<-1 ## for(j in 1:3) { ## for(i in 1:5) { ## EDF[ct1,1]<-sigmas[j] ## EDF[ct1,2]<-Ds[i] ## #objects for holding info temporarily ## hold<-numeric(qtd) ## holdEDF<-matrix(NA,ncol=6,nrow=qtd) ## for(k in 1:qtd) { ## currN<-roundNs[j,i] ## #create a data set ## inbet<-create.data(2*k,N=roundNs[j,i],sigma=sigmas[j],nzeroes=0,min=mins[j],max=maxs[j],delta=deltas[j]) ## #get the number of detected animals ## hold[k]<-inbet$nind ## #now getting the encounter frequency distribution ## curr.table<-table(c(rep(0,roundNs[j,i]-inbet$nind),apply(inbet$Y,1,sum)))/roundNs[j,i] ## ct2<-1 ## for(r in 0:5) { ## recaps<-as.numeric(names(curr.table)) ## present<-is.element(el=r,set=recaps) ## #if there were animals detected r times, get their frequency ## if(present) { ## holdEDF[k,r+1]<-curr.table[ct2] ## ct2<-ct2+1 ## } else { ## #if they were not, the frequency is 0 ## holdEDF[k,r+1]<-0 ## } ## } ## ## ## } ## #get E(n) ## EDF[ct1,9]<-mean(hold) ## #get E(p) ## EDF[ct1,10]<-(0.25*10^2+0*(Ss[j]-10^2))/Ss[j] ## EDF[ct1,3:8]<-apply(holdEDF,2,mean) ## ct1<-ct1+1 ## }} ## #sort rows to match RY2008 exactly ## tab1RY08<-EDF[c(11:15,6:10,1:5),] ## save(tab1RY08,file="tab1RY08.Rdata") ## rm(tab1RY08) ################################################### ### chunk number 15: ################################################### load("tab1RY08.Rdata") capt<-"Replica of table 1 in Royle and Young (2008)" t1RY08toprint<-xtable(tab1RY08,digits=c(0,0,rep(3,7),0,3), caption=capt, label="tab1RY2008") #print(t1RY08toprint,floating.environment='sidewaystable') ################################################### ### chunk number 16: cap.hist.by.scenario eval=FALSE ################################################### ## cap.hist.by.scenario<-matrix(NA,nrow=15,ncol=7,dimnames=list(c(1:15),c("sigma","Ds","1","2","3","4","5"))) ## cap.hist.by.scenario[,1:2]<-tab1RY08[,1:2] ## for(j in 1:15) { ## cap.hist.by.scenario[j,3:7]<-tab1RY08[j,4:8]/sum(tab1RY08[j,4:8])*tab1RY08[j,9] ## } ## cap.hist.by.scenario[,c(1,3:7)]<-round(cap.hist.by.scenario[,c(1,3:7)]) ################################################### ### chunk number 17: eval=FALSE ################################################### ## capt<-"Our table 1" ## cap.hist.by.scenariotoprint<-xtable(cap.hist.by.scenario,digits=c(0,0,3,0,0,0,0,0), ## caption=capt, ## label="cap.hist.by.scenario") ## print(cap.hist.by.scenariotoprint) ################################################### ### chunk number 18: table2 eval=FALSE ################################################### ## #rebuilding Royle's table 2 ## nam.r<-c("D0.234","D0.352","D0.469","D0.586","D0.781") ## nam.c<-c("Ds","Mean1","Med1","Sd1","IC1","Mean2","Med2","Sd2","IC2","Mean4","Med4","Sd4","IC4") ## tRY.tab2<-matrix(NA,nrow=5,ncol=13,dimnames=list(nam.r,nam.c)) ## tRY.tab2[,1]<-c(0.234,0.352,0.469,0.586,0.781) ## #this file path needs changing depending on the user's system ## file.path<-"table2/burn6000iter16000thin1/" ## roundNs<-round(NSs) ## sigmas<-c(1,2,4) ## for(i in 1:5) { ## for(j in 1:3) { ## #get the current true number of animals in S ## currN<-roundNs[j,i] ## load(paste(file.path,"RY",nam.r[i],"sigmat",sigmas[j],".RData",sep="")) ## curr.data<-eval(parse(text=paste("RY",nam.r[i],"sigmat",sigmas[j],sep=""))) ## tRY.tab2[i,2+(j-1)*4]<-mean(curr.data[5,1,]/Ss[j]) ## tRY.tab2[i,3+(j-1)*4]<-median(curr.data[5,1,]/Ss[j]) ## tRY.tab2[i,4+(j-1)*4]<-sd(curr.data[5,1,]/Ss[j]) ## tRY.tab2[i,5+(j-1)*4]<-sum(curr.data[5,3,]currN)/100 ## }} ## rm(nam.r,nam.c) ################################################### ### chunk number 19: eval=FALSE ################################################### ## capt<-"Our table 2: Replica of table 2 in Royle and Young (2008)" ## tRYtab<-xtable(tRY.tab2,digits=c(2,rep(3,4),2,rep(3,3),2,rep(3,3),2), ## caption=capt, ## label="tRYsims") ## print(tRYtab,floating.environment='sidewaystable') ################################################### ### chunk number 20: table3 eval=FALSE ################################################### ## #rebuilding Royle's table 2 -but using the mean of the mean, median and mode ## nam.r<-c("D0.234","D0.352","D0.469","D0.586","D0.781") ## nam.c<-c("Ds","Mean1","Med1","Mode1","IC1","Mean2","Med2","Mode2","IC2","Mean4","Med4","Mode4","IC4") ## tRY.tab2<-matrix(NA,nrow=5,ncol=13,dimnames=list(nam.r,nam.c)) ## tRY.tab2[,1]<-c(0.234,0.352,0.469,0.586,0.781) ## roundNs<-round(NSs) ## sigmas<-c(1,2,4) ## for(i in 1:5) { ## for(j in 1:3) { ## #get the current true number of animals in S ## currN<-roundNs[j,i] ## load(paste(file.path,"RY",nam.r[i],"sigmat",sigmas[j],".RData",sep="")) ## curr.data<-eval(parse(text=paste("RY",nam.r[i],"sigmat",sigmas[j],sep=""))) ## tRY.tab2[i,2+(j-1)*4]<-mean(curr.data[5,1,]/Ss[j]) ## tRY.tab2[i,3+(j-1)*4]<-mean(curr.data[5,8,]/Ss[j]) ## tRY.tab2[i,4+(j-1)*4]<-mean(curr.data[5,10,]/Ss[j]) ## tRY.tab2[i,5+(j-1)*4]<-sum(curr.data[5,3,]currN)/100 ## }} ## rm(nam.r,nam.c) ################################################### ### chunk number 21: eval=FALSE ################################################### ## capt<-"Our table 3" ## tRYtab2<-xtable(tRY.tab2,digits=c(2,rep(3,4),2,rep(3,3),2,rep(3,3),2), ## caption=capt, ## label="tRYsims2") ## print(tRYtab2,floating.environment='sidewaystable') ################################################### ### chunk number 22: table4 eval=FALSE ################################################### ## #building table 3 ## nam.r<-c("Original","sigma","p","T","sigma p","sigma T","p T","sigma p T") ## nam.c<-c("sigma","p","T","Mean","% Bias","Median","% Bias","Mode","% Bias") ## btw.sim.res<-matrix(NA,nrow=8,ncol=9,dimnames=list(nam.r,nam.c)) ## btw.sim.res[,1]<-c(1,0.5,1,1,0.5,0.5,1,0.5) ## btw.sim.res[,2]<-c(0.25,0.25,0.5,0.25,0.5,0.25,0.5,0.5) ## btw.sim.res[,3]<-c(5,5,5,10,5,10,10,10) ## #load(file = "brokentillworks/sigma0.5p0.25D0.234T5/sigma0.5p0.25D0.234T5.Rdata") ## #load(file = "brokentillworks/sigma1p0.5D0.234T5/sigma1p0.5D0.234T5.Rdata") ## #load(file = "brokentillworks/sigma1p0.25D0.234T10/sigma1p0.25D0.234T10.Rdata") ## #load(file = "brokentillworks/sigma0.5p0.5D0.234T5/sigma0.5p0.5D0.234T5.Rdata") ## #load(file = "brokentillworks/sigma1p0.5D0.234T10/sigma1p0.5D0.234T10.Rdata") ## #load(file = "brokentillworks/sigma0.5p0.25D0.234T10/sigma0.5p0.25D0.234T10.Rdata") ## #load(file = "brokentillworks/sigma0.5p0.5D0.234T10/sigma0.5p0.5D0.234T10.Rdata") ## true<-60 ## #Mean N estimate ## #Original ## btw.sim.res[1,4]<-round(mean(original[5,1,]),2) ## #-------------- less sigma ## btw.sim.res[2,4]<-round(mean(sigma0.5p0.25D0.234T5[5,1,]),2) ## #-------------- MORE p ## btw.sim.res[3,4]<-round(mean(sigma1p0.5D0.234T5[5,1,]),2) ## #-------------- MORE Times ## btw.sim.res[4,4]<-round(mean(sigma1p0.25D0.234T10[5,1,]),2) ## #-------------- Less sigma MORE p ## btw.sim.res[5,4]<-round(mean(sigma0.5p0.5D0.234T5[5,1,]),2) ## #-------------- Less sigma MORE Times ## btw.sim.res[6,4]<-round(mean(sigma0.5p0.25D0.234T10[5,1,]),2) ## #-------------- MORE p MORE Times ## btw.sim.res[7,4]<-round(mean(sigma1p0.5D0.234T10[5,1,]),2) ## #-------------- Less sigma MORE p MORE Times ## btw.sim.res[8,4]<-round(mean(sigma0.5p0.5D0.234T10[5,1,]),2) ## #Mean % Bias ## #Original ## btw.sim.res[1,5]<-round(100*(mean(original[5,1,])-true)/true,2) ## #-------------- less sigma ## btw.sim.res[2,5]<-round(100*(mean(sigma0.5p0.25D0.234T5[5,1,])-true)/true,2) ## #-------------- MORE p ## btw.sim.res[3,5]<-round(100*(mean(sigma1p0.5D0.234T5[5,1,])-true)/true,2) ## #-------------- MORE Times ## btw.sim.res[4,5]<-round(100*(mean(sigma1p0.25D0.234T10[5,1,])-true)/true,2) ## #-------------- Less sigma MORE p ## btw.sim.res[5,5]<-round(100*(mean(sigma0.5p0.5D0.234T5[5,1,])-true)/true,2) ## #-------------- Less sigma MORE Times ## btw.sim.res[6,5]<-round(100*(mean(sigma0.5p0.25D0.234T10[5,1,])-true)/true,2) ## #-------------- MORE p MORE Times ## btw.sim.res[7,5]<-round(100*(mean(sigma1p0.5D0.234T10[5,1,])-true)/true,2) ## #-------------- Less sigma MORE p MORE Times ## btw.sim.res[8,5]<-round(100*(mean(sigma0.5p0.5D0.234T10[5,1,])-true)/true,2) ## #MEDIAN ## #Original ## btw.sim.res[1,6]<-round(mean(original[5,8,]),2) ## #-------------- less sigma ## btw.sim.res[2,6]<-round(mean(sigma0.5p0.25D0.234T5[5,8,]),2) ## #-------------- MORE p ## btw.sim.res[3,6]<-round(mean(sigma1p0.5D0.234T5[5,8,]),2) ## #-------------- MORE Times ## btw.sim.res[4,6]<-round(mean(sigma1p0.25D0.234T10[5,8,]),2) ## #-------------- Less sigma MORE p ## btw.sim.res[5,6]<-round(mean(sigma0.5p0.5D0.234T5[5,8,]),2) ## #-------------- Less sigma MORE Times ## btw.sim.res[6,6]<-round(mean(sigma0.5p0.25D0.234T10[5,8,]),2) ## #-------------- MORE p MORE Times ## btw.sim.res[7,6]<-round(mean(sigma1p0.5D0.234T10[5,8,]),2) ## #-------------- Less sigma MORE p MORE Times ## btw.sim.res[8,6]<-round(mean(sigma0.5p0.5D0.234T10[5,8,]),2) ## #Mean % Bias ## #Original ## btw.sim.res[1,7]<-round(100*(mean(original[5,8,])-true)/true,2) ## #-------------- less sigma ## btw.sim.res[2,7]<-round(100*(mean(sigma0.5p0.25D0.234T5[5,8,])-true)/true,2) ## #-------------- MORE p ## btw.sim.res[3,7]<-round(100*(mean(sigma1p0.5D0.234T5[5,8,])-true)/true,2) ## #-------------- MORE Times ## btw.sim.res[4,7]<-round(100*(mean(sigma1p0.25D0.234T10[5,8,])-true)/true,2) ## #-------------- Less sigma MORE p ## btw.sim.res[5,7]<-round(100*(mean(sigma0.5p0.5D0.234T5[5,8,])-true)/true,2) ## #-------------- Less sigma MORE Times ## btw.sim.res[6,7]<-round(100*(mean(sigma0.5p0.25D0.234T10[5,8,])-true)/true,2) ## #-------------- MORE p MORE Times ## btw.sim.res[7,7]<-round(100*(mean(sigma1p0.5D0.234T10[5,8,])-true)/true,2) ## #-------------- Less sigma MORE p MORE Times ## btw.sim.res[8,7]<-round(100*(mean(sigma0.5p0.5D0.234T10[5,8,])-true)/true,2) ## #MODE ## #Original ## btw.sim.res[1,8]<-round(mean(original[5,10,]),2) ## #-------------- less sigma ## btw.sim.res[2,8]<-round(mean(sigma0.5p0.25D0.234T5[5,10,]),2) ## #-------------- MORE p ## btw.sim.res[3,8]<-round(mean(sigma1p0.5D0.234T5[5,10,]),2) ## #-------------- MORE Times ## btw.sim.res[4,8]<-round(mean(sigma1p0.25D0.234T10[5,10,]),2) ## #-------------- Less sigma MORE p ## btw.sim.res[5,8]<-round(mean(sigma0.5p0.5D0.234T5[5,10,]),2) ## #-------------- Less sigma MORE Times ## btw.sim.res[6,8]<-round(mean(sigma0.5p0.25D0.234T10[5,10,]),2) ## #-------------- MORE p MORE Times ## btw.sim.res[7,8]<-round(mean(sigma1p0.5D0.234T10[5,10,]),2) ## #-------------- Less sigma MORE p MORE Times ## btw.sim.res[8,8]<-round(mean(sigma0.5p0.5D0.234T10[5,10,]),2) ## #Mean % Bias ## #Original ## btw.sim.res[1,9]<-round(100*(mean(original[5,10,])-true)/true,2) ## #-------------- less sigma ## btw.sim.res[2,9]<-round(100*(mean(sigma0.5p0.25D0.234T5[5,10,])-true)/true,2) ## #-------------- MORE p ## btw.sim.res[3,9]<-round(100*(mean(sigma1p0.5D0.234T5[5,10,])-true)/true,2) ## #-------------- MORE Times ## btw.sim.res[4,9]<-round(100*(mean(sigma1p0.25D0.234T10[5,10,])-true)/true,2) ## #-------------- Less sigma MORE p ## btw.sim.res[5,9]<-round(100*(mean(sigma0.5p0.5D0.234T5[5,10,])-true)/true,2) ## #-------------- Less sigma MORE Times ## btw.sim.res[6,9]<-round(100*(mean(sigma0.5p0.25D0.234T10[5,10,])-true)/true,2) ## #-------------- MORE p MORE Times ## btw.sim.res[7,9]<-round(100*(mean(sigma1p0.5D0.234T10[5,10,])-true)/true,2) ## #-------------- Less sigma MORE p MORE Times ## btw.sim.res[8,9]<-round(100*(mean(sigma0.5p0.5D0.234T10[5,10,])-true)/true,2) ## ################################################### ### chunk number 23: eval=FALSE ################################################### ## capt<-"Our table 4" ## tbtw.sim.res<-xtable(btw.sim.res, ## caption=capt,digits=c(0,1,2,0,rep(2,6)), ## label="tbtw.sim.res") ## print(tbtw.sim.res)