Permalink
Cannot retrieve contributors at this time
Name already in use
A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
Francomano_et_al_2023_PAM_for_Penguins/Statistical_Analysis_and_Figures_Script.R
Go to fileThis commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
892 lines (793 sloc)
70.8 KB
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ###################################################################### | |
| # Code from Francomano et al. (IBIS 2023) # | |
| ###################################################################### | |
| # Statistical analysis and figures script # | |
| ###################################################################### | |
| # Created by: Dante Francomano (dfrancomano@alumni.purdue.edu) | |
| # Date: 20221113 | |
| ###################################################################### | |
| ###################################################################### | |
| # Load packages | |
| library(vegan) | |
| library(BiodiversityR) | |
| library(lme4) | |
| library(r2glmm) | |
| library(parallel) | |
| library(lmPerm) | |
| library(dplyr) | |
| library(suncalc) | |
| library(RColorBrewer) | |
| library(effects) | |
| library(pBrackets) | |
| # Read in data | |
| soundCountsRockhopperRaw<-read.csv("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Conteos_de_Sonidos/Datos/Datos_Completos/CSV/Forma_de_Datos_de_Cuentas_de_Sonidos_Estados_PPA_1.csv",na.strings="N/A") | |
| soundCountsMagellanicRaw<-lapply(list.files("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Conteos_de_Sonidos/Datos/Datos_Completos/CSV",full.names=TRUE)[3:6],read.csv,na.strings="N/A") | |
| soundCountsMagellanicRaw<-rbind(soundCountsMagellanicRaw[[1]],soundCountsMagellanicRaw[[2]],soundCountsMagellanicRaw[[3]],soundCountsMagellanicRaw[[4]]) | |
| acousticMetricsMagellanic<-read.csv("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Acoustic_Metrics/Acoustic_Metrics_PM.csv") | |
| acousticMetricsRockhopper<-read.csv("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Acoustic_Metrics/Acoustic_Metrics_PPA.csv") | |
| cameraTrapCountsMagellanic<-read.csv("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Datos_de_Cámaras_Trampa/PM_Data.csv") | |
| cameraTrapCountsRockhopper<-read.csv("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Datos_de_Cámaras_Trampa/PPA_Data.csv") | |
| populationDensitiesMagellanic<-read.csv("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Conteos_de_Pingüinos/PM_Densities_by_Radii.csv") | |
| populationDensitiesRockhopper<-read.csv("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Conteos_de_Pingüinos/PPA_Densities_by_Radii.csv") | |
| recorderLocations<-read.csv("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Site_Coordinates.csv") | |
| # Clean data | |
| cameraTrapCountsMagellanic$Site[which(cameraTrapCountsMagellanic$Site=="Zona_Tur\x92stica")]<-"Zona_Turistica" | |
| acousticMetricsMagellanic$Date<-as.ordered(acousticMetricsMagellanic$Date) | |
| acousticMetricsRockhopper$Date<-ordered(acousticMetricsRockhopper$Date,levels=levels(acousticMetricsMagellanic$Date)) | |
| acousticMetricsMagellanic$Time<-as.factor(acousticMetricsMagellanic$Time) | |
| acousticMetricsRockhopper$Time<-as.factor(acousticMetricsRockhopper$Time) | |
| cameraTrapCountsMagellanic$Site<-as.factor(cameraTrapCountsMagellanic$Site) | |
| cameraTrapCountsRockhopper$Site<-as.factor(cameraTrapCountsRockhopper$Site) | |
| cameraTrapCountsMagellanic$Breeding_Stage<-as.factor(cameraTrapCountsMagellanic$Breeding_Stage) | |
| cameraTrapCountsMagellanic$Date<-ordered(cameraTrapCountsMagellanic$Date,levels=levels(acousticMetricsMagellanic$Date)) | |
| cameraTrapCountsRockhopper$Date<-ordered(cameraTrapCountsRockhopper$Date,levels=levels(acousticMetricsMagellanic$Date)) | |
| cameraTrapCountsMagellanic<-cameraTrapCountsMagellanic[!is.na(cameraTrapCountsMagellanic$Count_2),] | |
| cameraTrapCountsRockhopper<-cameraTrapCountsRockhopper[!is.na(cameraTrapCountsRockhopper$Count_1),] | |
| populationDensitiesMagellanic$Site<-as.factor(c("PM_01","PM_02","19","Baliza","P53","Rata_Muerta","Zona_Erosionada","Zona_Nueva","Zona_Turistica")) | |
| populationDensitiesRockhopper$Site<-unique(cameraTrapCountsRockhopper$Site) | |
| ###################################################################### | |
| # Plot sound occurrence by breeding stage and species | |
| ###################################################################### | |
| # Remove unused points and observations with lost contact | |
| soundCountsRockhopper<-soundCountsRockhopperRaw[which(!is.na(soundCountsRockhopperRaw$F)&soundCountsRockhopperRaw$CP!=1&soundCountsRockhopperRaw$SPV!=1),] | |
| soundCountsMagellanic<-soundCountsMagellanicRaw[which(!is.na(soundCountsMagellanicRaw$F)&soundCountsMagellanicRaw$CP!=1&soundCountsMagellanicRaw$SPV!=1),] | |
| # Determine breeding stage | |
| soundCountsMagellanic$Breeding_Stage<-NA | |
| for(i in 1:dim(soundCountsMagellanic)[1]) | |
| { | |
| if(soundCountsMagellanic$F[i]<20181001) soundCountsMagellanic$Breeding_Stage[i]<-"arrival" | |
| if(soundCountsMagellanic$F[i]>20180930&soundCountsMagellanic$F[i]<20181021) soundCountsMagellanic$Breeding_Stage[i]<-"pair_formation" | |
| if(soundCountsMagellanic$F[i]>20181020&soundCountsMagellanic$F[i]<20181121) soundCountsMagellanic$Breeding_Stage[i]<-"incubation" | |
| if(soundCountsMagellanic$F[i]>20181120&soundCountsMagellanic$F[i]<20181201) soundCountsMagellanic$Breeding_Stage[i]<-"hatching" | |
| if(soundCountsMagellanic$F[i]>20181130&soundCountsMagellanic$F[i]<20190101) soundCountsMagellanic$Breeding_Stage[i]<-"early_care" | |
| if(soundCountsMagellanic$F[i]>20181231&soundCountsMagellanic$F[i]<20190215) soundCountsMagellanic$Breeding_Stage[i]<-"late_care" | |
| if(soundCountsMagellanic$F[i]>20190214&soundCountsMagellanic$F[i]<20190311) soundCountsMagellanic$Breeding_Stage[i]<-"pre-molt_feeding" | |
| if(soundCountsMagellanic$F[i]>20190310) soundCountsMagellanic$Breeding_Stage[i]<-"molting" | |
| } | |
| soundCountsRockhopper$Breeding_Stage<-"early_care" | |
| # Determine numbers of observations and numbers of vocalizations | |
| soundCountResults<-data.frame(Species=c("Magellanic","rockhopper","Magellanic","Magellanic","Magellanic"), | |
| Breeding_Stage=c("hatching/early_care","hatching/early_care","late care","pre-molt_feeding","molting"), | |
| Number_of_Observations=NA,Number_of_Observations_with_Vocalization=NA,Number_of_Observations_with_More_than_1_Vocalization=NA) | |
| soundCountResults$Number_of_Observations[1]<-dim(soundCountsMagellanic[which(soundCountsMagellanic$Breeding_Stage=="hatching"|soundCountsMagellanic$Breeding_Stage=="early_care"),])[1] | |
| soundCountResults$Number_of_Observations[2]<-dim(soundCountsRockhopper[which(soundCountsMagellanic$Breeding_Stage=="early_care"),])[1] | |
| soundCountResults$Number_of_Observations[3]<-dim(soundCountsMagellanic[which(soundCountsMagellanic$Breeding_Stage=="late_care"),])[1] | |
| soundCountResults$Number_of_Observations[4]<-dim(soundCountsMagellanic[which(soundCountsMagellanic$Breeding_Stage=="pre-molt_feeding"),])[1] | |
| soundCountResults$Number_of_Observations[5]<-dim(soundCountsMagellanic[which(soundCountsMagellanic$Breeding_Stage=="molting"),])[1] | |
| soundCountResults$Number_of_Observations_with_Vocalization[1]<-dim(soundCountsMagellanic[which((soundCountsMagellanic$Breeding_Stage=="hatching"|soundCountsMagellanic$Breeding_Stage=="early_care")&(soundCountsMagellanic$LDE>0|soundCountsMagellanic$LDM>0|soundCountsMagellanic$SHI>0)),])[1] | |
| soundCountResults$Number_of_Observations_with_Vocalization[2]<-dim(soundCountsRockhopper[which(soundCountsRockhopper$Breeding_Stage=="early_care"&(soundCountsRockhopper$LDE>0|soundCountsRockhopper$LDM>0|soundCountsRockhopper$GRU>0|soundCountsRockhopper$GRI>0)),])[1] | |
| soundCountResults$Number_of_Observations_with_Vocalization[3]<-dim(soundCountsMagellanic[which(soundCountsMagellanic$Breeding_Stage=="late_care"&(soundCountsMagellanic$LDE>0|soundCountsMagellanic$LDM>0|soundCountsMagellanic$SHI>0)),])[1] | |
| soundCountResults$Number_of_Observations_with_Vocalization[4]<-dim(soundCountsMagellanic[which(soundCountsMagellanic$Breeding_Stage=="pre-molt_feeding"&(soundCountsMagellanic$LDE>0|soundCountsMagellanic$LDM>0|soundCountsMagellanic$SHI>0)),])[1] | |
| soundCountResults$Number_of_Observations_with_Vocalization[5]<-dim(soundCountsMagellanic[which(soundCountsMagellanic$Breeding_Stage=="molting"&(soundCountsMagellanic$LDE>0|soundCountsMagellanic$LDM>0|soundCountsMagellanic$SHI>0)),])[1] | |
| soundCountResults$Number_of_Observations_with_More_than_1_Vocalization[1]<-dim(soundCountsMagellanic[which((soundCountsMagellanic$Breeding_Stage=="hatching"|soundCountsMagellanic$Breeding_Stage=="early_care")&(soundCountsMagellanic$LDE>1|soundCountsMagellanic$LDM>1|soundCountsMagellanic$SHI>1)),])[1] | |
| soundCountResults$Number_of_Observations_with_More_than_1_Vocalization[2]<-dim(soundCountsRockhopper[which(soundCountsRockhopper$Breeding_Stage=="early_care"&(soundCountsRockhopper$LDE>1|soundCountsRockhopper$LDM>1|soundCountsRockhopper$GRU>1|soundCountsRockhopper$GRI>1)),])[1] | |
| soundCountResults$Number_of_Observations_with_More_than_1_Vocalization[3]<-dim(soundCountsMagellanic[which(soundCountsMagellanic$Breeding_Stage=="late_care"&(soundCountsMagellanic$LDE>1|soundCountsMagellanic$LDM>1|soundCountsMagellanic$SHI>1)),])[1] | |
| soundCountResults$Number_of_Observations_with_More_than_1_Vocalization[4]<-dim(soundCountsMagellanic[which(soundCountsMagellanic$Breeding_Stage=="pre-molt_feeding"&(soundCountsMagellanic$LDE>1|soundCountsMagellanic$LDM>1|soundCountsMagellanic$SHI>1)),])[1] | |
| soundCountResults$Number_of_Observations_with_More_than_1_Vocalization[5]<-dim(soundCountsMagellanic[which(soundCountsMagellanic$Breeding_Stage=="molting"&(soundCountsMagellanic$LDE>1|soundCountsMagellanic$LDM>1|soundCountsMagellanic$SHI>1)),])[1] | |
| # Plot results | |
| pdf("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Figures/Sound_Count_Results.pdf",width=5,height=4,pointsize=6) | |
| par(mar=c(10,6,2,2)) | |
| plot(1,1,xlim=c(0.5,5.5),ylim=c(0,85),type="n",axes=FALSE,xlab="",ylab="") | |
| for(i in 1:5) polygon(x=c((i-0.4),(i+0.4),(i+0.4),(i-0.4)),y=c(0,0,rep(soundCountResults$Number_of_Observations_with_Vocalization[i],2)),border="grey60",col="grey60",lwd=1) | |
| for(i in 1:5) polygon(x=c((i-0.4),(i+0.4),(i+0.4),(i-0.4)),y=c(0,0,rep(soundCountResults$Number_of_Observations_with_More_than_1_Vocalization[i],2)),border="grey30",col="grey30",lwd=1) | |
| for(i in 1:5) polygon(x=c((i-0.4),(i+0.4),(i+0.4),(i-0.4)),y=c(0,0,rep(soundCountResults$Number_of_Observations[i],2)),border="grey30",lwd=1) | |
| for(i in 1:5) text(i,soundCountResults$Number_of_Observations_with_Vocalization[i]+3,paste(round(100*soundCountResults$Number_of_Observations_with_Vocalization[i]/soundCountResults$Number_of_Observations[i]),"%",sep=""),col="grey30",cex=1.5) | |
| axis(side=1,at=1:5,labels=c("Magellanic\nhatching/\nearly chick\nrearing","Rockhopper\nhatching/\nearly chick\nrearing","Magellanic\nlate chick\nrearing","Magellanic\npre-molt\nfeeding","Magellanic\nmolting"),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=0,line=-2,padj=1,las=1) | |
| axis(side=2,at=seq(0,80,10),labels=seq(0,80,10),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-0.75,padj=0.5,las=1) | |
| polygon(x=c(0.4,5.6,5.6,0.4),y=c(0,0,85,85),border="grey30",col=NA) | |
| mtext("species and breeding stage",side=1,line=7,col="grey30",cex=2) | |
| mtext("number of observations",side=2,line=3,col="grey30",cex=2) | |
| legend(5.49,82.6,xjust=1,yjust=1,c("no vocalizations observed","one vocalization observed","multiple vocalizations observed"),fill=c("white","grey60","grey30"),border="grey30",box.col="grey30",cex=1.5) | |
| dev.off() | |
| ###################################################################### | |
| # Examine dimensionality of acoustic data to determine acoustic variables to use | |
| ###################################################################### | |
| # Visually check distributions of values for approvimate normality | |
| apply(acousticMetricsMagellanic[,9:14],2,hist) | |
| apply(acousticMetricsRockhopper[,c(9,11,12,15)],2,hist) | |
| # Examine inter-metric correlations (using Spearman correlations for temporal occupancy and events per minute, given their distributions) | |
| interMetricCorrelationResultsMagellanic<-matrix(NA,6,6) | |
| rownames(interMetricCorrelationResultsMagellanic)<-colnames(interMetricCorrelationResultsMagellanic)<-colnames(acousticMetricsMagellanic[,9:14]) | |
| interMetricCorrelationResultsRockhopper<-matrix(NA,4,4) | |
| rownames(interMetricCorrelationResultsRockhopper)<-colnames(interMetricCorrelationResultsRockhopper)<-colnames(acousticMetricsRockhopper[,c(9,11,12,15)]) | |
| interMetricCorrelationResultsMagellanic[2,1]<-round(cor.test(acousticMetricsMagellanic$SPL,acousticMetricsMagellanic$BN_20,method="pearson")$estimate,2) | |
| interMetricCorrelationResultsMagellanic[3,1]<-round(cor.test(acousticMetricsMagellanic$SPL,acousticMetricsMagellanic$ACI,method="pearson")$estimate,2) | |
| interMetricCorrelationResultsMagellanic[4,1]<-round(cor.test(acousticMetricsMagellanic$SPL,acousticMetricsMagellanic$H,method="pearson")$estimate,2) | |
| interMetricCorrelationResultsMagellanic[5,1]<-round(cor.test(acousticMetricsMagellanic$SPL,acousticMetricsMagellanic$TO,method="spearman")$estimate,2) | |
| interMetricCorrelationResultsMagellanic[6,1]<-round(cor.test(acousticMetricsMagellanic$SPL,acousticMetricsMagellanic$EPM,method="spearman")$estimate,2) | |
| interMetricCorrelationResultsMagellanic[3,2]<-round(cor.test(acousticMetricsMagellanic$BN_20,acousticMetricsMagellanic$ACI,method="pearson")$estimate,2) | |
| interMetricCorrelationResultsMagellanic[4,2]<-round(cor.test(acousticMetricsMagellanic$BN_20,acousticMetricsMagellanic$H,method="pearson")$estimate,2) | |
| interMetricCorrelationResultsMagellanic[5,2]<-round(cor.test(acousticMetricsMagellanic$BN_20,acousticMetricsMagellanic$TO,method="spearman")$estimate,2) | |
| interMetricCorrelationResultsMagellanic[6,2]<-round(cor.test(acousticMetricsMagellanic$BN_20,acousticMetricsMagellanic$EPM,method="spearman")$estimate,2) | |
| interMetricCorrelationResultsMagellanic[4,3]<-round(cor.test(acousticMetricsMagellanic$ACI,acousticMetricsMagellanic$H,method="pearson")$estimate,2) | |
| interMetricCorrelationResultsMagellanic[5,3]<-round(cor.test(acousticMetricsMagellanic$ACI,acousticMetricsMagellanic$TO,method="spearman")$estimate,2) | |
| interMetricCorrelationResultsMagellanic[6,3]<-round(cor.test(acousticMetricsMagellanic$ACI,acousticMetricsMagellanic$EPM,method="spearman")$estimate,2) | |
| interMetricCorrelationResultsMagellanic[5,4]<-round(cor.test(acousticMetricsMagellanic$H,acousticMetricsMagellanic$TO,method="spearman")$estimate,2) | |
| interMetricCorrelationResultsMagellanic[6,4]<-round(cor.test(acousticMetricsMagellanic$H,acousticMetricsMagellanic$EPM,method="spearman")$estimate,2) | |
| interMetricCorrelationResultsMagellanic[6,5]<-round(cor.test(acousticMetricsMagellanic$TO,acousticMetricsMagellanic$EPM,method="spearman")$estimate,2) | |
| interMetricCorrelationResultsRockhopper[2,1]<-round(cor.test(acousticMetricsRockhopper$SPL,acousticMetricsRockhopper$ACI,method="pearson")$estimate,2) | |
| interMetricCorrelationResultsRockhopper[3,1]<-round(cor.test(acousticMetricsRockhopper$SPL,acousticMetricsRockhopper$H,method="pearson")$estimate,2) | |
| interMetricCorrelationResultsRockhopper[4,1]<-round(cor.test(acousticMetricsRockhopper$SPL,acousticMetricsRockhopper$BN_10,method="pearson")$estimate,2) | |
| interMetricCorrelationResultsRockhopper[3,2]<-round(cor.test(acousticMetricsRockhopper$ACI,acousticMetricsRockhopper$H,method="pearson")$estimate,2) | |
| interMetricCorrelationResultsRockhopper[4,2]<-round(cor.test(acousticMetricsRockhopper$ACI,acousticMetricsRockhopper$BN_10,method="pearson")$estimate,2) | |
| interMetricCorrelationResultsRockhopper[4,3]<-round(cor.test(acousticMetricsRockhopper$H,acousticMetricsRockhopper$BN_10,method="pearson")$estimate,2) | |
| interMetricCorrelationResultsMagellanic | |
| interMetricCorrelationResultsRockhopper | |
| # Conduct PCAs and examine axis loadings and significance of axes (dropping ACI, as it was just included for exploratory purposes) | |
| PCAMagellanic<-rda(acousticMetricsMagellanic[!is.na(acousticMetricsMagellanic$SPL),c(9,10,12:14)],scale=TRUE) | |
| PCARockhopper<-rda(acousticMetricsRockhopper[!is.na(acousticMetricsRockhopper$SPL),c(9,12,15)],scale=TRUE) | |
| summary(PCAMagellanic)$species | |
| summary(PCARockhopper)$species | |
| plot(1:5,summary(PCAMagellanic)$cont[[1]][2,],type="b",ylim=c(0,1)) | |
| plot(1:3,summary(PCARockhopper)$cont[[1]][2,],type="b",ylim=c(0,1)) | |
| PCAsignificance(PCAMagellanic) | |
| PCAsignificance(PCARockhopper) | |
| # Conduct and examine PCAs on a per-site basis to further assess dimensionality | |
| sitesMagellanic<-unique(acousticMetricsMagellanic$Site) | |
| sitesRockhopper<-unique(acousticMetricsRockhopper$Site) | |
| sitePCAsMagellanic<-lapply(sitesMagellanic,function(x) rda(acousticMetricsMagellanic[which(!is.na(acousticMetricsMagellanic$SPL)&acousticMetricsMagellanic$Site==x),c(9,10,12:14)],scale=TRUE)) | |
| sitePCAsRockhopper<-lapply(sitesRockhopper,function(x) rda(acousticMetricsRockhopper[which(!is.na(acousticMetricsRockhopper$SPL)&acousticMetricsRockhopper$Site==x),c(9,12,15)],scale=TRUE)) | |
| lapply(sitePCAsMagellanic,function(x) summary(x)$species) | |
| lapply(sitePCAsRockhopper,function(x) summary(x)$species) | |
| lapply(sitePCAsMagellanic,function(x) plot(1:5,summary(x)$cont[[1]][2,],type="b",ylim=c(0,1))) | |
| lapply(sitePCAsRockhopper,function(x) plot(1:3,summary(x)$cont[[1]][2,],type="b",ylim=c(0,1))) | |
| lapply(sitePCAsMagellanic,PCAsignificance) | |
| lapply(sitePCAsRockhopper,PCAsignificance) | |
| # Extract first axis of each PCA, flipping the sign for Magellanics, since loadings were negative | |
| acousticMetricsMagellanic$PC1<-NA | |
| acousticMetricsMagellanic[!is.na(acousticMetricsMagellanic$SPL),]$PC1<--summary(PCAMagellanic)$sites[,1] | |
| acousticMetricsRockhopper$PC1<-NA | |
| acousticMetricsRockhopper[!is.na(acousticMetricsRockhopper$SPL),]$PC1<-summary(PCARockhopper)$sites[,1] | |
| ###################################################################### | |
| # Conduct analysis based on camera trap data | |
| ###################################################################### | |
| # Obtain subsets of relevant data, including acoustic data buffers of 1, 3, and 5 files centered on the photo time (2-, 12-, and 22-minute periods) | |
| oneFilePhotoBufferMagellanic<-cameraTrapCountsMagellanic | |
| threeFilePhotoBufferMagellanic<-cameraTrapCountsMagellanic | |
| fiveFilePhotoBufferMagellanic<-cameraTrapCountsMagellanic | |
| for(i in 1:dim(cameraTrapCountsMagellanic)[1]) | |
| { | |
| acousticMetricIndex<-which(acousticMetricsMagellanic$Site==cameraTrapCountsMagellanic$Site[i]&acousticMetricsMagellanic$Date==cameraTrapCountsMagellanic$Date[i]&acousticMetricsMagellanic$Time==cameraTrapCountsMagellanic$Time[i]) | |
| oneFilePhotoBufferMagellanic$PC1[i]<-acousticMetricsMagellanic[acousticMetricIndex,]$PC1 | |
| threeFilePhotoBufferMagellanic$PC1[i]<-mean(c(acousticMetricsMagellanic[acousticMetricIndex,]$PC1,acousticMetricsMagellanic[acousticMetricIndex-1,]$PC1,acousticMetricsMagellanic[acousticMetricIndex+1,]$PC1)) | |
| fiveFilePhotoBufferMagellanic$PC1[i]<-mean(c(acousticMetricsMagellanic[acousticMetricIndex,]$PC1,acousticMetricsMagellanic[acousticMetricIndex-2,]$PC1,acousticMetricsMagellanic[acousticMetricIndex-1,]$PC1,acousticMetricsMagellanic[acousticMetricIndex+1,]$PC1,acousticMetricsMagellanic[acousticMetricIndex+2,]$PC1)) | |
| } | |
| oneFilePhotoBufferRockhopper<-cameraTrapCountsRockhopper | |
| threeFilePhotoBufferRockhopper<-cameraTrapCountsRockhopper | |
| fiveFilePhotoBufferRockhopper<-cameraTrapCountsRockhopper | |
| for(i in 1:dim(cameraTrapCountsRockhopper)[1]) | |
| { | |
| acousticMetricIndex<-which(acousticMetricsRockhopper$Site==cameraTrapCountsRockhopper$Site[i]&acousticMetricsRockhopper$Date==cameraTrapCountsRockhopper$Date[i]&acousticMetricsRockhopper$Time==cameraTrapCountsRockhopper$Time[i]) | |
| oneFilePhotoBufferRockhopper$PC1[i]<-acousticMetricsRockhopper[acousticMetricIndex,]$PC1 | |
| threeFilePhotoBufferRockhopper$PC1[i]<-mean(c(acousticMetricsRockhopper[acousticMetricIndex,]$PC1,acousticMetricsRockhopper[acousticMetricIndex-1,]$PC1,acousticMetricsRockhopper[acousticMetricIndex+1,]$PC1)) | |
| fiveFilePhotoBufferRockhopper$PC1[i]<-mean(c(acousticMetricsRockhopper[acousticMetricIndex,]$PC1,acousticMetricsRockhopper[acousticMetricIndex-2,]$PC1,acousticMetricsRockhopper[acousticMetricIndex-1,]$PC1,acousticMetricsRockhopper[acousticMetricIndex+1,]$PC1,acousticMetricsRockhopper[acousticMetricIndex+2,]$PC1)) | |
| } | |
| # Construct and test global mixed effects models for single files, visually checking normality of error (with Shapiro test as well), homogeneity of variance, and linearity | |
| # Treat site as a fixed effect due to the low number of levels | |
| oneFilePhotoModel1Magellanic2<-lmer(PC1~Count_2+Site+(1|Date),oneFilePhotoBufferMagellanic) | |
| plot(oneFilePhotoModel1Magellanic2) | |
| plot(PC1~Count_2,data=oneFilePhotoBufferMagellanic) | |
| shapiro.test(resid(oneFilePhotoModel1Magellanic2)) | |
| oneFilePhotoModelNullMagellanic2<-lmer(PC1~Site+(1|Date),oneFilePhotoBufferMagellanic) | |
| anova(oneFilePhotoModelNullMagellanic2,oneFilePhotoModel1Magellanic2) | |
| r2beta(oneFilePhotoModel1Magellanic2) | |
| oneFilePhotoModel1Magellanic4<-lmer(PC1~Count_4+Site+(1|Date),oneFilePhotoBufferMagellanic[which(!is.na(oneFilePhotoBufferMagellanic$Count_4)),]) | |
| plot(oneFilePhotoModel1Magellanic4) | |
| plot(PC1~Count_4,data=oneFilePhotoBufferMagellanic[which(!is.na(oneFilePhotoBufferMagellanic$Count_4)),]) | |
| shapiro.test(resid(oneFilePhotoModel1Magellanic4)) | |
| oneFilePhotoModelNullMagellanic4<-lmer(PC1~Site+(1|Date),oneFilePhotoBufferMagellanic[which(!is.na(oneFilePhotoBufferMagellanic$Count_4)),]) | |
| anova(oneFilePhotoModelNullMagellanic4,oneFilePhotoModel1Magellanic4) | |
| r2beta(oneFilePhotoModel1Magellanic4) | |
| oneFilePhotoModel1Rockhopper<-lmer(PC1~Count_1+Site+(1|Date),oneFilePhotoBufferRockhopper) | |
| oneFilePhotoModel2Rockhopper<-lm(PC1~Count_1+Site+Date,oneFilePhotoBufferRockhopper) | |
| plot(oneFilePhotoModel2Rockhopper) | |
| # | |
| # | |
| # | |
| # | |
| shapiro.test(resid(oneFilePhotoModel2Rockhopper)) | |
| oneFilePhotoModelNullRockhopper<-lm(PC1~Site+Date,oneFilePhotoBufferRockhopper) | |
| anova(oneFilePhotoModelNullRockhopper,oneFilePhotoModel2Rockhopper) | |
| r2beta(oneFilePhotoModel2Rockhopper) | |
| # Construct and test global mixed effects models for three files, visually checking normality of error (with Shapiro test as well), homogeneity of variance, and linearity | |
| # Treat site as a fixed effect due to the low number of levels | |
| threeFilePhotoModel1Magellanic2<-lmer(PC1~Count_2+Site+(1|Date),threeFilePhotoBufferMagellanic) | |
| plot(threeFilePhotoModel1Magellanic2) | |
| plot(PC1~Count_2,data=threeFilePhotoBufferMagellanic) | |
| shapiro.test(resid(threeFilePhotoModel1Magellanic2)) | |
| threeFilePhotoModelNullMagellanic2<-lmer(PC1~Site+(1|Date),threeFilePhotoBufferMagellanic) | |
| anova(threeFilePhotoModelNullMagellanic2,threeFilePhotoModel1Magellanic2) | |
| r2beta(threeFilePhotoModel1Magellanic2) | |
| threeFilePhotoModel1Magellanic4<-lmer(PC1~Count_4+Site+(1|Date),threeFilePhotoBufferMagellanic[which(!is.na(threeFilePhotoBufferMagellanic$Count_4)),]) | |
| plot(threeFilePhotoModel1Magellanic4) | |
| plot(PC1~Count_4,data=threeFilePhotoBufferMagellanic[which(!is.na(threeFilePhotoBufferMagellanic$Count_4)),]) | |
| shapiro.test(resid(threeFilePhotoModel1Magellanic4)) | |
| threeFilePhotoModelNullMagellanic4<-lmer(PC1~Site+(1|Date),threeFilePhotoBufferMagellanic[which(!is.na(threeFilePhotoBufferMagellanic$Count_4)),]) | |
| anova(threeFilePhotoModelNullMagellanic4,threeFilePhotoModel1Magellanic4) | |
| r2beta(threeFilePhotoModel1Magellanic4) | |
| threeFilePhotoModel1Rockhopper<-lmer(PC1~Count_1+Site+(1|Date),threeFilePhotoBufferRockhopper) | |
| plot(threeFilePhotoModel1Rockhopper) | |
| plot(PC1~Count_1,data=threeFilePhotoBufferRockhopper) | |
| shapiro.test(resid(threeFilePhotoModel1Rockhopper)) | |
| threeFilePhotoModelNullRockhopper<-lmer(PC1~Site+(1|Date),threeFilePhotoBufferRockhopper) | |
| anova(threeFilePhotoModelNullRockhopper,threeFilePhotoModel1Rockhopper) | |
| r2beta(threeFilePhotoModel1Rockhopper) | |
| # Construct and test global mixed effects models for five files, visually checking normality of error (with Shapiro test as well), homogeneity of variance, and linearity | |
| # Treat site as a fixed effect due to the low number of levels | |
| fiveFilePhotoModel1Magellanic2<-lmer(PC1~Count_2+Site+(1|Date),fiveFilePhotoBufferMagellanic) | |
| plot(fiveFilePhotoModel1Magellanic2) | |
| plot(PC1~Count_2,data=fiveFilePhotoBufferMagellanic) | |
| shapiro.test(resid(fiveFilePhotoModel1Magellanic2)) | |
| fiveFilePhotoModelNullMagellanic2<-lmer(PC1~Site+(1|Date),fiveFilePhotoBufferMagellanic) | |
| anova(fiveFilePhotoModelNullMagellanic2,fiveFilePhotoModel1Magellanic2) | |
| r2beta(fiveFilePhotoModel1Magellanic2) | |
| fiveFilePhotoModel1Magellanic4<-lmer(PC1~Count_4+Site+(1|Date),fiveFilePhotoBufferMagellanic[which(!is.na(fiveFilePhotoBufferMagellanic$Count_4)),]) | |
| plot(fiveFilePhotoModel1Magellanic4) | |
| plot(PC1~Count_4,data=fiveFilePhotoBufferMagellanic[which(!is.na(fiveFilePhotoBufferMagellanic$Count_4)),]) | |
| shapiro.test(resid(fiveFilePhotoModel1Magellanic4)) | |
| fiveFilePhotoModelNullMagellanic4<-lmer(PC1~Site+(1|Date),fiveFilePhotoBufferMagellanic[which(!is.na(fiveFilePhotoBufferMagellanic$Count_4)),]) | |
| anova(fiveFilePhotoModelNullMagellanic4,fiveFilePhotoModel1Magellanic4) | |
| r2beta(fiveFilePhotoModel1Magellanic4) | |
| fiveFilePhotoModel1Rockhopper<-lmer(PC1~Count_1+Site+(1|Date),fiveFilePhotoBufferRockhopper) | |
| plot(fiveFilePhotoModel1Rockhopper) | |
| plot(PC1~Count_1,data=fiveFilePhotoBufferRockhopper) | |
| shapiro.test(resid(fiveFilePhotoModel1Rockhopper)) | |
| fiveFilePhotoModelNullRockhopper<-lmer(PC1~Site+(1|Date),fiveFilePhotoBufferRockhopper) | |
| anova(fiveFilePhotoModelNullRockhopper,fiveFilePhotoModel1Rockhopper) | |
| r2beta(fiveFilePhotoModel1Rockhopper) | |
| # Bootstrap to generate vectors of F statistics and semi-partial R-squareds for each breeding stage for Magellanics (excluding IDLE data) | |
| oneFilePhotoBufferMagellanicNoIDLE<-oneFilePhotoBufferMagellanic[which(oneFilePhotoBufferMagellanic$Site!="PM_01"&oneFilePhotoBufferMagellanic$Site!="PM_02"),] | |
| allowableSampleLength<-min(sapply(unique(oneFilePhotoBufferMagellanicNoIDLE[,4]),function(x) dim(oneFilePhotoBufferMagellanicNoIDLE[which(oneFilePhotoBufferMagellanicNoIDLE$Breeding_Stage==x),])[1])) | |
| bootstrappedIMPhotoTests<-data.frame(Breeding_Stage=ordered(rep(unique(oneFilePhotoBufferMagellanicNoIDLE$Breeding_Stage),each=1000),levels=c("pair formation","incubation","early care","late care","molting")),F_Statistic=NA,R_Squared=NA,Message=NA) | |
| bootstrapFunction<-function(x) | |
| { | |
| stageRestrictedData<-oneFilePhotoBufferMagellanicNoIDLE[which(oneFilePhotoBufferMagellanicNoIDLE$Breeding_Stage==i),] | |
| sampledRows<-sample(1:allowableSampleLength,replace=TRUE) | |
| bootstrappedSample<-stageRestrictedData[sampledRows,] | |
| if(length(unique(bootstrappedSample$Site))<2|length(unique(bootstrappedSample$Date))<2|length(unique(bootstrappedSample$PC1))<2|length(unique(bootstrappedSample$Count_4))<2) | |
| { | |
| message<-"inadequate variability in sample" | |
| bootstrappedF<-NA | |
| RSquared<-NA | |
| }else{ | |
| theoreticalNumberOfCoefficients<-length(unique(bootstrappedSample$Site))+length(unique(bootstrappedSample$Date))+1+1-2 | |
| bootstrappedModel<-try(lm(PC1~Count_4+Site+Date,bootstrappedSample[which(!is.na(bootstrappedSample$Count_4)),])) | |
| message<-ifelse(class(bootstrappedModel)=="lm","model ran","error") | |
| if(sum(!is.na(bootstrappedModel$coefficients))!=theoreticalNumberOfCoefficients) | |
| { | |
| message<-"term(s) dropped" | |
| bootstrappedF<-NA | |
| RSquared<-NA | |
| }else{ | |
| bootstrappedModelNull<-try(lm(PC1~Site+Date,bootstrappedSample[which(!is.na(bootstrappedSample$Count_4)),])) | |
| bootstrappedF<-anova(bootstrappedModelNull,bootstrappedModel)$`F`[2] | |
| RSquared<-r2beta(bootstrappedModel)[which(r2beta(bootstrappedModel)$Effect=="Count_4"),]$Rsq | |
| } | |
| } | |
| return(cbind(bootstrappedF,RSquared,message)) | |
| } | |
| set.seed(3) | |
| startTime<-Sys.time() | |
| for(i in levels(bootstrappedIMPhotoTests$Breeding_Stage)) | |
| { | |
| cluster<-makeCluster(6) | |
| clusterExport(cluster,c("oneFilePhotoBufferMagellanicNoIDLE","allowableSampleLength","i","r2beta")) | |
| bootstrapConditionResults<-parSapply(cluster,1:1000,bootstrapFunction) | |
| stopCluster(cluster) | |
| bootstrappedIMPhotoTests[which(bootstrappedIMPhotoTests$Breeding_Stage==i),][,2:4]<-t(bootstrapConditionResults) | |
| } | |
| bootstrappedIMPhotoTests$F_Statistic<-as.numeric(bootstrappedIMPhotoTests$F_Statistic) | |
| bootstrappedIMPhotoTests$R_Squared<-as.numeric(bootstrappedIMPhotoTests$R_Squared) | |
| bootstrappedIMPhotoTests$Message<-as.factor(bootstrappedIMPhotoTests$Message) | |
| Sys.time()-startTime | |
| summary(bootstrappedIMPhotoTests) | |
| # Test for differences in F-statistics and semi-partial R-squareds between breeding stages, | |
| # passing to permutational tests due to failure to meet assumptions | |
| FStatisticModel1<-lm(F_Statistic~Breeding_Stage,bootstrappedIMPhotoTests) | |
| plot(FStatisticModel1) | |
| # | |
| # | |
| # | |
| # | |
| FStatisticModel2<-lm(sqrt(bootstrappedIMPhotoTests$F_Statistic+1)~Breeding_Stage,bootstrappedIMPhotoTests) | |
| plot(FStatisticModel2) | |
| # | |
| # | |
| # | |
| # | |
| FStatisticModel3<-lm(log(bootstrappedIMPhotoTests$F_Statistic+1)~Breeding_Stage,bootstrappedIMPhotoTests) | |
| plot(FStatisticModel3) | |
| # | |
| # | |
| # | |
| # | |
| shapiro.test(resid(FStatisticModel3)) | |
| hist(bootstrappedIMPhotoTests$F_Statistic) | |
| FStatisticModel4<-aovp(F_Statistic~Breeding_Stage,bootstrappedIMPhotoTests[!is.na(bootstrappedIMPhotoTests$F_Statistic),],perm="Prob") | |
| summary(FStatisticModel4) | |
| TukeyHSD(FStatisticModel4) | |
| R2Model1<-aovp(R_Squared~Breeding_Stage,bootstrappedIMPhotoTests[!is.na(bootstrappedIMPhotoTests$F_Statistic),],perm="Prob") | |
| summary(R2Model1) | |
| TukeyHSD(R2Model1) | |
| ###################################################################### | |
| # Plot overall count-acoustic activity relationships | |
| ###################################################################### | |
| # Plot results | |
| grey90Value<-col2rgb("grey90")[1,1] | |
| grey50Value<-col2rgb("grey50")[1,1] | |
| transparentGrey90<-rgb(grey90Value,grey90Value,grey90Value,alpha=50,maxColorValue=255) | |
| transparentGrey50<-rgb(grey50Value,grey50Value,grey50Value,alpha=100,maxColorValue=255) | |
| pdf("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Figures/Photo_Count-Acoustic_Activity_Relationships_Points_by_Colony.pdf",width=15,height=4,pointsize=6) | |
| set.seed(1) | |
| par(mfcol=c(1,3),mar=c(12,14,4,2)) | |
| plot(1,1,xlim=c(0,5),ylim=c(-0.15,0.15),type="n",axes=FALSE,xlab="",ylab="") | |
| axis(side=1,at=0:5,labels=0:5,col="grey30",col.axis="grey30",cex.axis=3,lwd=0,tck=0.02,lwd.ticks=1,line=1,padj=0.5,las=1) | |
| axis(side=2,at=seq(-0.15,0.15,0.05),labels=round(seq(-0.15,0.15,0.05),2),col="grey30",col.axis="grey30",cex.axis=3,lwd=0,tck=0.02,lwd.ticks=1,line=1,padj=0.5,las=1) | |
| box(col="grey30") | |
| title(xlab="2-meter Magellanic camera trap count",col.lab="grey30",cex.lab=3,line=7) | |
| title(ylab="acoustic activity",col.lab="grey30",cex.lab=3,line=10) | |
| oneFilePhotoBufferMagellanicShuffled<-oneFilePhotoBufferMagellanic[sample(1:nrow(oneFilePhotoBufferMagellanic)),] | |
| oneFilePhotoBufferMagellanicPoints<-ifelse(grepl("PM",oneFilePhotoBufferMagellanicShuffled$Site),24,21) | |
| points(jitter(oneFilePhotoBufferMagellanicShuffled$Count_2,amount=0.05),oneFilePhotoBufferMagellanicShuffled$PC1,pch=oneFilePhotoBufferMagellanicPoints,cex=2.5,col=transparentGrey50,bg=transparentGrey90) | |
| lines(x=predictorEffect("Count_2",oneFilePhotoModel1Magellanic2)$variables$Count_2$levels,y=predictorEffect("Count_2",oneFilePhotoModel1Magellanic2)$fit,col="grey30",lwd=2) | |
| lines(x=predictorEffect("Count_2",oneFilePhotoModel1Magellanic2)$variables$Count_2$levels,y=predictorEffect("Count_2",oneFilePhotoModel1Magellanic2)$lower,col="grey60",lwd=2,lty=2) | |
| lines(x=predictorEffect("Count_2",oneFilePhotoModel1Magellanic2)$variables$Count_2$levels,y=predictorEffect("Count_2",oneFilePhotoModel1Magellanic2)$upper,col="grey60",lwd=2,lty=2) | |
| text(0.1,0.15,"A.",col="grey30",cex=4,pos=1) | |
| legend(5.1, -0.1,xjust=1,yjust=1,c("Martillo Island","Staten Island"),border="grey30",pch=c(21,24),col=transparentGrey50,pt.bg=transparentGrey90,box.col="grey30",cex=2.5) | |
| plot(1,1,xlim=c(0,12),ylim=c(-0.15,0.15),type="n",axes=FALSE,xlab="",ylab="") | |
| axis(side=1,at=seq(0,12,2),labels=seq(0,12,2),col="grey30",col.axis="grey30",cex.axis=3,lwd=0,tck=0.02,lwd.ticks=1,line=1,padj=0.5,las=1) | |
| axis(side=2,at=seq(-0.15,0.15,0.05),labels=round(seq(-0.15,0.15,0.05),2),col="grey30",col.axis="grey30",cex.axis=3,lwd=0,tck=0.02,lwd.ticks=1,line=1,padj=0.5,las=1) | |
| box(col="grey30") | |
| title(xlab="4-meter Magellanic camera trap count",col.lab="grey30",cex.lab=3,line=7) | |
| title(ylab="acoustic activity",col.lab="grey30",cex.lab=3,line=10) | |
| points(jitter(oneFilePhotoBufferMagellanic$Count_4,amount=0.05*12/5),oneFilePhotoBufferMagellanic$PC1,pch=21,cex=2.5,col=transparentGrey50,bg=transparentGrey90) | |
| lines(x=predictorEffect("Count_4",oneFilePhotoModel1Magellanic4)$variables$Count_4$levels,y=predictorEffect("Count_4",oneFilePhotoModel1Magellanic4)$fit,col="grey30",lwd=2) | |
| lines(x=predictorEffect("Count_4",oneFilePhotoModel1Magellanic4)$variables$Count_4$levels,y=predictorEffect("Count_4",oneFilePhotoModel1Magellanic4)$lower,col="grey60",lwd=2,lty=2) | |
| lines(x=predictorEffect("Count_4",oneFilePhotoModel1Magellanic4)$variables$Count_4$levels,y=predictorEffect("Count_4",oneFilePhotoModel1Magellanic4)$upper,col="grey60",lwd=2,lty=2) | |
| text(0.1*12/5,0.15,"B.",col="grey30",cex=4,pos=1) | |
| plot(1,1,xlim=c(0,14),ylim=c(-0.4,0.4),type="n",axes=FALSE,xlab="",ylab="") | |
| axis(side=1,at=seq(0,14,2),labels=seq(0,14,2),col="grey30",col.axis="grey30",cex.axis=3,lwd=0,tck=0.02,lwd.ticks=1,line=1,padj=0.5,las=1) | |
| axis(side=2,at=seq(-0.4,0.4,0.1),labels=round(seq(-0.4,0.4,0.1),2),col="grey30",col.axis="grey30",cex.axis=3,lwd=0,tck=0.02,lwd.ticks=1,line=1,padj=0.5,las=1) | |
| box(col="grey30") | |
| title(xlab="1-meter Rockhopper camera trap count",col.lab="grey30",cex.lab=3,line=7) | |
| title(ylab="acoustic activity",col.lab="grey30",cex.lab=3,line=8) | |
| points(jitter(oneFilePhotoBufferRockhopper$Count_1,amount=0.05*14/5),oneFilePhotoBufferRockhopper$PC1,pch=24,cex=2.5,col=transparentGrey50,bg=transparentGrey90) | |
| lines(x=predictorEffect("Count_1",oneFilePhotoModel2Rockhopper)$variables$Count_1$levels,y=predictorEffect("Count_1",oneFilePhotoModel2Rockhopper)$fit,col="grey30",lwd=2) | |
| lines(x=predictorEffect("Count_1",oneFilePhotoModel2Rockhopper)$variables$Count_1$levels,y=predictorEffect("Count_1",oneFilePhotoModel2Rockhopper)$lower,col="grey60",lwd=2,lty=2) | |
| lines(x=predictorEffect("Count_1",oneFilePhotoModel2Rockhopper)$variables$Count_1$levels,y=predictorEffect("Count_1",oneFilePhotoModel2Rockhopper)$upper,col="grey60",lwd=2,lty=2) | |
| text(0.1*14/5,0.4,"C.",col="grey30",cex=4,pos=1) | |
| dev.off() | |
| ###################################################################### | |
| # Plot breeding stage differences in F-statistics and semi-partial R-squareds | |
| ###################################################################### | |
| # Extract distributional information | |
| FQuartiles<-sapply(split(bootstrappedIMPhotoTests[!is.na(bootstrappedIMPhotoTests$F_Statistic),],bootstrappedIMPhotoTests$Breeding_Stage[!is.na(bootstrappedIMPhotoTests$F_Statistic)]),function(x) summary(x$F_Statistic)) | |
| R2Quartiles<-sapply(split(bootstrappedIMPhotoTests[!is.na(bootstrappedIMPhotoTests$F_Statistic),],bootstrappedIMPhotoTests$Breeding_Stage[!is.na(bootstrappedIMPhotoTests$F_Statistic)]),function(x) summary(x$R_Squared)) | |
| contrastLetters<-c("a","a","b","c","d") | |
| # Plot results | |
| pdf("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Figures/Breeding_Stage_Fs_and_R2s.pdf",width=5,height=4,pointsize=6) | |
| par(mfcol=c(2,1),mar=c(4,6,2,2)) | |
| plot(1,1,xlim=c(0.5,5.5),ylim=c(0,85),type="n",axes=FALSE,xlab="",ylab="") | |
| for(i in 1:5) polygon(x=c((i-0.4),(i+0.4),(i+0.4),(i-0.4)),y=c(rep(FQuartiles[2,i],2),rep(FQuartiles[5,i],2)),border="grey30",lwd=1) | |
| for(i in 1:5) lines(x=c((i-0.4),(i+0.4)),y=c(rep(FQuartiles[3,i],2)),col="grey30",lwd=3,lend=1) | |
| for(i in 1:5) lines(x=rep(i,2),y=c(FQuartiles[5,i],FQuartiles[6,i]),col="grey30",lwd=1) | |
| for(i in 1:5) lines(x=rep(i,2),y=c(FQuartiles[2,i],FQuartiles[1,i]),col="grey30",lwd=1) | |
| for(i in 1:5) text(i,FQuartiles[6,i]+8.5,contrastLetters[i],col="grey30",cex=1.5) | |
| axis(side=1,at=1:5,labels=c("pair formation","incubation","early chick\nrearing","late chick\nrearing","molting"),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=0,line=-2,padj=1,las=1) | |
| axis(side=2,at=seq(0,80,10),labels=seq(0,80,10),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-0.75,padj=0.5,las=1) | |
| polygon(x=c(0.4,5.6,5.6,0.4),y=c(0,0,85,85),border="grey30",col=NA) | |
| mtext(expression(paste(italic("F"),"-statistic")),side=2,line=3,col="grey30",cex=2) | |
| text(0.6,76.5,"A.",col="grey30",cex=2) | |
| par(mar=c(6,6,0,2)) | |
| plot(1,1,xlim=c(0.5,5.5),ylim=c(0,0.7),type="n",axes=FALSE,xlab="",ylab="") | |
| for(i in 1:5) polygon(x=c((i-0.4),(i+0.4),(i+0.4),(i-0.4)),y=c(rep(R2Quartiles[2,i],2),rep(R2Quartiles[5,i],2)),border="grey30",lwd=1) | |
| for(i in 1:5) lines(x=c((i-0.4),(i+0.4)),y=c(rep(R2Quartiles[3,i],2)),col="grey30",lwd=3,lend=1) | |
| for(i in 1:5) lines(x=rep(i,2),y=c(R2Quartiles[5,i],R2Quartiles[6,i]),col="grey30",lwd=1) | |
| for(i in 1:5) lines(x=rep(i,2),y=c(R2Quartiles[2,i],R2Quartiles[1,i]),col="grey30",lwd=1) | |
| for(i in 1:5) text(i,R2Quartiles[6,i]+0.06,contrastLetters[i],col="grey30",cex=1.5) | |
| axis(side=1,at=1:5,labels=c("pair formation","incubation","early chick\nrearing","late chick\nrearing","molting"),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=0,line=-2,padj=1,las=1) | |
| axis(side=2,at=seq(0,0.7,0.1),labels=seq(0,0.7,0.1),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-0.75,padj=0.5,las=1) | |
| polygon(x=c(0.4,5.6,5.6,0.4),y=c(0,0,0.7,0.7),border="grey30",col=NA) | |
| mtext(expression("semi-partial"~italic("R")^"2"),side=2,line=3,col="grey30",cex=2) | |
| mtext("breeding stage",side=1,line=-2,col="grey30",cex=2,outer=TRUE,adj=0.55) | |
| text(0.6,0.54,"B.",col="grey30",cex=2) | |
| dev.off() | |
| ###################################################################### | |
| # Test relationships between population counts at varying distances and acoustic activity for various breeding stages at varying times | |
| ###################################################################### | |
| # Create container template for results and assign appropriate dates for breeding periods | |
| uniqueTimes<-acousticMetricsMagellanic$Time[1:288] | |
| testDistances<-1:125 | |
| populationAcousticRelationshipsResultsTemplate<-matrix(NA,125,288,dimnames=list(testDistances,uniqueTimes)) | |
| validCameraTrapDatesAndBreedingStages<-unique(cameraTrapCountsMagellanic[,3:4]) | |
| validPopulationDatesAndBreedingStages<-rbind(validCameraTrapDatesAndBreedingStages[which(validCameraTrapDatesAndBreedingStages$Breeding_Stage!="late care"),], | |
| data.frame(Date=c(20190127:20190131,20190201:20190205),Breeding_Stage=rep("late care",10))) | |
| # Create function to calculate Pearson's Rs and Spearman Ps | |
| populationAcousticRelationshipsCorrelationFunction<-function(populationData,acousticData,breedingStage) | |
| { | |
| relevantDates<-validPopulationDatesAndBreedingStages$Date[which(validPopulationDatesAndBreedingStages$Breeding_Stage==breedingStage)] | |
| relevantAcousticData<-arrange(acousticData[acousticData$Date%in%relevantDates,],Site) | |
| outputPs<-populationAcousticRelationshipsResultsTemplate | |
| outputRs<-populationAcousticRelationshipsResultsTemplate | |
| for(i in testDistances) | |
| { | |
| relevantDensities<-populationData[,(i+2)] | |
| for(j in uniqueTimes) | |
| { | |
| relevantAcousticDataAtTime<-relevantAcousticData[which(relevantAcousticData$Time==j&!is.na(relevantAcousticData$PC1)),] | |
| relevantPC1Averages<-sapply(split(relevantAcousticDataAtTime,relevantAcousticDataAtTime$Site),function(x) mean(x$PC1)) | |
| outputPs[i,j]<-cor.test(relevantPC1Averages,relevantDensities,method="spearman")$p.value | |
| outputRs[i,j]<-cor.test(relevantPC1Averages,relevantDensities,method="pearson")$estimate | |
| } | |
| } | |
| return(list(outputPs,outputRs)) | |
| } | |
| # Run tests | |
| earlyCareSitesMagellanic<-unique(acousticMetricsMagellanic$Site)[-c(1:3,6)] | |
| lateCareSitesMagellanic<-unique(acousticMetricsMagellanic$Site)[-(4:5)] | |
| moltingSitesMagellanic<-unique(acousticMetricsMagellanic$Site)[-(4:5)] | |
| populationAcousticRelationshipsRockhopperResults<-populationAcousticRelationshipsCorrelationFunction(populationDensitiesRockhopper,acousticMetricsRockhopper,"early care") | |
| populationAcousticRelationshipsMagellanicResultsEarlyCare<-populationAcousticRelationshipsCorrelationFunction(populationDensitiesMagellanic[-(3:6),],acousticMetricsMagellanic[acousticMetricsMagellanic$Site%in%earlyCareSitesMagellanic,],"early care") | |
| populationAcousticRelationshipsMagellanicResultsLateCare<-populationAcousticRelationshipsCorrelationFunction(populationDensitiesMagellanic[-(1:2),],acousticMetricsMagellanic[acousticMetricsMagellanic$Site%in%lateCareSitesMagellanic,],"late care") | |
| populationAcousticRelationshipsMagellanicResultsMolting<-populationAcousticRelationshipsCorrelationFunction(populationDensitiesMagellanic[-(1:2),],acousticMetricsMagellanic[acousticMetricsMagellanic$Site%in%moltingSitesMagellanic,],"molting") | |
| ###################################################################### | |
| # Plot relationships between population counts at varying distances and acoustic activity for various breeding stages at varying times | |
| ###################################################################### | |
| # Find regional sunrise and sunset times for the center date of each period | |
| centralCoordinatesAndDatesForPopulation<-data.frame(Latitude=c(mean(recorderLocations$Latitude[3:7]),mean(recorderLocations$Latitude[c(1:2,12:14)]),rep(mean(recorderLocations$Latitude[8:14]),2)), | |
| Longitude=c(mean(recorderLocations$Longitude[3:7]),mean(recorderLocations$Longitude[c(1:2,12:14)]),rep(mean(recorderLocations$Longitude[8:14]),2)), | |
| Date=c(rep("20181212",2),"20190131","20190326")) | |
| sunriseHoursForPopulation<-sapply(1:4,function(x) as.numeric(substr(format(round(getSunlightTimes(date=as.Date(centralCoordinatesAndDatesForPopulation$Date[x],"%Y%m%d"),lat=centralCoordinatesAndDatesForPopulation$Latitude[x],lon=centralCoordinatesAndDatesForPopulation$Longitude[x],tz="America/Argentina/Ushuaia",keep=c("sunrise"))$sunrise,units="mins"),format="%H%M"),1,2))) | |
| sunriseMinutesForPopulation<-sapply(1:4,function(x) as.numeric(substr(format(round(getSunlightTimes(date=as.Date(centralCoordinatesAndDatesForPopulation$Date[x],"%Y%m%d"),lat=centralCoordinatesAndDatesForPopulation$Latitude[x],lon=centralCoordinatesAndDatesForPopulation$Longitude[x],tz="America/Argentina/Ushuaia",keep=c("sunrise"))$sunrise,units="mins"),format="%H%M"),3,4))) | |
| sunriseTimesInMinutesForPopulation<-sunriseHoursForPopulation*60+sunriseMinutesForPopulation | |
| sunrisePointsForPopulation<-288*sunriseTimesInMinutesForPopulation/(60*24) | |
| sunsetHoursForPopulation<-sapply(1:4,function(x) as.numeric(substr(format(round(getSunlightTimes(date=as.Date(centralCoordinatesAndDatesForPopulation$Date[x],"%Y%m%d"),lat=centralCoordinatesAndDatesForPopulation$Latitude[x],lon=centralCoordinatesAndDatesForPopulation$Longitude[x],tz="America/Argentina/Ushuaia",keep=c("sunset"))$sunset,units="mins"),format="%H%M"),1,2))) | |
| sunsetMinutesForPopulation<-sapply(1:4,function(x) as.numeric(substr(format(round(getSunlightTimes(date=as.Date(centralCoordinatesAndDatesForPopulation$Date[x],"%Y%m%d"),lat=centralCoordinatesAndDatesForPopulation$Latitude[x],lon=centralCoordinatesAndDatesForPopulation$Longitude[x],tz="America/Argentina/Ushuaia",keep=c("sunset"))$sunset,units="mins"),format="%H%M"),3,4))) | |
| sunsetTimesInMinutesForPopulation<-sunsetHoursForPopulation*60+sunsetMinutesForPopulation | |
| sunsetPointsForPopulation<-288*sunsetTimesInMinutesForPopulation/(60*24) | |
| # Plot original figure | |
| png("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Figures/Population_Acoustic_Relationships_Original.png",units="in",width=6,height=8,pointsize=6,res=900) | |
| par(mfcol=c(4,1),mar=c(1,1,1,0),oma=c(6,6,2,14)) | |
| plot(1,1,xlim=c(0,288),ylim=c(0,125),type="n",axes=FALSE,xlab="",ylab="") | |
| clip(0.5,288.5,0.5,125.5) | |
| image(1:288,20:125,t(populationAcousticRelationshipsRockhopperResults[[2]])[,20:125],col=colorRampPalette(brewer.pal(11,"RdBu"))(200),breaks=seq(-1,1,0.01),add=TRUE) | |
| image(1:288,20:125,t(populationAcousticRelationshipsRockhopperResults[[1]])[,20:125]<=0.05,col=c(rgb(col2rgb("grey10")[1]/288,col2rgb("grey10")[2]/288,col2rgb("grey10")[3]/288,0.8),NA),add=TRUE) | |
| polygon(c(0.5,288.5,288.5,0.5),c(0.5,0.5,19.5,19.5),col="grey90",border=NA) | |
| lines(x=c(0,sunrisePointsForPopulation[1]),y=rep(125,2),col="grey90",lwd=8) | |
| lines(x=c(sunsetPointsForPopulation[1],288),y=rep(125,2),col="grey90",lwd=8) | |
| abline(v=c(73,145,217),col="grey60",lwd=0.5,lty=2) | |
| text(2.5,112.5,"A. rockhopper, early chick rearing",cex=2,col="grey90",pos=4) | |
| axis(side=1,at=seq(1,287,12),labels=seq(0,23),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-1,padj=0.5,las=1) | |
| axis(side=2,at=seq(10,125,10),labels=seq(10,125,10),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-2.6,padj=0.5,las=1) | |
| polygon(c(0.5,288.5,288.5,0.5),c(0.5,0.5,125.5,125.5),col=NA,border="grey30") | |
| for(j in 1:200) polygon(x=c(292,292,300,300),y=c((j)/200*125,(j+1)/200*125,(j+1)/200*125,j/200*125),col=colorRampPalette(brewer.pal(11,"RdBu"))(200)[j],border=NA,xpd=NA) | |
| polygon(x=c(292,292,300,300),y=c(0.5,125.5,125.5,0.5),col=NA,border="grey30",lwd=1,xpd=NA) | |
| text(x=317,y=62.5,labels="correlation\ncoefficient",cex=2,pos=4,col="grey30",xpd=NA) | |
| sapply(seq(-1,1,0.2),function(x) lines(c(300,302),rep((x+1)*125/2+0.5,2),col="grey30",xpd=NA)) | |
| sapply(seq(-1,1,0.2),function(x) text(302.5,(x+1)*125/2+0.5,x,col="grey30",cex=1.5,pos=4,xpd=NA)) | |
| plot(1,1,xlim=c(0,288),ylim=c(0,125),type="n",axes=FALSE,xlab="",ylab="") | |
| clip(0.5,288.5,0.5,125.5) | |
| image(1:288,1:20,t(populationAcousticRelationshipsMagellanicResultsEarlyCare[[2]])[,1:20],col=colorRampPalette(brewer.pal(11,"RdBu"))(200),breaks=seq(-1,1,0.01),add=TRUE) | |
| image(1:288,1:20,t(populationAcousticRelationshipsMagellanicResultsEarlyCare[[1]])[,1:20]<=0.05,col=c(rgb(col2rgb("grey10")[1]/288,col2rgb("grey10")[2]/288,col2rgb("grey10")[3]/288,0.8),NA),add=TRUE) | |
| polygon(c(0.5,288.5,288.5,0.5),c(20.5,20.5,125.5,125.5),col="grey90",border=NA) | |
| lines(x=c(0,sunrisePointsForPopulation[2]),y=rep(125,2),col="grey30",lwd=8) | |
| lines(x=c(sunsetPointsForPopulation[2],288),y=rep(125,2),col="grey30",lwd=8) | |
| abline(v=c(73,145,217),col="grey60",lwd=0.5,lty=2) | |
| text(2.5,112.5,"B. Magellanic, early chick rearing",cex=2,col="grey30",pos=4) | |
| axis(side=1,at=seq(1,287,12),labels=seq(0,23),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-1,padj=0.5,las=1) | |
| axis(side=2,at=seq(10,125,10),labels=seq(10,125,10),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-2.6,padj=0.5,las=1) | |
| polygon(c(0.5,288.5,288.5,0.5),c(0.5,0.5,125.5,125.5),col=NA,border="grey30") | |
| plot(1,1,xlim=c(0,288),ylim=c(0,125),type="n",axes=FALSE,xlab="",ylab="") | |
| clip(0.5,288.5,0.5,125.5) | |
| image(1:288,1:125,t(populationAcousticRelationshipsMagellanicResultsLateCare[[2]])[,1:125],col=colorRampPalette(brewer.pal(11,"RdBu"))(200),breaks=seq(-1,1,0.01),add=TRUE) | |
| image(1:288,1:125,t(populationAcousticRelationshipsMagellanicResultsLateCare[[1]])[,1:125]<=0.05,col=c(rgb(col2rgb("grey10")[1]/288,col2rgb("grey10")[2]/288,col2rgb("grey10")[3]/288,0.8),NA),add=TRUE) | |
| lines(x=c(0,sunrisePointsForPopulation[3]),y=rep(125,2),col="grey90",lwd=8) | |
| lines(x=c(sunsetPointsForPopulation[3],288),y=rep(125,2),col="grey90",lwd=8) | |
| abline(v=c(73,145,217),col="grey60",lwd=0.5,lty=2) | |
| text(2.5,112.5,"C. Magellanic, late chick rearing",cex=2,col="grey90",pos=4) | |
| axis(side=1,at=seq(1,287,12),labels=seq(0,23),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-1,padj=0.5,las=1) | |
| axis(side=2,at=seq(10,125,10),labels=seq(10,125,10),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-2.6,padj=0.5,las=1) | |
| polygon(c(0.5,288.5,288.5,0.5),c(0.5,0.5,125.5,125.5),col=NA,border="grey30") | |
| plot(1,1,xlim=c(0,288),ylim=c(0,125),type="n",axes=FALSE,xlab="",ylab="") | |
| clip(0.5,288.5,0.5,125.5) | |
| image(1:288,1:125,t(populationAcousticRelationshipsMagellanicResultsMolting[[2]])[,1:125],col=colorRampPalette(brewer.pal(11,"RdBu"))(200),breaks=seq(-1,1,0.01),add=TRUE) | |
| image(1:288,1:125,t(populationAcousticRelationshipsMagellanicResultsMolting[[1]])[,1:125]<=0.05,col=c(rgb(col2rgb("grey10")[1]/288,col2rgb("grey10")[2]/288,col2rgb("grey10")[3]/288,0.8),NA),add=TRUE) | |
| lines(x=c(0,sunrisePointsForPopulation[4]),y=rep(125,2),col="grey90",lwd=8) | |
| lines(x=c(sunsetPointsForPopulation[4],288),y=rep(125,2),col="grey90",lwd=8) | |
| abline(v=c(73,145,217),col="grey60",lwd=0.5,lty=2) | |
| text(2.5,112.5,"D. Magellanic, molting",cex=2,col="grey90",pos=4) | |
| axis(side=1,at=seq(1,287,12),labels=seq(0,23),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-1,padj=0.5,las=1) | |
| axis(side=2,at=seq(10,125,10),labels=seq(10,125,10),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-2.6,padj=0.5,las=1) | |
| polygon(c(0.5,288.5,288.5,0.5),c(0.5,0.5,125.5,125.5),col=NA,border="grey30") | |
| mtext("hour of day",side=1,line=3,col="grey30",cex=2,outer=TRUE) | |
| mtext("radius around recorder for density measurement (m)",side=2,line=1,col="grey30",cex=2,outer=TRUE) | |
| dev.off() | |
| # Plot figure with no transparent grey | |
| png("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Figures/Population_Acoustic_Relationships_No_Transparent.png",units="in",width=6,height=8,pointsize=6,res=900) | |
| par(mfcol=c(4,1),mar=c(1,1,1,0),oma=c(6,6,2,14)) | |
| plot(1,1,xlim=c(0,288),ylim=c(0,125),type="n",axes=FALSE,xlab="",ylab="") | |
| clip(0.5,288.5,0.5,125.5) | |
| image(1:288,1:20,t(populationAcousticRelationshipsMagellanicResultsEarlyCare[[2]])[,1:20],col=colorRampPalette(brewer.pal(11,"RdBu"))(200),breaks=seq(-1,1,0.01),add=TRUE) | |
| image(1:288,1:20,t(populationAcousticRelationshipsMagellanicResultsEarlyCare[[1]])[,1:20]<=0.05,col=c("grey60",NA),add=TRUE) | |
| polygon(c(0.5,288.5,288.5,0.5),c(20.5,20.5,125.5,125.5),col="grey90",border=NA) | |
| lines(x=c(0,sunrisePointsForPopulation[2]),y=rep(125,2),col="grey30",lwd=8) | |
| lines(x=c(sunsetPointsForPopulation[2],288),y=rep(125,2),col="grey30",lwd=8) | |
| abline(v=c(73,145,217),col="grey30",lwd=0.5,lty=2) | |
| text(2.5,112.5,"A. Magellanic, early chick rearing",cex=2,col="grey30",pos=4) | |
| axis(side=1,at=seq(1,287,12),labels=seq(0,23),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-1,padj=0.5,las=1) | |
| axis(side=2,at=seq(10,125,10),labels=seq(10,125,10),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-2.6,padj=0.5,las=1) | |
| polygon(c(0.5,288.5,288.5,0.5),c(0.5,0.5,125.5,125.5),col=NA,border="grey30") | |
| for(j in 1:200) polygon(x=c(292,292,300,300),y=c((j)/200*125,(j+1)/200*125,(j+1)/200*125,j/200*125),col=colorRampPalette(brewer.pal(11,"RdBu"))(200)[j],border=NA,xpd=NA) | |
| polygon(x=c(292,292,300,300),y=c(0.5,125.5,125.5,0.5),col=NA,border="grey30",lwd=1,xpd=NA) | |
| text(x=317,y=62.5,labels="correlation\ncoefficient",cex=2,pos=4,col="grey30",xpd=NA) | |
| sapply(seq(-1,1,0.2),function(x) lines(c(300,302),rep((x+1)*125/2+0.5,2),col="grey30",xpd=NA)) | |
| sapply(seq(-1,1,0.2),function(x) text(302.5,(x+1)*125/2+0.5,x,col="grey30",cex=1.5,pos=4,xpd=NA)) | |
| plot(1,1,xlim=c(0,288),ylim=c(0,125),type="n",axes=FALSE,xlab="",ylab="") | |
| clip(0.5,288.5,0.5,125.5) | |
| image(1:288,20:125,t(populationAcousticRelationshipsRockhopperResults[[2]])[,20:125],col=colorRampPalette(brewer.pal(11,"RdBu"))(200),breaks=seq(-1,1,0.01),add=TRUE) | |
| image(1:288,20:125,t(populationAcousticRelationshipsRockhopperResults[[1]])[,20:125]<=0.05,col=c("grey60",NA),add=TRUE) | |
| polygon(c(0.5,288.5,288.5,0.5),c(0.5,0.5,19.5,19.5),col="grey90",border=NA) | |
| lines(x=c(0,sunrisePointsForPopulation[1]),y=rep(125,2),col="grey90",lwd=8) | |
| lines(x=c(sunsetPointsForPopulation[1],288),y=rep(125,2),col="grey90",lwd=8) | |
| abline(v=c(73,145,217),col="grey30",lwd=0.5,lty=2) | |
| text(2.5,112.5,"B. rockhopper, early chick rearing",cex=2,col="grey90",pos=4) | |
| axis(side=1,at=seq(1,287,12),labels=seq(0,23),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-1,padj=0.5,las=1) | |
| axis(side=2,at=seq(10,125,10),labels=seq(10,125,10),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-2.6,padj=0.5,las=1) | |
| polygon(c(0.5,288.5,288.5,0.5),c(0.5,0.5,125.5,125.5),col=NA,border="grey30") | |
| plot(1,1,xlim=c(0,288),ylim=c(0,125),type="n",axes=FALSE,xlab="",ylab="") | |
| clip(0.5,288.5,0.5,125.5) | |
| image(1:288,1:125,t(populationAcousticRelationshipsMagellanicResultsLateCare[[2]])[,1:125],col=colorRampPalette(brewer.pal(11,"RdBu"))(200),breaks=seq(-1,1,0.01),add=TRUE) | |
| image(1:288,1:125,t(populationAcousticRelationshipsMagellanicResultsLateCare[[1]])[,1:125]<=0.05,col=c("grey60",NA),add=TRUE) | |
| lines(x=c(0,sunrisePointsForPopulation[3]),y=rep(125,2),col="grey90",lwd=8) | |
| lines(x=c(sunsetPointsForPopulation[3],288),y=rep(125,2),col="grey90",lwd=8) | |
| abline(v=c(73,145,217),col="grey30",lwd=0.5,lty=2) | |
| text(2.5,112.5,"C. Magellanic, late chick rearing",cex=2,col="grey90",pos=4) | |
| axis(side=1,at=seq(1,287,12),labels=seq(0,23),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-1,padj=0.5,las=1) | |
| axis(side=2,at=seq(10,125,10),labels=seq(10,125,10),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-2.6,padj=0.5,las=1) | |
| polygon(c(0.5,288.5,288.5,0.5),c(0.5,0.5,125.5,125.5),col=NA,border="grey30") | |
| plot(1,1,xlim=c(0,288),ylim=c(0,125),type="n",axes=FALSE,xlab="",ylab="") | |
| clip(0.5,288.5,0.5,125.5) | |
| image(1:288,1:125,t(populationAcousticRelationshipsMagellanicResultsMolting[[2]])[,1:125],col=colorRampPalette(brewer.pal(11,"RdBu"))(200),breaks=seq(-1,1,0.01),add=TRUE) | |
| image(1:288,1:125,t(populationAcousticRelationshipsMagellanicResultsMolting[[1]])[,1:125]<=0.05,col=c("grey60",NA),add=TRUE) | |
| lines(x=c(0,sunrisePointsForPopulation[4]),y=rep(125,2),col="grey90",lwd=8) | |
| lines(x=c(sunsetPointsForPopulation[4],288),y=rep(125,2),col="grey90",lwd=8) | |
| abline(v=c(73,145,217),col="grey30",lwd=0.5,lty=2) | |
| text(2.5,112.5,"D. Magellanic, molting",cex=2,col="grey90",pos=4) | |
| axis(side=1,at=seq(1,287,12),labels=seq(0,23),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-1,padj=0.5,las=1) | |
| axis(side=2,at=seq(10,125,10),labels=seq(10,125,10),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-2.6,padj=0.5,las=1) | |
| polygon(c(0.5,288.5,288.5,0.5),c(0.5,0.5,125.5,125.5),col=NA,border="grey30") | |
| mtext("Hour of day",side=1,line=3,col="grey30",cex=2,outer=TRUE) | |
| mtext("Radius around recorder for density measurement (m)",side=2,line=1,col="grey30",cex=2,outer=TRUE) | |
| dev.off() | |
| # Plot figure so functional in greyscale | |
| png("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Figures/Population_Acoustic_Relationships_No_Transparent_for_Greyscale.png",units="in",width=6,height=8,pointsize=6,res=900) | |
| par(mfcol=c(4,1),mar=c(1,1,1,0),oma=c(6,6,2,14)) | |
| plot(1,1,xlim=c(0,288),ylim=c(0,125),type="n",axes=FALSE,xlab="",ylab="") | |
| clip(0.5,288.5,0.5,125.5) | |
| image(1:288,1:20,t(populationAcousticRelationshipsMagellanicResultsEarlyCare[[2]])[,1:20],col=colorRampPalette(brewer.pal(11,"RdBu"))(200),breaks=seq(-1,1,0.01),add=TRUE) | |
| image(1:288,1:20,t(populationAcousticRelationshipsMagellanicResultsEarlyCare[[1]])[,1:20]<=0.05,col=c("white",NA),add=TRUE) | |
| patternFill1<-data.frame(X=rep(c(seq(0,288,10),seq(5,288,10)),11),Y=rep(seq(21,126,5),each=29),Text=".") | |
| text(patternFill1$X,patternFill1$Y,patternFill1$Text,col="grey60",cex=4) | |
| abline(v=c(73,145,217),col="grey30",lwd=0.5,lty=2) | |
| lines(x=c(0,sunrisePointsForPopulation[2]),y=rep(125,2),col="grey60",lwd=8) | |
| lines(x=c(sunsetPointsForPopulation[2],288),y=rep(125,2),col="grey60",lwd=8) | |
| brackets(114,3,222,3,h=3,col="grey30") | |
| brackets(212,11,123,11,h=3,col="grey30") | |
| legend(2.5,118,"A. Magellanic, early chick rearing",cex=2,adj=0.05,col="grey30",box.col=NA,bg=rgb(153,153,153,200,maxColorValue=255)) | |
| axis(side=1,at=seq(1,287,12),labels=seq(0,23),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-1,padj=0.5,las=1) | |
| axis(side=2,at=seq(10,125,10),labels=seq(10,125,10),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-2.6,padj=0.5,las=1) | |
| polygon(c(0.5,288.5,288.5,0.5),c(0.5,0.5,125.5,125.5),col=NA,border="grey30") | |
| for(j in 1:200) polygon(x=c(292,292,300,300),y=c((j)/200*125,(j+1)/200*125,(j+1)/200*125,j/200*125),col=colorRampPalette(brewer.pal(11,"RdBu"))(200)[j],border=NA,xpd=NA) | |
| polygon(x=c(292,292,300,300),y=c(0.5,125.5,125.5,0.5),col=NA,border="grey30",lwd=1,xpd=NA) | |
| text(x=317,y=62.5,labels="correlation\ncoefficient",cex=2,pos=4,col="grey30",xpd=NA) | |
| sapply(seq(-1,1,0.2),function(x) lines(c(300,302),rep((x+1)*125/2+0.5,2),col="grey30",xpd=NA)) | |
| sapply(seq(-1,1,0.2),function(x) text(302.5,(x+1)*125/2+0.5,x,col="grey30",cex=1.5,pos=4,xpd=NA)) | |
| plot(1,1,xlim=c(0,288),ylim=c(0,125),type="n",axes=FALSE,xlab="",ylab="") | |
| clip(0.5,288.5,0.5,125.5) | |
| image(1:288,20:125,t(populationAcousticRelationshipsRockhopperResults[[2]])[,20:125],col=colorRampPalette(brewer.pal(11,"RdBu"))(200),breaks=seq(-1,1,0.01),add=TRUE) | |
| image(1:288,20:125,t(populationAcousticRelationshipsRockhopperResults[[1]])[,20:125]<=0.05,col=c("white",NA),add=TRUE) | |
| patternFill2<-data.frame(X=rep(c(seq(0,288,10),seq(5,288,10)),4),Y=rep(seq(1,19,5),each=29),Text=".") | |
| text(patternFill2$X,patternFill2$Y,patternFill2$Text,col="grey60",cex=4) | |
| abline(v=c(73,145,217),col="grey30",lwd=0.5,lty=2) | |
| lines(x=c(0,sunrisePointsForPopulation[1]),y=rep(125,2),col="grey60",lwd=8) | |
| lines(x=c(sunsetPointsForPopulation[1],288),y=rep(125,2),col="grey60",lwd=8) | |
| legend(2.5,118,"B. Rockhopper, early chick rearing",cex=2,adj=0.05,col="grey30",box.col=NA,bg=rgb(153,153,153,200,maxColorValue=255)) | |
| axis(side=1,at=seq(1,287,12),labels=seq(0,23),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-1,padj=0.5,las=1) | |
| axis(side=2,at=seq(10,125,10),labels=seq(10,125,10),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-2.6,padj=0.5,las=1) | |
| polygon(c(0.5,288.5,288.5,0.5),c(0.5,0.5,125.5,125.5),col=NA,border="grey30") | |
| plot(1,1,xlim=c(0,288),ylim=c(0,125),type="n",axes=FALSE,xlab="",ylab="") | |
| clip(0.5,288.5,0.5,125.5) | |
| image(1:288,1:125,t(populationAcousticRelationshipsMagellanicResultsLateCare[[2]])[,1:125],col=colorRampPalette(brewer.pal(11,"RdBu"))(200),breaks=seq(-1,1,0.01),add=TRUE) | |
| image(1:288,1:125,t(populationAcousticRelationshipsMagellanicResultsLateCare[[1]])[,1:125]<=0.05,col=c("white",NA),add=TRUE) | |
| abline(v=c(73,145,217),col="grey30",lwd=0.5,lty=2) | |
| lines(x=c(0,sunrisePointsForPopulation[3]),y=rep(125,2),col="grey60",lwd=8) | |
| lines(x=c(sunsetPointsForPopulation[3],288),y=rep(125,2),col="grey60",lwd=8) | |
| legend(2.5,118,"C. Magellanic, late chick rearing",cex=2,adj=0.05,col="grey30",box.col=NA,bg=rgb(153,153,153,200,maxColorValue=255)) | |
| axis(side=1,at=seq(1,287,12),labels=seq(0,23),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-1,padj=0.5,las=1) | |
| axis(side=2,at=seq(10,125,10),labels=seq(10,125,10),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-2.6,padj=0.5,las=1) | |
| polygon(c(0.5,288.5,288.5,0.5),c(0.5,0.5,125.5,125.5),col=NA,border="grey30") | |
| plot(1,1,xlim=c(0,288),ylim=c(0,125),type="n",axes=FALSE,xlab="",ylab="") | |
| clip(0.5,288.5,0.5,125.5) | |
| image(1:288,1:125,t(populationAcousticRelationshipsMagellanicResultsMolting[[2]])[,1:125],col=colorRampPalette(brewer.pal(11,"RdBu"))(200),breaks=seq(-1,1,0.01),add=TRUE) | |
| image(1:288,1:125,t(populationAcousticRelationshipsMagellanicResultsMolting[[1]])[,1:125]<=0.05,col=c("white",NA),add=TRUE) | |
| abline(v=c(73,145,217),col="grey30",lwd=0.5,lty=2) | |
| lines(x=c(0,sunrisePointsForPopulation[4]),y=rep(125,2),col="grey60",lwd=8) | |
| lines(x=c(sunsetPointsForPopulation[4],288),y=rep(125,2),col="grey60",lwd=8) | |
| brackets(88,38,96,38,h=3,col="grey30") | |
| brackets(96,95,88,95,h=3,col="grey30") | |
| legend(2.5,118,"D. Magellanic, molting",cex=2,adj=0.08,col="grey30",box.col=NA,bg=rgb(153,153,153,200,maxColorValue=255)) | |
| axis(side=1,at=seq(1,287,12),labels=seq(0,23),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-1,padj=0.5,las=1) | |
| axis(side=2,at=seq(10,125,10),labels=seq(10,125,10),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-2.6,padj=0.5,las=1) | |
| polygon(c(0.5,288.5,288.5,0.5),c(0.5,0.5,125.5,125.5),col=NA,border="grey30") | |
| mtext("hour of day",side=1,line=3,col="grey30",cex=2,outer=TRUE) | |
| mtext("radius around recorder for density measurement (m)",side=2,line=1,col="grey30",cex=2,outer=TRUE) | |
| dev.off() | |
| # Plot figure with just early chick rearing and no transparent grey | |
| png("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Figures/Population_Acoustic_Relationships_No_Transparent_Early_Only.png",units="in",width=7,height=4,pointsize=6,res=900) | |
| par(mfcol=c(2,1),mar=c(1,1,1,0),oma=c(4,6,1,12)) | |
| plot(1,1,xlim=c(0,288),ylim=c(0,125),type="n",axes=FALSE,xlab="",ylab="") | |
| clip(0.5,288.5,0.5,125.5) | |
| image(1:288,1:20,t(populationAcousticRelationshipsMagellanicResultsEarlyCare[[2]])[,1:20],col=colorRampPalette(brewer.pal(11,"RdBu"))(200),breaks=seq(-1,1,0.01),add=TRUE) | |
| image(1:288,1:20,t(populationAcousticRelationshipsMagellanicResultsEarlyCare[[1]])[,1:20]<=0.05,col=c("grey60",NA),add=TRUE) | |
| polygon(c(0.5,288.5,288.5,0.5),c(20.5,20.5,125.5,125.5),col="grey90",border=NA) | |
| lines(x=c(0,sunrisePointsForPopulation[2]),y=rep(125,2),col="grey30",lwd=8) | |
| lines(x=c(sunsetPointsForPopulation[2],288),y=rep(125,2),col="grey30",lwd=8) | |
| abline(v=c(73,145,217),col="grey30",lwd=0.5,lty=2) | |
| text(2.5,112.5,"A. Magellanic, early chick rearing",cex=1.5,col="grey30",pos=4) | |
| axis(side=1,at=seq(1,287,12),labels=seq(0,23),col="grey30",col.axis="grey30",cex.axis=1,lwd=0,lwd.ticks=1,line=-0.6,padj=0,las=1) | |
| axis(side=2,at=seq(10,125,10),labels=seq(10,125,10),col="grey30",col.axis="grey30",cex.axis=1,lwd=0,lwd.ticks=1,line=-1.9,padj=0.5,las=1) | |
| polygon(c(0.5,288.5,288.5,0.5),c(0.5,0.5,125.5,125.5),col=NA,border="grey30") | |
| for(j in 1:200) polygon(x=c(292,292,300,300),y=c((j)/200*125,(j+1)/200*125,(j+1)/200*125,j/200*125),col=colorRampPalette(brewer.pal(11,"RdBu"))(200)[j],border=NA,xpd=NA) | |
| polygon(x=c(292,292,300,300),y=c(0.5,125.5,125.5,0.5),col=NA,border="grey30",lwd=1,xpd=NA) | |
| text(x=316,y=62.5,labels="correlation\ncoefficient",cex=1.5,pos=4,col="grey30",xpd=NA) | |
| sapply(seq(-1,1,0.2),function(x) lines(c(300,302),rep((x+1)*125/2+0.5,2),col="grey30",xpd=NA)) | |
| sapply(seq(-1,1,0.2),function(x) text(302.5,(x+1)*125/2+0.5,x,col="grey30",cex=1,pos=4,xpd=NA)) | |
| plot(1,1,xlim=c(0,288),ylim=c(0,125),type="n",axes=FALSE,xlab="",ylab="") | |
| clip(0.5,288.5,0.5,125.5) | |
| image(1:288,20:125,t(populationAcousticRelationshipsRockhopperResults[[2]])[,20:125],col=colorRampPalette(brewer.pal(11,"RdBu"))(200),breaks=seq(-1,1,0.01),add=TRUE) | |
| image(1:288,20:125,t(populationAcousticRelationshipsRockhopperResults[[1]])[,20:125]<=0.05,col=c("grey60",NA),add=TRUE) | |
| polygon(c(0.5,288.5,288.5,0.5),c(0.5,0.5,19.5,19.5),col="grey90",border=NA) | |
| lines(x=c(0,sunrisePointsForPopulation[1]),y=rep(125,2),col="grey30",lwd=8) | |
| lines(x=c(sunsetPointsForPopulation[1],288),y=rep(125,2),col="grey30",lwd=8) | |
| abline(v=c(73,145,217),col="grey30",lwd=0.5,lty=2) | |
| text(2.5,112.5,"B. rockhopper, early chick rearing",cex=1.5,col="grey30",pos=4) | |
| axis(side=1,at=seq(1,287,12),labels=seq(0,23),col="grey30",col.axis="grey30",cex.axis=1,lwd=0,lwd.ticks=1,line=-0.6,padj=0,las=1) | |
| axis(side=2,at=seq(10,125,10),labels=seq(10,125,10),col="grey30",col.axis="grey30",cex.axis=1,lwd=0,lwd.ticks=1,line=-1.9,padj=0.5,las=1) | |
| polygon(c(0.5,288.5,288.5,0.5),c(0.5,0.5,125.5,125.5),col=NA,border="grey30") | |
| mtext("hour of day",side=1,line=2,col="grey30",cex=2,outer=TRUE) | |
| mtext("radius around recorder for\ndensity measurement (m)",side=2,line=1,col="grey30",cex=2,outer=TRUE) | |
| dev.off() | |
| ###################################################################### | |
| # Create figure showing daily temporal dynamics | |
| ###################################################################### | |
| # Create dataset | |
| temporalDynamicsData<-rep(list(data.frame(Species=c(rep("Magellanic",3),"rockhopper",rep("Magellanic",2)),Breeding_Stage=c("pair formation","incubation",rep("early care",2),"late care","molting"), | |
| Mean_PC1=rep(NA,6),Standard_Error=rep(NA,6))),288) | |
| names(temporalDynamicsData)<-uniqueTimes | |
| for(i in 1:288) | |
| { | |
| for(j in 1:6) | |
| { | |
| if(temporalDynamicsData[[i]]$Species[j]=="Magellanic") speciesData<-acousticMetricsMagellanic | |
| if(temporalDynamicsData[[i]]$Species[j]=="rockhopper") speciesData<-acousticMetricsRockhopper | |
| relevantDates<-validPopulationDatesAndBreedingStages$Date[validPopulationDatesAndBreedingStages$Breeding_Stage==temporalDynamicsData[[i]]$Breeding_Stage[j]] | |
| speciesDataFromBreedingStage<-speciesData[speciesData$Date%in%relevantDates,] | |
| speciesDataFromBreedingStageAndTime<-speciesDataFromBreedingStage[speciesDataFromBreedingStage$Time==uniqueTimes[i],] | |
| siteMeans<-sapply(split(speciesDataFromBreedingStageAndTime,speciesDataFromBreedingStageAndTime$Site),function(x) mean(x$PC1,na.rm=TRUE)) | |
| temporalDynamicsData[[i]]$Mean_PC1[j]<-mean(siteMeans[!is.na(siteMeans)]) | |
| temporalDynamicsData[[i]]$Standard_Error[j]<-sd(siteMeans[!is.na(siteMeans)])/sqrt(length(siteMeans[!is.na(siteMeans)])) | |
| } | |
| } | |
| # Reshape dataset | |
| temporalDynamicsMeans<-matrix(NA,6,288,dimnames=list(apply(temporalDynamicsData[[1]],1,function(x) paste(x[1:2],collapse=" ")),uniqueTimes)) | |
| temporalDynamicsStandardErrors<-matrix(NA,6,288,dimnames=list(apply(temporalDynamicsData[[1]],1,function(x) paste(x[1:2],collapse=" ")),uniqueTimes)) | |
| for(i in 1:6) temporalDynamicsMeans[i,]<-sapply(temporalDynamicsData,function(x) x[i,3]) | |
| for(i in 1:6) temporalDynamicsStandardErrors[i,]<-sapply(temporalDynamicsData,function(x) x[i,4]) | |
| # Obtain sunrise and sunset times | |
| centralCoordinatesAndDatesForDailyDynamics<-data.frame(Latitude=c(rep(mean(recorderLocations$Latitude[12:14]),2),mean(recorderLocations$Latitude[c(1:2,12:14)]),mean(recorderLocations$Latitude[3:7]),rep(mean(recorderLocations$Latitude[8:14]),2)), | |
| Longitude=c(rep(mean(recorderLocations$Longitude[12:14]),2),mean(recorderLocations$Longitude[c(1:2,12:14)]),mean(recorderLocations$Longitude[3:7]),rep(mean(recorderLocations$Longitude[8:14]),2)), | |
| Date=c("20181016","20181029",rep("20181212",2),"20190131","20190326")) | |
| sunriseHoursForDailyDynamics<-sapply(1:6,function(x) as.numeric(substr(format(round(getSunlightTimes(date=as.Date(centralCoordinatesAndDatesForDailyDynamics$Date[x],"%Y%m%d"),lat=centralCoordinatesAndDatesForDailyDynamics$Latitude[x],lon=centralCoordinatesAndDatesForDailyDynamics$Longitude[x],tz="America/Argentina/Ushuaia",keep=c("sunrise"))$sunrise,units="mins"),format="%H%M"),1,2))) | |
| sunriseMinutesForDailyDynamics<-sapply(1:6,function(x) as.numeric(substr(format(round(getSunlightTimes(date=as.Date(centralCoordinatesAndDatesForDailyDynamics$Date[x],"%Y%m%d"),lat=centralCoordinatesAndDatesForDailyDynamics$Latitude[x],lon=centralCoordinatesAndDatesForDailyDynamics$Longitude[x],tz="America/Argentina/Ushuaia",keep=c("sunrise"))$sunrise,units="mins"),format="%H%M"),3,4))) | |
| sunriseTimesInMinutesForDailyDynamics<-sunriseHoursForDailyDynamics*60+sunriseMinutesForDailyDynamics | |
| sunrisePointsForDailyDynamics<-288*sunriseTimesInMinutesForDailyDynamics/(60*24) | |
| sunsetHoursForDailyDynamics<-sapply(1:6,function(x) as.numeric(substr(format(round(getSunlightTimes(date=as.Date(centralCoordinatesAndDatesForDailyDynamics$Date[x],"%Y%m%d"),lat=centralCoordinatesAndDatesForDailyDynamics$Latitude[x],lon=centralCoordinatesAndDatesForDailyDynamics$Longitude[x],tz="America/Argentina/Ushuaia",keep=c("sunset"))$sunset,units="mins"),format="%H%M"),1,2))) | |
| sunsetMinutesForDailyDynamics<-sapply(1:6,function(x) as.numeric(substr(format(round(getSunlightTimes(date=as.Date(centralCoordinatesAndDatesForDailyDynamics$Date[x],"%Y%m%d"),lat=centralCoordinatesAndDatesForDailyDynamics$Latitude[x],lon=centralCoordinatesAndDatesForDailyDynamics$Longitude[x],tz="America/Argentina/Ushuaia",keep=c("sunset"))$sunset,units="mins"),format="%H%M"),3,4))) | |
| sunsetTimesInMinutesForDailyDynamics<-sunsetHoursForDailyDynamics*60+sunsetMinutesForDailyDynamics | |
| sunsetPointsForDailyDynamics<-288*sunsetTimesInMinutesForDailyDynamics/(60*24) | |
| # Plot | |
| dailyDynamicsPlotLabels<-c("A. Magellanic, pair formation","B. Magellanic, incubation","C. Magellanic, early chick rearing","D. Rockhopper, early chick rearing","E. Magellanic, late chick rearing","F. Magellanic, molting") | |
| pdf("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Figures/Daily_Temporal_Dynamics.pdf",width=5,height=8,pointsize=6) | |
| par(mfcol=c(6,1),mar=c(2,2,2,2),oma=c(4,5,0,0)) | |
| for(i in 1:6) | |
| { | |
| if(i!=4) | |
| { | |
| plot(1,1,xlim=c(1,288),ylim=c(-0.12,0.12),type="n",axes=FALSE,xlab="",ylab="") | |
| clip(1,288,-0.12,0.12) | |
| polygon(c(1:288,288:1),c((temporalDynamicsMeans[i,]-temporalDynamicsStandardErrors[i,]),rev((temporalDynamicsMeans[i,]+temporalDynamicsStandardErrors[i,]))),col="grey90",border="grey90") | |
| abline(v=c(73,145,217),col="grey60",lwd=0.5,lty=2) | |
| lines(1:288,temporalDynamicsMeans[i,],col="grey30") | |
| lines(x=c(1,sunrisePointsForDailyDynamics[i]),y=rep(0.12,2),col="grey30",lwd=8) | |
| lines(x=c(sunsetPointsForDailyDynamics[i],288),y=rep(0.12,2),col="grey30",lwd=8) | |
| axis(side=1,at=seq(1,287,12),labels=seq(0,23),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-0.55,padj=0.5,las=1) | |
| axis(side=2,at=seq(-0.1,0.1,0.05),labels=seq(-0.1,0.1,0.05),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-2.45,padj=0.5,las=1) | |
| text(2.5,0.088,dailyDynamicsPlotLabels[i],cex=2,col="grey30",pos=4) | |
| polygon(c(1,288,288,1),c(rep(-0.12,2),rep(0.12,2)),col=NA,border="grey30") | |
| } | |
| if(i==4) | |
| { | |
| plot(1,1,xlim=c(1,288),ylim=c(-0.3,0.3),type="n",axes=FALSE,xlab="",ylab="") | |
| clip(1,288,-0.3,0.3) | |
| polygon(c(1:288,288:1),c((temporalDynamicsMeans[i,]-temporalDynamicsStandardErrors[i,]),rev((temporalDynamicsMeans[i,]+temporalDynamicsStandardErrors[i,]))),col="grey90",border="grey90") | |
| abline(v=c(73,145,217),col="grey60",lwd=0.5,lty=2) | |
| lines(1:288,temporalDynamicsMeans[i,],col="grey30") | |
| lines(x=c(1,sunrisePointsForDailyDynamics[i]),y=rep(0.3,2),col="grey30",lwd=8) | |
| lines(x=c(sunsetPointsForDailyDynamics[i],288),y=rep(0.3,2),col="grey30",lwd=8) | |
| axis(side=1,at=seq(1,287,12),labels=seq(0,23),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-0.55,padj=0.5,las=1) | |
| axis(side=2,at=seq(-0.3,0.3,0.1),labels=c(-0.3,-0.2,-0.1,0,0.1,0.2,0.3),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-2.45,padj=0.5,las=1) | |
| text(2.5,0.22,dailyDynamicsPlotLabels[i],cex=2,col="grey30",pos=4) | |
| polygon(c(1,288,288,1),c(rep(-0.3,2),rep(0.3,2)),col=NA,border="grey30") | |
| } | |
| } | |
| mtext("hour of day",side=1,line=2,col="grey30",cex=2,outer=TRUE) | |
| mtext("acoustic activity",side=2,line=1,col="grey30",cex=2,outer=TRUE) | |
| dev.off() | |
| # Plot just early chick rearing | |
| dailyDynamicsPlotLabelsEarlyOnly<-c(NA,NA,"A. Magellanic, early chick rearing","B. rockhopper, early chick rearing") | |
| pdf("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Figures/Daily_Temporal_Dynamics_Early_Only.pdf",width=8,height=4,pointsize=6) | |
| par(mfcol=c(2,1),mar=c(2,2,2,2),oma=c(4,5,0,0)) | |
| for(i in 3:4) | |
| { | |
| if(i!=4) | |
| { | |
| plot(1,1,xlim=c(1,288),ylim=c(-0.12,0.12),type="n",axes=FALSE,xlab="",ylab="") | |
| clip(1,288,-0.12,0.12) | |
| polygon(c(1:288,288:1),c((temporalDynamicsMeans[i,]-temporalDynamicsStandardErrors[i,]),rev((temporalDynamicsMeans[i,]+temporalDynamicsStandardErrors[i,]))),col="grey90",border="grey90") | |
| abline(v=c(73,145,217),col="grey60",lwd=0.5,lty=2) | |
| lines(1:288,temporalDynamicsMeans[i,],col="grey30") | |
| lines(x=c(1,sunrisePointsForDailyDynamics[i]),y=rep(0.12,2),col="grey30",lwd=8) | |
| lines(x=c(sunsetPointsForDailyDynamics[i],288),y=rep(0.12,2),col="grey30",lwd=8) | |
| axis(side=1,at=seq(1,287,12),labels=seq(0,23),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-0.55,padj=0.5,las=1) | |
| axis(side=2,at=seq(-0.1,0.1,0.05),labels=seq(-0.1,0.1,0.05),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-2.6,padj=0.5,las=1) | |
| text(2.5,0.088,dailyDynamicsPlotLabelsEarlyOnly[i],cex=2,col="grey30",pos=4) | |
| polygon(c(1,288,288,1),c(rep(-0.12,2),rep(0.12,2)),col=NA,border="grey30") | |
| } | |
| if(i==4) | |
| { | |
| plot(1,1,xlim=c(1,288),ylim=c(-0.3,0.3),type="n",axes=FALSE,xlab="",ylab="") | |
| clip(1,288,-0.3,0.3) | |
| polygon(c(1:288,288:1),c((temporalDynamicsMeans[i,]-temporalDynamicsStandardErrors[i,]),rev((temporalDynamicsMeans[i,]+temporalDynamicsStandardErrors[i,]))),col="grey90",border="grey90") | |
| abline(v=c(73,145,217),col="grey60",lwd=0.5,lty=2) | |
| lines(1:288,temporalDynamicsMeans[i,],col="grey30") | |
| lines(x=c(1,sunrisePointsForDailyDynamics[i]),y=rep(0.3,2),col="grey30",lwd=8) | |
| lines(x=c(sunsetPointsForDailyDynamics[i],288),y=rep(0.3,2),col="grey30",lwd=8) | |
| axis(side=1,at=seq(1,287,12),labels=seq(0,23),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-0.55,padj=0.5,las=1) | |
| axis(side=2,at=seq(-0.3,0.3,0.1),labels=c(-0.3,-0.2,-0.1,0,0.1,0.2,0.3),col="grey30",col.axis="grey30",cex.axis=1.5,lwd=0,lwd.ticks=1,line=-2.6,padj=0.5,las=1) | |
| text(2.5,0.22,dailyDynamicsPlotLabelsEarlyOnly[i],cex=2,col="grey30",pos=4) | |
| polygon(c(1,288,288,1),c(rep(-0.3,2),rep(0.3,2)),col=NA,border="grey30") | |
| } | |
| } | |
| mtext("hour of day",side=1,line=2,col="grey30",cex=2,outer=TRUE) | |
| mtext("acoustic activity",side=2,line=1,col="grey30",cex=2,outer=TRUE) | |
| dev.off() |