### R script to extract survival data require(stringr) require(plyr) require(ggplot2) require(grid) require(scales) require(gridExtra) ##################### ### DATA PREPARATION ### ##################### ### Load file - recorded from the tab 2007-2013 and save as *.csv data1 <- read.csv("C:/Users/ygager/Desktop/20140527_Molmol_caps.csv") #################################### ### Cleaning of the raw datafile ### #################################### # A lot of code that are not of interest and therefore removed ############################################################## ### Raw data filtered to obtain only variables of interest ### ############################################################## data2f <- data1 data2f$date <- as.Date(data2f$date,"%m/%d/%Y") data2f$year <- year(data2f$date) # A lot of code consisting of renaming data for age of bats, site, was also removed ### Filtering site data keeproost <- c("100","147","150","152A","154A","154B","157","164","165","257","270","271","275","451") # Roosts with several recaptures and possibility to determine group size ### Automatic function to keep dataset with roosts with identified group sizes which(data2f[,"site"] == keeproost[1]) ### List to ID filtered roosts and create a filtered dataframe possiteL <- list() for(i in 1:length(keeproost)){ possiteL[[i]] <- which(data2f[,"site"] == keeproost[i]) #position of transponder in the CMR data } siteL <- unlist(possiteL) data2f <- data2f[siteL,] sort(unique(data2f$site)) ## Plot site and dates selected (after filtering) ggplot(data=data2f) + geom_point(aes(x=date,y=1,color=site),size=3) + ylab("Capture event") + xlab("Date") + geom_line(aes(x=date,y=1,color=site),linetype="dashed") + facet_grid(site ~ .) + theme(legend.position="none") ### Calculation number of times of transponder is recaptured fqtransp <- rle(sort(data2f$transponder)) fqtransp$lengths sort(fqtransp$lengths) freq = freq(fqtransp$lengths) plot(freq,xlab="# captures",ylab="Individuals",main="#captures per individual (20140108)") ############################ ### Determine group size ### ############################ ### Once the transponders are all 5, creation of a vector of unique transponders utransp <- unique(data2f$transponder) utransp <- as.factor(utransp) ### Calculation number of times of transponder is recaptured fqtransp <- rle(sort(data2f$transponder)) fqtransp$lengths ### Exchange between roosts trcheck <- paste(data2f$transponder,data2f$site) disptransp <- length(unique(trcheck)) disptransp - length(utransp) #Number of transponders that changed roosts ### Debriefing of the filtered dataset (adults) length(utransp) #number of unique transponders freq(fqtransp$lengths) #Percent of capture/recapture per transponder sort(unique(data2f$date)) #vector with the nights of capture length(sort(unique(data2f$date))) #number of night of captures ### Column with date & site data2f$match <- paste(data2f$date,data2f$site) ucatch <- unique(data2f$match) #sort(ucatch) sort(ucatch) ### Calculation of group size per capture dates length(which(data2f$match == ucatch[1])) ### Calculation of group size per capture date GSL <- list() for(i in 1: length(ucatch)){ GSL[[i]] <- length(which(data2f$match == ucatch[i])) } names(GSL) <- ucatch GSL ### Infer group size to the filtered data frame posroostL <- list() roostL <- list() data2f$gs <- NA for(i in 1:length(ucatch)){ posroostL[[i]] <- which(data2f$match == ucatch[i]) roostL[[i]] <- unique(data2f[posroostL[[i]],"site"]) data2f[posroostL[[i]],"gs"] <- GSL[[i]] } ### Plot the group variation per site over time (not season but real date) #unique date and site uGSL <- unlist(GSL) dfGSL <- as.data.frame(uGSL) colnames(dfGSL)[1] <- "gs" dfGSL$datesite <- row.names(dfGSL) dfGSL$site <- str_sub(dfGSL$datesite, 12) dfGSL$date <- str_sub(dfGSL$datesite, 1,10) dfGSL$date <- as.Date(dfGSL$date) dfGSL$month <- month(dfGSL$date) ### Mean of median group size over all captures and all sites gs_meds <- ddply(dfGSL, .(site), summarise, med = median(gs)) meanmed <- mean(gs_meds$med) ### Histogram with normal curve x <- data.frame(aa <- dfGSL$gs) b <- ggplot(x, aes(aa)) + geom_histogram(aes(y=..density..),binwidth = 1) + stat_function(fun=dnorm, lwd=1,args=list(mean=mean(x$aa), sd=sd(x$aa))) + xlab("Group size") + ylab("Count") b ### Density curve d <- density(dfGSL$gs) # returns the density data plot(d,xlab="Group size", main =NULL) ### Plot GS ~ date (grid site + points + line) ggplot(data=dfGSL) + geom_point(aes(x=date,y=gs,color=site),size=3) + ylab("Group size") + xlab("Date") + geom_line(aes(x=date,y=gs,color=site),linetype="dashed") + facet_grid(site ~ .) + theme(legend.position="none") ### dfsite <- as.data.frame(unique(dfGSL$site)) vlength <- 1:nrow(uniksite) dfsite$site2 <- vlength colnames(dfsite) <- c("site","site2") msites <- match(dfGSL$site,dfsite$site) dfGSL$Site <- as.factor(dfsite[msites,"site2"]) ### Plot GS ~ month (points + boxplot + smooth loess - all sites combined) = good plot for survival paper tabmon <- as.data.frame(unlist(table(dfGSL$month))) month <- c(1:12) mfreq <- tabmon[match(month,tabmon$Var1),"Freq"] tabmon <- cbind(month,mfreq) tabmon[is.na(tabmon)] <- 0 tabmon <- as.data.frame(tabmon) ### Max gs per month monthL <- dlply(dfGSL,.(month)) maxL <- list() for(i in 1:length(monthL)){ maxL[[i]] <- max(monthL[[i]]$gs) } vnojune <- unlist(maxL) tabmon$max <- c(4,13,13,18,15,0,26,22,16,10,14,11) ggplot(data=dfGSL, aes(x=factor(month,levels=1:12),gs)) + scale_x_discrete(breaks = seq(min(dfGSL$month), max(dfGSL$month), by = 1),drop=FALSE, labels=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec")) + stat_boxplot(geom ="errorbar",width=0.5) + geom_boxplot(fill="lightgray") + xlab("Month") + ylab("Total group size") + theme_classic() + geom_text(data=tabmon,aes(x=month,y=max+2,label=mfreq),size=3) + theme(axis.ticks.length=unit(-0.12, "cm"), axis.ticks.margin=unit(0.3, "cm"),text = element_text(size=12),axis.title.y = element_text(vjust=1.2), axis.title.x = element_text(vjust=-0.2)) #+ annotate("text",x=3,y=25,label="81 capture events") #+ geom_point(aes(x=factor(month),y=gs,color=site),size=3) table(dfGSL$site) ### Proportion of adults and juveniles # 656 adults marked, including 490 length(unique(data2f$transponder)) table(data2f$age) #http://stackoverflow.com/questions/12993545/how-to-put-whisker-ends-on-ggplot2-boxplot ### Save figure for publication #ggsave(file="Totalgs_Molmol_Gamboa.pdf") ggsave(file="Totalgs_Molmol_Gamboa.pdf",width=14,height=8,unit="cm") ### Save table for data accessibility dataTGS <- dfGSL[,c("site","date","month","gs")] colnames(dataTGS) <- c("group","date","month","totalgs") #write.csv(dataTGS,file="Totalgs_Molossus_Gamboa.csv",row.names=FALSE) #ggplot(data=dfGSL, aes(x=factor(month,levels=1:12),gs)) + scale_x_discrete(breaks = seq(min(dfGSL$month), max(dfGSL$month), by = 1),drop=FALSE) + geom_boxplot(fill="lightgray") + xlab("Month") + ylab("Total group size") + theme_bw() + annotate("text",x=3,y=25,label="81 capture events") ### Measurements of median round(median(dfGSL$gs),1) ### Measurements of mean round(mean(dfGSL$gs),1) round(sd(dfGSL$gs),1) ### Plot GS ~ month (facet site + points + boxplot) ggplot(data=dfGSL, aes(as.factor(month),gs)) + geom_point(aes(x=as.factor(month),y=gs,color=site),size=3) + geom_boxplot(aes(as.factor(month),gs)) + xlab("Month") + ylab("Group size") #+ facet_wrap(~ site) ### Plot GS ~ site (points + boxplot) - Order by median ggplot(data=dfGSL, aes(x = reorder(site, gs, FUN=median),y=gs)) + geom_boxplot(fill="lightgrey") + ylab("Group size") + xlab("Roost #") + geom_point(aes(x = reorder(site, gs, FUN=median),y=gs,color=site),size=3) + theme(legend.position="none") + theme_bw() + theme(legend.position=none)# Treshold mean of medians ! + geom_hline(aes(yintercept=meanmed),colour="gray") ### Information on group size data round(mean(dfGSL$gs),1) round(sd(dfGSL$gs),1) ### Attribute group size to datadf mdates <- match(data2f$date,dfGSL$date) data2f$gs <- dfGSL[mdates,"gs"] ### Check of the single individuals which(dfGSL$gs == 1) data2f[data2f$date == "2014-05-13",] #152A f data2f[data2f$date == "2010-03-25",] #164 m data2f[data2f$date == "2011-12-16",] #164 m