#created by Mark Christie, 5/14/2012 #Please send email to redpath.christie@gmail.com for questions or wanted modifications #Please cite Carrol et al. (2012) Molec Ecol, if you use this script. #This script determines the probability of observing 10 paternities with known sample sizes and population size estimates from Australia and NZ. #Assuming equal reproductive success, this script calculates the probability of observing 10 paternities if BOTH Australian and NZ males sired offspring #Note: if you change the Australian contribution to be very small (demogrpahic closure), then then probability of observing 10 paternities is not small at all (consistent with hypothesis of demographic closure). OUT=NULL for (i in 1:10000) { #10000 equals the number of iterations Australian=2900/2 #estimated number of Australian whales Atypes=rep("a",Australian) Ntypes=rep("n",1085) #estimated number of NZ whales pop=c(Ntypes,Atypes) #combined population s1=pop[sample(1:length(pop),34,replace=TRUE)] #sample this combined population randomly, 34 times, to create fathers for all of the offspring sampled. s2=length(which(s1=="n")) #determine how many Fathers would have come from NZ s3=rep("dad",s2) #create NZ dads pop2=Ntypes[-(1:s2)] #remove NZ individuals to be replaced with known "dads" pop2=c(s3,pop2) #add known "dads" to NZ sample dads2=314 #number of males sampled from NZ population asp=sample(pop2,dads2,replace=FALSE) #take a sample of the NZ males asp2=length(which(asp=="dad")) #count how many "dads" there are OUT=rbind(OUT,asp2) #write to output } dads=OUT[,1] mean(dads) range(dads) p=length(which(dads>=10))/length(dads) #caluclate p value, here using 10 as the number of assigned paternities p