#Script for Communication biology manuscript: # Marine fish larvae consistently use external cues for orientation # Igal Berenshtein, Robin Faillettaz, Jean-Oliver Irisson, Moshe Kiflawi, Ulrike E. Siebeckg, Jeffery M. Leis, Claire B. Paris #######section 1: compute the CRW null for the CRW-VM method rm(list=ls())#clears the workspace graphics.off()#closes figures # setwd('/path/to/working/directory/') setwd('C:/Users/ixb260/Downloads/') #install the following libraries if not installed: library("circular") library("reshape2") library("akima") ###computing the expected by CRW #this code should be run based on the number of observations per trial (Nobs=21, 80, 180, 300), here we run only Nobs=300 kappMax=400 # maximal kappa value (von misses distribution concentration parameter) nStep=21 # number of observations per trial (Nobs=21, 80, 180, 300) nRep=1000 # number of permutations Rvec_traj=matrix(0,nrow=nRep,ncol =kappMax) # matrix for bearings R values of simulated trajectories Rvec_ta=matrix(0,nrow=nRep,ncol =kappMax)# matrix for R values of "angle changes" (turning angles) for simulated trajectories kappaEstHeading=matrix(0,nrow=nRep,ncol =kappMax)# matrix for von-Mises estimated kappa of bearings kappaTA=matrix(0,nrow=nRep,ncol =kappMax)# matrix for von-Mises estimated kappa of turning angles # This may take a few minuets to run for (j in 1:nRep){ for (i in 1:kappMax){ KAPPA=(i-1) ta=rvonmises(nStep, 0, KAPPA) # create virtual turning angles startheading = runif(1, min=0, max=(2*pi)) # initial heading (0, 2pi) headingsSwim <- cumsum(c(startheading, ta[-1]))%%(2*pi) # compute successive headings Rvec_traj[j,i]=as.numeric(rho.circular(headingsSwim))#compute R of virtual bearings Rvec_ta[j,i]=as.numeric(rho.circular(ta))#compute R of virtual turning angles z<-mle.vonmises(headingsSwim)#compute von-Mises kappa parameter of virtual bearings kappaEstHeading[j,i]=z$kappa kappaTA[j,i]=KAPPA#record VM turning angles kappa } print(j) } Data_21_step<-list(Rvec_ta,Rvec_traj,kappaEstHeading,kappaTA)#save all this data into a single list variable save(Data_21_step,file = "Data_Nobs_21_step.RData")#save so that this code will not need to be run over and over again ############################ load("Data_Nobs_21_step.RData") #load the data saved above ### define data frame that records orientation statistics per individual larva SuperMetaData=NULL SuperMetaData<-matrix(NA,nrow=1000,ncol=11) SuperMetaData<-as.data.frame(SuperMetaData) names(SuperMetaData)=c("species","method","originalFidhId","DepNum","MeanAz","Rc","R_TA","KappaHeading","est_Kapp_TA","Quant", "Ray_Pval") ##Caesio_cuning - example DataCaesio<-read.csv('Caesio_cuning.csv')# orientation trial bearings data indFish<-DataCaesio$Fish..[!duplicated(DataCaesio$Fish..)]#get vector of individual trials/larvae DataCaesioMeta<-matrix(NA,nrow=length(indFish),ncol=8)#create a dataframe per species DataCaesioMeta<-as.data.frame(DataCaesioMeta) names(DataCaesioMeta)=c("originalFidhId","DepNum","MeanAz","Rc","KappaHeading","est_Kapp_TA","Quant", "Ray_Pval") DataCaesioMeta$originalFidhId=indFish#assign individual trials/larvae as rows for (i in 1:length(indFish)){ #run through the trials and compute statistics at the individual level #CRW-von Mises method tempData=NULL # define data of the individual tempData$headings=DataCaesio$Bearing[which(DataCaesio$Fish..==indFish[i])] # bearing data of the individual if(length(tempData$headings)>20){ tempData$TA=NA tempData$TA[2:length(tempData$headings)]=diff(tempData$headings) # turning angles of the individual tempData$TA<-as.circular(tempData$TA, type="angles", units="degrees",template="geographics", zero=0, rotation="clock") tempData$headings<-as.circular(tempData$headings, type="angles", units="degrees",template="geographics", zero=0, rotation="clock") vm<-mle.vonmises(tempData$headings) #estimate Von-Mises dsitribution kappa parameter for bearings DataCaesioMeta$KappaHeading[i]<-vm$kappa ray<-rayleigh.test(tempData$headings)#compute Rayleigh's test for the bearings sequence DataCaesioMeta$Ray_Pval[i]=ray$p.value#record p-value of the Rayleigh's test vm_ta<-mle.vonmises(tempData$TA)#estimate Von-Mises dsitribution kappa parameter for turning angles DataCaesioMeta$est_Kapp_TA[i]<-vm_ta$kappa DataCaesioMeta$MeanAz[i]<-mean.circular(tempData$headings)#compute the mean bearing DataCaesioMeta$Rc[i]<-rho.circular(tempData$headings)#compute the R of the bearings DataCaesioMeta$R_TA[i]<-rho.circular(tempData$TA[2:length(tempData$headings)])#compute the R of the turning angles DataCaesioMeta$DepNum[i]=i#record the trial's number #### simulated vs. observed Rhos CRW-resampling method (CRW-r) newDirVec=NULL #temporary variable for virtual headings sequence based on resampled turning angles rhoVec=NULL #vector of Rs computed on virtual headings sequence based on resampled turning angles testTA<-tempData$TA #get the empirical turning angles of the trial for (iSim1 in 1:100){ samp=sample(testTA[2:length(testTA)],replace = FALSE)#resample the empirical turning angles without replacement newDirVec=cumsum(c(tempData$headings[1],samp))#compute virtual bearings sequence newDirVec=newDirVec %% 360#compute bearings sequence newDirVec<- as.circular(newDirVec, type="angles", #define variable as circular units="degrees",template="geographics", zero=0, rotation="clock") rhoVec[iSim1]=rho.circular(newDirVec,na.rm = T)#compute R of virtual bearings sequence } Qvec=quantile(rhoVec, probs = c(1:100, NA)/100)#compute the 1-100 quantiles of all Rs obtained from the virual bearings sequences Qvec=Qvec[!is.na(Qvec)]#remove NAs quantDel=abs(DataCaesioMeta$Rc[i]-Qvec)#find the quantile of the observed R quantDel<-quantDel[!is.na(quantDel)] Q<-which.min(quantDel) DataCaesioMeta$Quant[i]=Q[[1]]#find the quantile of the observed R quantDel2=abs(rhoVec[iSim1]-Qvec) #computing the the null distribution of quantiles under a randomly sampled sequence quantDel2<-quantDel2[!is.na(quantDel2)] Q2<-which.min(quantDel2) DataCaesioMeta$Quant_null[i]=Q2[[1]] #### } } #record all information in sepcies specific meta-dataframe DataCaesioMeta<-DataCaesioMeta[complete.cases(DataCaesioMeta$MeanAz),] DataCaesioMeta<-DataCaesioMeta[which(DataCaesioMeta$Ray_Pval<0.05),] #keep only directional trials #record all information in on general meta-dataframe SuperMetaData$species[1:length(DataCaesioMeta$Rc)]="Caesio_cuning" SuperMetaData$method[1:length(DataCaesioMeta$Rc)]="Following" SuperMetaData$originalFidhId[1:length(DataCaesioMeta$Rc)]=DataCaesioMeta$originalFidhId SuperMetaData$DepNum[1:length(DataCaesioMeta$Rc)]=DataCaesioMeta$DepNum SuperMetaData$MeanAz[1:length(DataCaesioMeta$Rc)]=DataCaesioMeta$MeanAz SuperMetaData$Rc[1:length(DataCaesioMeta$Rc)]=DataCaesioMeta$Rc SuperMetaData$R_TA[1:length(DataCaesioMeta$Rc)]=DataCaesioMeta$R_TA SuperMetaData$KappaHeading[1:length(DataCaesioMeta$Rc)]=DataCaesioMeta$KappaHeading SuperMetaData$est_Kapp_TA[1:length(DataCaesioMeta$Rc)]=DataCaesioMeta$est_Kapp_TA SuperMetaData$Quant[1:length(DataCaesioMeta$Rc)]=DataCaesioMeta$Quant SuperMetaData$Ray_Pval[1:length(DataCaesioMeta$Rc)]=DataCaesioMeta$Ray_Pval SuperMetaData$Quant_null[1:length(DataCaesioMeta$Rc)]=DataCaesioMeta$Quant_null ##########################Here we apply the same protocol for the rest of the species SuperMetaData<-SuperMetaData[complete.cases(SuperMetaData$Rc),]#after inserting information from all species we keep only the data and remove the empty lines ########## plotting following data############# load("Data_21_step.RData") #crw Rvec_ta_21<-Data_21_step[[1]] #load CRW-VM null Rho values of turning angles Rvec_traj_21<-Data_21_step[[2]] #load CRW-VM null Rho values of bearings SuperMetaData$Steps_Num=21 #Scuba-following method SuperMetaData$Steps_Num[which(SuperMetaData$method=="DISC")]=180#Main DISC SuperMetaData$Steps_Num[which(SuperMetaData$species=="Chromis_atripectoralis_ind")]=90#Chromis_Atrip.. irisson et al. 2015 c("Chromis_atripectoralis_ind","Chromis_atripectoralis_group") SuperMetaData$Steps_Num[which(SuperMetaData$species=="Chromis_atripectoralis_group")]=90#Chromis_Atrip.. irisson et al. 2015 c("Chromis_atripectoralis_ind","Chromis_atripectoralis_group") SuperMetaData$Steps_Num[which(SuperMetaData$species=="Premnas_biaculeatus")]=300#Berenshtein et al. 2012 ###computing the delta Rh (from CRW expected) fit_21<-smooth.spline(Rvec_ta_21,Rvec_traj_21,df=16) #16 degrees of freedom ########### #21 dim(Rvec_traj_21) Q_21=NULL#upper confidence interval band lower_21=NULL#lower confidence interval band mid_21=NULL#median confidence interval band mid_TA_21=NULL#vector of CRW-VM TA values (X-axis) Q_21_05=Q_21_95=Q_21_01=Q_21_1=Q_21_2=Q_21_3=Q_21_4=Q_21_99=Q_21_9=Q_21_8=Q_21_7=Q_21_6=Q_21_5=NULL #define CRW-VM quantile vectors for (i in 1:dim(Rvec_traj_21)[2]){ Q_21_99[i]=quantile(Rvec_traj_21[,i],0.99) #compute quantile vector based on the desired quantile Q_21_95[i]=quantile(Rvec_traj_21[,i],0.95) Q_21_9[i]=quantile(Rvec_traj_21[,i],0.90) Q_21_8[i]=quantile(Rvec_traj_21[,i],0.80) Q_21_7[i]=quantile(Rvec_traj_21[,i],0.70) Q_21_6[i]=quantile(Rvec_traj_21[,i],0.60) Q_21_5[i]=quantile(Rvec_traj_21[,i],0.50) Q_21_4[i]=quantile(Rvec_traj_21[,i],0.40) Q_21_3[i]=quantile(Rvec_traj_21[,i],0.30) Q_21_2[i]=quantile(Rvec_traj_21[,i],0.20) Q_21_1[i]=quantile(Rvec_traj_21[,i],0.10) Q_21_05[i]=quantile(Rvec_traj_21[,i],0.05) Q_21_01[i]=quantile(Rvec_traj_21[,i],0.01) lower_21[i]=quantile(Rvec_traj_21[,i],0.1) mid_TA_21[i]=quantile(Rvec_ta_21[,i],0.5) } for (i in 1:dim(Rvec_ta_21)[2]){#vector of CRW-VM TA values (X-axis) mid_TA_21[i]=quantile(Rvec_ta_21[,i],0.5) } #############################################applying Quant-CRW-VM ###########2D linear interpolation to obtain the CRW-VM quantiles #Nobs=21 smooth the quantile vectors fit_21_01<-smooth.spline(mid_TA_21,Q_21_01,df=16) #16 degrees of freedom fit_21_05<-smooth.spline(mid_TA_21,Q_21_05,df=16) #16 degrees of freedom fit_21_10<-smooth.spline(mid_TA_21,Q_21_1,df=16) fit_21_20<-smooth.spline(mid_TA_21,Q_21_2,df=16) fit_21_30<-smooth.spline(mid_TA_21,Q_21_3,df=16) fit_21_40<-smooth.spline(mid_TA_21,Q_21_4,df=16) fit_21_50<-smooth.spline(mid_TA_21,Q_21_5,df=16) fit_21_60<-smooth.spline(mid_TA_21,Q_21_6,df=16) fit_21_70<-smooth.spline(mid_TA_21,Q_21_7,df=16) fit_21_80<-smooth.spline(mid_TA_21,Q_21_8,df=16) fit_21_90<-smooth.spline(mid_TA_21,Q_21_9,df=16) fit_21_95<-smooth.spline(mid_TA_21,Q_21_95,df=16) fit_21_99<-smooth.spline(mid_TA_21,Q_21_99,df=16) data2d_interp<-matrix(NA,nrow=length(fit_21_01$x), ncol=15) #create matrix for 2-D interpolation data2d_interp[,1]=mid_TA_21 data2d_interp[,2]=fit_21_01$y data2d_interp[,3]=fit_21_10$y data2d_interp[,4]=fit_21_20$y data2d_interp[,5]=fit_21_30$y data2d_interp[,6]=fit_21_40$y data2d_interp[,7]=fit_21_50$y data2d_interp[,8]=fit_21_60$y data2d_interp[,9]=fit_21_70$y data2d_interp[,10]=fit_21_80$y data2d_interp[,11]=fit_21_90$y data2d_interp[,12]=fit_21_95$y data2d_interp[,13]=fit_21_99$y data2d_interp[,14]=seq(0.9,1,length.out = length(fit_21_99$y))#max 100% quantile data2d_interp[,15]=seq(0,0.01,length.out = length(fit_21_99$y))#min 1% quantile data2d_interp<-as.data.frame(data2d_interp) names(data2d_interp)=c("TA","01","10","20","30","40","50","60","70","80","90","95","99","100","0")#,"100","0" data2d_interp2=melt(data2d_interp,id.vars=c("TA")) #re-arrange the data so that it would be prepared for interpolation data2d_interp2$variable2<-as.numeric(as.character(data2d_interp2$variable)) names(data2d_interp2)=c("R_TA","Quant_name","R_c","Qaunt_num") ## linear interpolation meta.li <- interp(data2d_interp2$R_TA, data2d_interp2$R_c, data2d_interp2$Qaunt_num) li.zmin <- min(meta.li$z,na.rm=TRUE) li.zmax <- max(meta.li$z,na.rm=TRUE) breaks <- pretty(c(li.zmin,li.zmax),10) colors <- heat.colors(length(breaks)-1) meta<-data2d_interp2 ## increase smoothness (using finer grid): names(meta) meta.smooth <- with(meta, interp(R_TA, R_c, Qaunt_num, nx=100, ny=100)) si.zmin <- min(meta.smooth$z,na.rm=TRUE) si.zmax <- max(meta.smooth$z,na.rm=TRUE) breaks <- pretty(c(si.zmin,si.zmax),10) colors <- heat.colors(length(breaks)-1) #plot the intepolated data # image (meta.smooth, main = "interp(, *) on finer grid", # breaks=breaks, col=colors) # contour(meta.smooth, add = TRUE, levels=breaks, col = "thistle") # points(meta$R_TA,meta$R_c, pch = 3, cex = 2, col = "blue") #################################inserting the CRW-VM quantile data to the meta data frame SuperMetaData$Qunat_CRW_VM=NA for (i in 1:length(SuperMetaData$species)){ if (SuperMetaData$Steps_Num[i]==21){ meta_lin_temp <- with(meta,interp(R_TA,R_c,Qaunt_num,SuperMetaData$R_TA[i],SuperMetaData$Rc[i],linear = FALSE, extrap = TRUE)) SuperMetaData$Qunat_CRW_VM[i]=meta_lin_temp$z print(i) } } SuperMetaData$Qunat_CRW_VM[SuperMetaData$Qunat_CRW_VM>100]=100 #bound the quantiles between 1 and 100 SuperMetaData$Qunat_CRW_VM[SuperMetaData$Qunat_CRW_VM<1]=1 ###########manual correction of quantiles after visual inspection id=which(SuperMetaData$Qunat_CRW_VM==100) SuperMetaData$Qunat_CRW_VM[id]=75 id2=which(SuperMetaData$Qunat_CRW_VM==1) SuperMetaData$Qunat_CRW_VM[id2]=100 #########################plotting ##plotting CRW baseline SuperMetaDataFollowing=SuperMetaData[SuperMetaData$method=="Following",] #plotting crosses spVec=SuperMetaData$species[!duplicated(SuperMetaData$species)] #species vector spVec=spVec[complete.cases(spVec)] spVecFollowing=SuperMetaDataFollowing$species[!duplicated(SuperMetaDataFollowing$species)]#species vector for the following experiment spVecFollowing=spVecFollowing[complete.cases(spVecFollowing)]#remove NAs colorVec=c("#81B8D9", "#9E43E6", "#9A729C", "#78E0C9", "#DEB2D6", "#8486DE", "#D4DBCF", "#C5DA99", "#DA64D0", "#D6DF4F", "#D6A066","#E26474", "#70DF71") par(mfrow=c(1,2)) #############Quantile Quality Control and visual inspection plot(c(0,1),c(0,1),xlim=c(0.2,1),ylim=c(0.2,1),pch=20,cex=0.5,xlab="R_d_theta",ylab="R_theta") yy=colMeans(Rvec_traj_21) xx=colMeans(Rvec_ta_21) #plot the CRW-VM quantiles lines(fit_21_05$x,fit_21_05$y, col="grey",lwd=3) lines(fit_21_10$x,fit_21_10$y, col="grey") lines(fit_21_20$x,fit_21_20$y, col="grey") lines(fit_21_30$x,fit_21_30$y, col="grey") lines(fit_21_40$x,fit_21_40$y, col="grey") lines(fit_21_50$x,fit_21_50$y, col="grey") lines(fit_21_60$x,fit_21_60$y, col="grey") lines(fit_21_70$x,fit_21_70$y, col="grey") lines(fit_21_80$x,fit_21_80$y, col="grey") lines(fit_21_90$x,fit_21_90$y, col="grey") lines(fit_21_95$x,fit_21_95$y, col="grey",lwd=3) fit_21_mean<-smooth.spline(mid_TA_21,colMeans(Rvec_traj_21),df=16) #16 degrees of freedom lines(fit_21_mean$x,fit_21_mean$y, col="blue",lwd=4)#plot the smoothed vector of CRW-VM R_thetha means #delta Rc - ditance from CRW curve SP_DataFollowing<-matrix(0,nrow=length(spVecFollowing),ncol=5) SP_DataFollowing<-as.data.frame(SP_DataFollowing) names(SP_DataFollowing)<-c("Species","mean_R_TA","mean_Rc","Rc_expected_CRW","Delta_Rc") SP_DataFollowing$Species=spVecFollowing # rankVec<-(0.25+round(rank(SP_DataFollowing$mean_Rc)/length(SP_DataFollowing$mean_Rc),2))/(max((0.2+round(rank(SP_DataFollowing$mean_Rc)/length(SP_DataFollowing$mean_Rc),2)))) for (i in 1:length(spVecFollowing)){ #compute error bars based on 95% Confidence Intervals Following_TA_ci<- qnorm(0.975)*sd(SuperMetaDataFollowing$R_TA[which(SuperMetaDataFollowing$species==spVecFollowing[i])])/sqrt(length(SuperMetaDataFollowing$R_TA[which(SuperMetaDataFollowing$species==spVecFollowing[i])])) Following_Rc_ci<- qnorm(0.975)*sd(SuperMetaDataFollowing$Rc[which(SuperMetaDataFollowing$species==spVecFollowing[i])])/sqrt(length(SuperMetaDataFollowing$Rc[which(SuperMetaDataFollowing$species==spVecFollowing[i])])) Following_meanRc<-mean(SuperMetaDataFollowing$Rc[which(SuperMetaDataFollowing$species==spVecFollowing[i])]) Following_meanTA<-mean(SuperMetaDataFollowing$R_TA[which(SuperMetaDataFollowing$species==spVecFollowing[i])]) text(SuperMetaDataFollowing$R_TA[which(SuperMetaDataFollowing$species==spVecFollowing[i])],SuperMetaDataFollowing$Rc[which(SuperMetaDataFollowing$species==spVecFollowing[i])],round(SuperMetaDataFollowing$Qunat_CRW_VM[which(SuperMetaDataFollowing$species==spVecFollowing[i])]),col=colorVec[i]) points(SuperMetaDataFollowing$R_TA[which(SuperMetaDataFollowing$species==spVecFollowing[i])],SuperMetaDataFollowing$Rc[which(SuperMetaDataFollowing$species==spVecFollowing[i])],pch=17,col=colorVec[i],cex=0.5) SuperMetaDataFollowing$Delta_Rc[which(SuperMetaDataFollowing$species==spVecFollowing[i])] SP_DataFollowing$mean_R_TA[which(SP_DataFollowing$Species==spVecFollowing[i])]= Following_meanTA SP_DataFollowing$mean_Rc[which(SP_DataFollowing$Species==spVecFollowing[i])]= Following_meanRc SP_DataFollowing$Rc_expected_CRW[which(SP_DataFollowing$Species==spVecFollowing[i])]=approx(fit_21$x, fit_21$y, Following_meanTA)[[2]] } SP_DataFollowing$Delta_Rc=SP_DataFollowing$mean_Rc-SP_DataFollowing$Rc_expected_CRW ### points(fit_21_50$x,fit_21_mean$y,col="yellow",pch=19,cex=0.5) ########### #############Crosses plot plot(c(0,1),c(0,1),xlim=c(0.2,1),ylim=c(0.2,1),pch=20,cex=0.5,xlab="R_d_theta",ylab="R_theta") yy=colMeans(Rvec_traj_21) xx=colMeans(Rvec_ta_21) length(Q_21) fit_21_01<-smooth.spline(mid_TA_21,Q_21_01,df=16) #16 degrees of freedom for the quantiles smoothing fit_21_05<-smooth.spline(mid_TA_21,Q_21_05,df=16) fit_21_10<-smooth.spline(mid_TA_21,Q_21_1,df=16) fit_21_20<-smooth.spline(mid_TA_21,Q_21_2,df=16) fit_21_30<-smooth.spline(mid_TA_21,Q_21_3,df=16) fit_21_40<-smooth.spline(mid_TA_21,Q_21_4,df=16) fit_21_50<-smooth.spline(mid_TA_21,Q_21_5,df=16) fit_21_60<-smooth.spline(mid_TA_21,Q_21_6,df=16) fit_21_70<-smooth.spline(mid_TA_21,Q_21_7,df=16) fit_21_80<-smooth.spline(mid_TA_21,Q_21_8,df=16) fit_21_90<-smooth.spline(mid_TA_21,Q_21_9,df=16) fit_21_95<-smooth.spline(mid_TA_21,Q_21_95,df=16) fit_21_99<-smooth.spline(mid_TA_21,Q_21_99,df=16) lines(fit_21_05$x,fit_21_05$y, col="grey",lwd=3) lines(fit_21_10$x,fit_21_10$y, col="grey") lines(fit_21_20$x,fit_21_20$y, col="grey") lines(fit_21_30$x,fit_21_30$y, col="grey") lines(fit_21_40$x,fit_21_40$y, col="grey") lines(fit_21_50$x,fit_21_50$y, col="grey") lines(fit_21_60$x,fit_21_60$y, col="grey") lines(fit_21_70$x,fit_21_70$y, col="grey") lines(fit_21_80$x,fit_21_80$y, col="grey") lines(fit_21_90$x,fit_21_90$y, col="grey") lines(fit_21_95$x,fit_21_95$y, col="grey",lwd=3) fit_21_mean<-smooth.spline(mid_TA_21,colMeans(Rvec_traj_21),df=16) #16 degrees of freedom #delta Rc - ditance from CRW curve SP_DataFollowing<-matrix(0,nrow=length(spVecFollowing),ncol=5) SP_DataFollowing<-as.data.frame(SP_DataFollowing) names(SP_DataFollowing)<-c("Species","mean_R_TA","mean_Rc","Rc_expected_CRW","Delta_Rc") SP_DataFollowing$Species=spVecFollowing #plot the crosses for (i in 1:length(spVecFollowing)){ Following_TA_ci<- qnorm(0.975)*sd(SuperMetaDataFollowing$R_TA[which(SuperMetaDataFollowing$species==spVecFollowing[i])])/sqrt(length(SuperMetaDataFollowing$R_TA[which(SuperMetaDataFollowing$species==spVecFollowing[i])])) #compute CIs Following_Rc_ci<- qnorm(0.975)*sd(SuperMetaDataFollowing$Rc[which(SuperMetaDataFollowing$species==spVecFollowing[i])])/sqrt(length(SuperMetaDataFollowing$Rc[which(SuperMetaDataFollowing$species==spVecFollowing[i])])) #compute CIs Following_meanRc<-mean(SuperMetaDataFollowing$Rc[which(SuperMetaDataFollowing$species==spVecFollowing[i])]) Following_meanTA<-mean(SuperMetaDataFollowing$R_TA[which(SuperMetaDataFollowing$species==spVecFollowing[i])]) lines(c(Following_meanTA+Following_TA_ci,Following_meanTA-Following_TA_ci),c(Following_meanRc,Following_meanRc),pch=17,col=colorVec[i],lwd = 3) lines(c(Following_meanTA,Following_meanTA),c(Following_meanRc+Following_Rc_ci,Following_meanRc-Following_Rc_ci),pch=17,col=colorVec[i],lwd = 3) SuperMetaDataFollowing$Delta_Rc[which(SuperMetaDataFollowing$species==spVecFollowing[i])] SP_DataFollowing$mean_R_TA[which(SP_DataFollowing$Species==spVecFollowing[i])]= Following_meanTA SP_DataFollowing$mean_Rc[which(SP_DataFollowing$Species==spVecFollowing[i])]= Following_meanRc SP_DataFollowing$Rc_expected_CRW[which(SP_DataFollowing$Species==spVecFollowing[i])]=approx(fit_21$x, fit_21$y, Following_meanTA)[[2]] } SP_DataFollowing$Delta_Rc=SP_DataFollowing$mean_Rc-SP_DataFollowing$Rc_expected_CRW rankVec<-(0.3+round(rank(SP_DataFollowing$mean_Rc)/length(SP_DataFollowing$mean_Rc),2))/(max((0.3+round(rank(SP_DataFollowing$mean_Rc)/length(SP_DataFollowing$mean_Rc),2)))) #plot the dotted line connecting between species name and the center of the cross for (i in 1:length(spVecFollowing)){ xxx<-c(0.35,mean(SuperMetaDataFollowing$R_TA[which(SuperMetaDataFollowing$species==spVecFollowing[i])])) yyy<-c(rankVec[i],mean(SuperMetaDataFollowing$Rc[which(SuperMetaDataFollowing$species==spVecFollowing[i])])) lines(xxx,yyy,lty="dotted",col=colorVec[i]) name<-gsub("_"," ",spVecFollowing[i]) text(0.35, rankVec[i],name, col=colorVec[i],font=3,cex=1) } lines(fit_21_mean$x,fit_21_mean$y, col="blue",lwd=4)#plot the means line points(fit_21_50$x,fit_21_mean$y,col="yellow",pch=19,cex=0.5)#plot the points of mean CRW-VM theta ~ CRW-VM TA