### Determine first and last lat and long on migration ############################################################## library(rgdal) data3 <- subset(data2,abs(data2$long.m.spd) < 50) ## First and last fix per bird first.dts <- aggregate(data3$dt,by=list(data3$name),FUN='min') colnames(first.dts) <- c("name","start") attributes(first.dts$start)$tzone <- "UTC" last.dts <- aggregate(data3$dt,by=list(data3$name),FUN='max') colnames(last.dts) <- c("name","end") attributes(last.dts$end)$tzone <- "UTC" # Merge starts and ends to df startends <- merge(first.dts,last.dts) data3 <- merge(data3,startends,all.x=T) ## Determine end lat in degrees end.lats <- data3[which(data3$dt == data3$end),c("name","lat")] colnames(end.lats)[2] <- c("end.lat") data3 <- merge(data3,end.lats,all.x=T) ## Determine end long in degrees end.longs <- data3[which(data3$dt == data3$end),c("name","long")] colnames(end.longs)[2] <- c("end.long") data3 <- merge(data3,end.longs,all.x=T) ## Determine start long in degrees st.longs <- data3[which(data3$dt == data3$start),c("name","long")] colnames(st.longs)[2] <- c("st.long") data3 <- merge(data3,st.longs,all.x=T) #subset data to include only flight: speed threshold 5km per hour over sea and during daylight hours data3 <- data3[which(data3$spd >=5 &((data3$hr >= 7 & data3$hr < 21)| data3$landmask == 'sea')),] # set spd in scale m/s instead of km/h data3$spd.m <- data3$spd / 3.6 # remove outliers in terms of spd data3 <- subset(data3,data3$spd.m <=25) data3 <- data3[order(data3$migr,data3$dt),] #### Run Multiple Linear Regression Models #################################################### m1 <- lm(long.m.spd~ u ,data=data3) m2 <- lm(long.m.spd~ u + v,data=data3) m3 <- lm(long.m.spd~ u * v,data=data3) (m<-anova(m1,m2)) (m<-anova(m1,m3)) (m<-anova(m2,m3)) data3$res <- data3$long.m.spd - predict(m3) outliers <- which(abs(data3$res)>14) data3 <- data3[-outliers,] m1 <- lm(long.m.spd~ u ,data=data3) m2 <- lm(long.m.spd~ u + v,data=data3) m3 <- lm(long.m.spd~ u * v,data=data3) (m<-anova(m1,m2)) (m<-anova(m1,m3)) (m<-anova(m2,m3)) #### PLOT 3D RESPONSE SURFACE AND RESIDUAL VALUES ####################################################### library(rgl) myColorRamp <- function(colors, values) { v <- (values - min(values))/diff(range(values)) x <- colorRamp(colors)(v) rgb(x[,1], x[,2], x[,3], maxColorValue = 255) } cols <- myColorRamp(c("blue", "yellow","red"), data3$res) newdat <- expand.grid(u=seq(-20,20,by=.5),v=seq(-15,15,by=.5)) newdat$long.m.spd <- predict(m3,newdat) dev.new() with(data3,plot3d(u,v,long.m.spd, col=cols, size=.5,type="s", main="")) with(newdat,surface3d(unique(u),unique(v),long.m.spd, alpha=0.3,front="line", back="line")) rgl.snapshot("Fig2A.png", fmt = "png", top = TRUE ) ### PLOT FIGURE 2B #################################################### ## To do this replace U-wind component in script for figure 1 by model residuals ### MODEL DELTA-LONG AS A FUNCTION OF WIND EN ROUTE ####################################################### meanu <- aggregate(data3$u,by=list(data3$name),FUN='mean') colnames(meanu)[1:2] <- c("name","meanu") meanv <- aggregate(data3$v,by=list(data3$name),FUN='mean') colnames(meanv)[1:2] <- c("name","meanv") finlong <- unique(data3[,c("name","st.long","end.long","end.lat","cat")]) birds <- merge(meanu,finlong) birds <- merge(birds,meanv) # Subset individuals that survived until end of migration based on final lat survivors <- subset(birds,birds$end.lat.d < 20) m1 <- lm((end.long-st.long)~meanu,data=survivors) m2 <- lm((end.long-st.long)~meanu+meanv,data=survivors) m3 <- lm((end.long-st.long)~meanu*meanv,data=survivors) m4 <- lm((end.long-st.long)~meanu+meanv+cat,data=survivors) m5 <- lm((end.long-st.long)~meanu*meanv+cat,data=survivors) (m<-anova(m1,m2)) (m<-anova(m2,m3)) (m<-anova(m2,m4)) (m<-anova(m3,m5)) survivors$pred <- predict(m2) survivors$obs <- survivors$end.long - survivors$st.long #### PLOT FIGURE 3 ############################################### # determine which outleirs need to be labelled outliers <- c("Matti","Hans","Venus") # Generate plot ggplot() + geom_point(data=survivors,aes(x=obs,y=pred))+ theme_bw() + ylab('Predicted ?long') + xlab('Observed ?long') + geom_abline(intercept=0,slope=1)+ geom_text(data=survivors[which(survivors$name %in% outliers),],aes(obs,pred,label=name),hjust = 0, nudge_x = 0.25)+ theme(legend.position = 'bottom', legend.direction= 'horizontal', plot.background = element_rect(fill='white'), panel.grid.major = element_line(colour='grey50',linetype='dashed'), panel.grid.minor = element_line(colour='grey70',linetype='dashed'), legend.background = element_rect(fill='white'), panel.background = element_rect(fill='white'), strip.background = element_rect(fill='white'), legend.title = element_text(size=10,colour='black'), legend.text = element_text(size=10,colour='black'), strip.text = element_text(size=14,colour='black',face='bold'), axis.text = element_text(size=12,colour='black',face='italic'), axis.title = element_text(size=12,colour='black',face='bold')) ggsave('Fig3.tiff',width=4,height=4,dpi=300)