################################### # This file records the capture-recapture diversity calculation procedure ################################### library(RMark) load("baltic_squares_04_02.RData") tbins <- sort(unique(baltic_binnings$bb)) tbins_upper <- as.numeric(c("484.5", "480", "472.5", "468.7", "460", "453.5", "448.5", "442.7", "438.5", "433.4", "427.4", "423", "419.2")) tbins_lower <- as.numeric(c("492.3", "484.5", "480", "472.5", "468.7", "460", "453.5", "448.5", "442.7", "438.5", "433.4", "427.4", "423")) tbins_mid <-(tbins_lower+tbins_upper)/2 tbins_dur <- tbins_lower-tbins_upper tbins_upper = tbins_upper*-1 tbins_lower = tbins_lower*-1 tbins_mid = tbins_mid*-1 baltic_stages <- c("Pakerort", "Varangu", "Hunneberg", "Billingen","Volkhov", "Kunda", "Aseri", "Lasnamägi", "Uhaku", "Kukruse", "Haljala", "Keila", "Oandu", "Rakvere", "Nabala", "Vormsi", "Pirgu", "Porkuni","Juuru", "Raikküla", "Adavere", "Jaani", "Jaagarahu", "Rootsiküla", "Paadla", "Kuuressaare", "Kaugatuma", "Ohesaare") baltic_lower <- c(485.4, 481.5, 479.7, 475, 472.5, 468.7, 464.7, 463.0, 460.0, 458.4, 457.5, 453.5, 452.5, 451.8, 450.3, 449.5, 448.5, 443.1, 442.7, 441.6, 438.8, 433.4, 432.2, 428.6, 427.4, 423.2, 423, 421.2) # baltic_upper <- c(481.5, 479.7, 475, 472.5, 468.7, 464.7, 463.0, 460.0, 458.4, 457.5, 453.5, 452.5, 451.8, 450.3, 449.5, 448.5, 443.1, 442.7, 441.6, 438.8, 433.4, 432.2, 428.6, 427.4, 423.2, 423, 421.2, 419.2) baltic_mid <- (baltic_lower + baltic_upper)/2 baltic_duration <- baltic_lower - baltic_upper baltic_stage_times <- cbind(baltic_stages, baltic_lower, baltic_upper, baltic_mid, baltic_duration) baltic_stage_times <-as.data.frame(baltic_stage_times) bstages <- cbind(baltic_stages, baltic_lower, baltic_upper) ### gamma baltic stages level bb_a <- brachs_binned_es[,c("genus", "stage1")] colnames(bb_a) <- c("genus", "bb") bb_b <- brachs_binned_es[,c("genus", "stage2")] colnames(bb_b) <- c("genus", "bb") bb_t <- unique(rbind(bb_a, bb_b)) occs.sb <- table(bb_t[,c("genus","bb")]) occs.sb <- as.data.frame.matrix(occs.sb) occs.sb <- subset(occs.sb, rowSums(occs.sb[,])>0 ) bin.name <- array() rx <- occs.sb rxh <- matrix(,nrow(occs.sb), length(baltic_stages)) for (k in 1:(length(baltic_stages))) { cns <- as.character(baltic_stages[k]) bin.name[k] <- cns if (cns %in% colnames(rx) == TRUE) { rxh[,k] <- rx[,cns] } } colnames(rxh) <- bin.name rxh <- rxh[,-(1:3)] rxh[is.na(rxh)] <- 0 rxh.full <- rxh rxh[rxh > 0] <- 1 rxh <- subset(rxh, rowSums(rxh)>0 ) colSums(rxh) colSums(rxh.full) cmrarr <- as.data.frame(matrix(,nrow(rxh), 1)) for(kk in 1:nrow(rxh)) { cmrarr[kk,] <- paste(rxh[kk,1], rxh[kk,2], rxh[kk,3], rxh[kk,4], rxh[kk,5], rxh[kk,6], rxh[kk,7], rxh[kk,8], rxh[kk,9], rxh[kk,10], rxh[kk,11], rxh[kk,12], rxh[kk,13], rxh[kk,14], rxh[kk,15], rxh[kk,16], rxh[kk,17], rxh[kk,18], rxh[kk,19], rxh[kk,20], rxh[kk,21], rxh[kk,22], rxh[kk,23], rxh[kk,24], rxh[kk,25], sep="") # with Sil rxh[kk,26], rxh[kk,27], rxh[kk,28] cmrarr[kk,] <- as.character(cmrarr[kk,]) } colnames(cmrarr) <- "ch" osq.in <- cmrarr # an absence/presence matrix of all genus observations (Sobs) per time bin prepared for RMark # Models Phitime=list(formula=~time) ptime=list(formula=~time) Penttime=list(formula=~time) Gammatime=list(formula=~time) #POPAN all.processed <- process.data(osq.in, model="POPAN", time.intervals = baltic_duration[-(1:4)]) all.ddl <- make.design.data(all.processed) all.time.ch <- mark(all.processed,all.ddl, model.parameters=list(Phi=Phitime, p=ptime, pent=Penttime)) div.osq <- all.time.ch$results$derived$N plot(-baltic_mid[-(1:3)], div.osq$estimate, type="o", ylim=c(0,max(div.osq$ucl))) arrows(-baltic_mid[-(1:3)], (div.osq$lcl), -baltic_mid[-(1:3)], (div.osq$ucl), length=0.05, angle=90, code=3) gamma_resis <- div.osq gamma_resis$gen_count <- colSums(rxh) gamma_resis$mid_age <- -as.numeric(as.character(baltic_stage_times$baltic_mid[-(1:3)])) gamma_resis$duration <- baltic_stage_times$baltic_duration[-(1:3)] gamma_resis$stages <- colnames(rxh) gamma_stage_resis <- gamma_resis ############ ########### gamma time bin resolution bb_a <- brachs_binned_es[,c("genus", "bb1")] colnames(bb_a) <- c("genus", "stage") bb_b <- brachs_binned_es[,c("genus", "bb2")] colnames(bb_b) <- c("genus", "stage") bb_t <- unique(rbind(bb_a, bb_b)) occs.sb <- table(bb_t[,c("genus","stage")]) occs.sb <- as.data.frame.matrix(occs.sb) occs.sb <- subset(occs.sb, rowSums(occs.sb[,])>0 ) bin.name <- array() rx <- occs.sb rxh <- matrix(,nrow(occs.sb), length(tbins)) for (k in 1:(length(tbins))) { cns <- as.character(tbins[k]) bin.name[k] <- cns if (cns %in% colnames(rx) == TRUE) { rxh[,k] <- rx[,cns] } } colnames(rxh) <- bin.name rxh <- rxh[,-(1:2)] rxh[is.na(rxh)] <- 0 rxh.full <- rxh rxh[rxh > 0] <- 1 rxh <- subset(rxh, rowSums(rxh)>0 ) colSums(rxh) colSums(rxh.full) cmrarr <- as.data.frame(matrix(,nrow(rxh), 1)) for(kk in 1:nrow(rxh)) { cmrarr[kk,] <- paste(rxh[kk,1], rxh[kk,2], rxh[kk,3], rxh[kk,4], rxh[kk,5], rxh[kk,6], rxh[kk,7], rxh[kk,8], rxh[kk,9], rxh[kk,10], rxh[kk,11], sep="") # with Sil , rxh[kk,13] , rxh[kk,12] cmrarr[kk,] <- as.character(cmrarr[kk,]) } colnames(cmrarr) <- "ch" osq.in <- cmrarr # an absence/presence matrix of all genus observations (Sobs) per time bin prepared for RMark # Models Phitime=list(formula=~time) ptime=list(formula=~time) Penttime=list(formula=~time) Gammatime=list(formula=~time) #POPAN all.processed <- process.data(osq.in, model="POPAN", time.intervals = tbins_dur[-(1:3)]) all.ddl <- make.design.data(all.processed) all.time.ch <- mark(all.processed,all.ddl, model.parameters=list(Phi=Phitime, p=ptime, pent=Penttime)) div.osq <- all.time.ch$results$derived$N plot(-tbins_mid[-(1:2)], div.osq$estimate,col="blue", type="o", ylim=c(0,max(div.osq$ucl))) arrows(-tbins_mid[-(1:2)], (div.osq$lcl), -tbins_mid[-(1:2)], (div.osq$ucl), col="blue",length=0.05, angle=90, code=3) lines(-baltic_mid[-(1:4)], gamma_stage_resis$estimate, col="red") arrows(-baltic_mid[-(1:4)], (gamma_stage_resis$lcl), -baltic_mid[-(1:4)], (gamma_stage_resis$ucl), col="red",length=0.05, angle=90, code=3) abline(v = -tbins_lower[], col="blue", lwd=1, lty=2) abline(v = -baltic_mid[c(3,4,5,6)], col="red", lwd=1, lty=2) gamma_resis <- div.osq gamma_resis$gen_count <- colSums(rxh) gamma_resis$mid_age <- -tbins_mid[-(1:2)] gamma_resis$bins <- colnames(rxh) gamma_bins_resis <- gamma_resis save(gamma_bins_resis, gamma_stage_resis, baltic_stage_times,file="cr_results_06_02_21corrected.Rdata")