### The R script for Väli et al. 'Genetic determination of migration strategy in large soaring birds: evidence from hybrid eagles'
### Completed on July 5th, 2018

### The source data files are available as Excel sheets 
### in the Electronic Supplementary Material and must be converted to csv-files.

#Data for the analysis of timing the onset of migration
data_timing<-read.csv("Datatable_timing.csv", sep = ";")

data_timing$fNo.of.bird <- as.factor(data_timing$No.of.bird)
data_timing$fNo..of.breeding.territory <- as.factor(data_timing$No..of.breeding.territory)
data_timing$Species_1<-factor(data_timing$Species, 
                              levels = c("C. clanga", "Hybrids", "C. pomarina"))
summary(data_timing$Species)
table (data_timing$Tracking.year, data_timing$Species)

#Data for the analysis of location of wintering ranges (original and ordered groups)
data_latlong<-read.csv("Datatable_latlong.csv", sep = ";")
data_latlong$Species_1<-factor(data_latlong$Species, 
                                levels = c("C. clanga", "Hybrids", "C. pomarina"))
summary(data_latlong$Species)
table (data_latlong$Tracking.year, data_latlong$Species)

#Data for the analysis of the wintering range
data_range<-read.csv("Datatable_range.csv", sep = ";")
data_range$Species_1<-factor(data_range$Species, 
                                levels = c("C. clanga", "Hybrids", "C. pomarina"))

summary(data_range$Species)
table (data_range$Tracking.year, data_range$Species)
library(nlme)
library(multcomp)


###Onset of migration
fligner.test(Julian.start.autumn ~ Species, data=data_timing)
fligner.test(Julian.start.autumn ~ Species, data=data_timing[data_timing$Species=="C. clanga" | 
                                                               data_timing$Species=="C. pomarina",])
fligner.test(Julian.start.autumn ~ Species, data=data_timing[data_timing$Species=="C. clanga" | 
                                                               data_timing$Species=="Hybrids",])
fligner.test (Julian.start.autumn ~ Species, data=data_timing[data_timing$Species=="Hybrids" | 
                                                                data_timing$Species=="C. pomarina",])

lme_onset <- lme(Julian.start.autumn ~ Species,
           data=data_timing, random = list (~ 1|fNo.of.bird, ~ 1|fNo..of.breeding.territory, ~ 1|Tracking.year),
           weights=varIdent(form=~1|Species), method = "ML")
summary (lme_onset)
lme_onset_0<-lme(Julian.start.autumn ~ 1, 
              data=data_timing, random = list (~ 1|fNo.of.bird, ~ 1|fNo..of.breeding.territory, ~ 1|Tracking.year),
              weights=varIdent(form=~1|Species), method = "ML")
summary(lme_onset_0)
anova(lme_onset, lme_onset_0)
summary(glht(lme_onset, mcp(Species ="Tukey")))

plot (Julian.start.autumn ~ Species_1, data=data_timing,
      ylab="Julian date", xlab=" ", 
      font.lab=2, las = 1, lab=c(10,5,1), ylim = c(240,300), 
      outline = FALSE)
#GLMM Tukey results
text(1.5, 295, "***", cex=1)
text(2, 300, "***", cex=1)
text(2.5, 295, "ns", cex=0.75)
#FK test
text(1.5, 247, "ns", cex=0.75)
text(2, 242, "***", cex=1)
text(2.5, 247, "**", cex=1)
#Sample sizes
text(1, 290, "46", cex=1)
text(2, 279, "23", cex=1)
text(3, 271, "49", cex=1)

###Wintering latitude
fligner.test (Winter.latitude..wgs84. ~ Species, data=data_latlong)
fligner.test (Winter.latitude..wgs84. ~ Species, data=data_latlong[data_latlong$Species=="C. clanga" | 
                                                                     data_latlong$Species=="C. pomarina",])
fligner.test (Winter.latitude..wgs84. ~ Species, data=data_latlong[data_latlong$Species=="C. clanga" | 
                                                                     data_latlong$Species=="Hybrids",])
fligner.test (Winter.latitude..wgs84. ~ Species, data=data_latlong[data_latlong$Species=="Hybrids" | 
                                                                     data_latlong$Species=="C. pomarina",])

lme_lat<-lme(Winter.latitude..wgs84. ~ Species,
             random = list (~ 1|No.of.bird, ~ 1|No..of.breeding.territory), 
             weights=varIdent(form=~1|Species),
             method = "ML",
             data=data_latlong)
summary(lme_lat)
lme_lat_0<-lme(Winter.latitude..wgs84. ~ 1,
               random = list (~ 1|No.of.bird, ~ 1|No..of.breeding.territory), 
               weights=varIdent(form=~1|Species), 
               method = "ML",
               data=data_latlong)
anova(lme_lat, lme_lat_0)
summary(glht(lme_lat, mcp(Species ="Tukey")))

summary(data_latlong$Species)

plot (Winter.latitude..wgs84. ~ Species_1, data=data_latlong, 
      ylab="Latitude (degrees)", xlab=" ", 
      font.lab=2, las = 1, lab=c(10,5,1), ylim= c(-40,70), outline = FALSE)
#GLMM Tukey results
text(1.5, 62, "*", cex=1)
text(2, 70, "***", cex=1)
text(2.5, 62, "***", cex=1)
#F-K test results
text(1.5, -30, "ns", cex=0.75)
text(2, -38, "*", cex=1)
text(2.5, -30, "*", cex=1)
#Sample sizes
text(1, 49, "28", cex=1)
text(2, 48, "17", cex=1)
text(3, -6, "23", cex=1)
      
###Longitude
fligner.test (Winter.longitude..wgs84. ~ Species, data=data_latlong)
fligner.test (Winter.longitude..wgs84. ~ Species, data=data_latlong[data_latlong$Species=="C. clanga" | 
                                                                      data_latlong$Species=="C. pomarina",])
fligner.test (Winter.longitude..wgs84. ~ Species, data=data_latlong[data_latlong$Species=="C. clanga" | 
                                                                      data_latlong$Species=="Hybrids",])
fligner.test (Winter.longitude..wgs84. ~ Species, data=data_latlong[data_latlong$Species=="Hybrids" | 
                                                                      data_latlong$Species=="C. pomarina",])

lme_long<-lme(Winter.longitude..wgs84. ~ Species,
              random = list (~ 1|No.of.bird, ~ 1|No..of.breeding.territory), 
              weights=varIdent(form=~1|Species),
              method = "ML",
              data=data_latlong)
summary(lme_long)
summary(glht(lme_long, mcp(Species ="Tukey")))
lme_long_0<-lme(Winter.longitude..wgs84. ~ 1,
                random = list (~ 1|No.of.bird, ~ 1|No..of.breeding.territory),  
                weights=varIdent(form=~1|Species),
                method = "ML",
                data=data_latlong)
anova(lme_long, lme_long_0)

plot (Winter.longitude..wgs84. ~ Species_1, data=data_latlong,  
      ylab="Longitude (degrees)", xlab=" ", 
      font.lab=2, las = 1, lab=c(10,5,1), ylim= c(10,45), outline = FALSE)

#GLMM Tukey results
text(1.5, 42, "ns", cex=0.75)
text(2, 45, "ns", cex=0.75)
text(2.5, 42, "ns", cex=0.75)
#F-K test results
text(1.5, 13, "ns", cex=0.75)
text(2, 11, "***", cex=1)
text(2.5, 13, "**", cex=1)
#Sample sizes
text(1, 37.2, "28", cex=1)
text(2, 39, "17", cex=1)
text(3, 31.3, "23", cex=1)

###Wintering range size
fligner.test (Winter.range..km2. ~ Species, data=data_range)
fligner.test (Winter.range..km2. ~ Species, data=data_range[data_range$Species=="C. clanga" | 
                                                                    data_range$Species=="C. pomarina",])
fligner.test (Winter.range..km2. ~ Species, data=data_range[data_range$Species=="C. clanga" | 
                                                                    data_range$Species=="Hybrids",])
fligner.test (Winter.range..km2. ~ Species, data=data_range[data_range$Species=="Hybrids" | 
                                                                    data_range$Species=="C. pomarina",])

lme_winterarea<-lme(Winter.range..km2. ~ Species,
                    random = list (~ 1|No.of.bird, ~ 1|No..of.breeding.territory),  
                    weights=varIdent(form=~1|Species),
                    method = "ML",
                    data=data_range)
summary(lme_winterarea)
summary(glht(lme_winterarea, mcp(Species ="Tukey")))
lme_winterarea_0<-lme(Winter.range..km2. ~ 1,
                      random = list (~ 1|No.of.bird, ~ 1|No..of.breeding.territory),  
                      weights=varIdent(form=~1|Species),
                      method = "ML",
                      data=data_range)
anova(lme_winterarea, lme_winterarea_0)

plot (Winter.range..km2. ~ Species_1, data=data_range,   
      ylab="Wintering range (km2)", xlab="",
      font.lab=2, las = 1, lab=c(10,5,1), ylim= c(-3000,15000), outline = FALSE)


#GLMM Tukey results
text(1.5, 13500, "ns", cex=0.75)
text(2, 15000, "***", cex=1)
text(2.5, 13500, "**", cex=1)

#F-K test results
text(1.5, -1500, "**", cex=0.75)
text(2, -2500, "***", cex=1)
text(2.5, -1500, "*", cex=1)
#Sample sizes
text(1, 700, "28", cex=1)
text(2, 3300, "17", cex=1)
text(3, 11100, "23", cex=1)