####################################################################################################################
###################models for Temporal changes to migratory fuel-loading in migratory birds across Europe###########
####################################################################################################################
#R version 4.2.0 (2022-04-22 ucrt)
#Platform: x86_64-w64-mingw32/x64 (64-bit)
#Running under: Windows 10 x64 (build 19045)

library(data.table) ##data.table_1.16.0
library(dplyr) ## dplyr_1.1.4
library(plyr) ## plyr_1.8.7
library(tidyr) ## tidyr_1.3.1
library(ggplot2) ##ggplot2_3.5.1 
library(glmmTMB) ## glmmTMB_1.1.3  ##Matrix_1.4-1   ##lme4_1.1-35.1
library(emmeans) ## emmeans_1.8.4-1

path<-"F:/weight_analysis/"


SMI_max_per_yr<-readRDS(file=paste0(path, "data/SMI_max_year_winglength3rdprimary.rds"))

##remove rustic bunting as not an afro-palearctic migrant
SMI_max_per_yr<-subset(SMI_max_per_yr, !Eng_Name=="Rustic Bunting")

##remove countries not on the western flyway
SMI_max_per_yr<-subset(SMI_max_per_yr, !Country%in% c("Kazakhstan", "Ukraine", "Cyprus"))

##scale northing
SMI_max_per_yr$Northing_10km_scal<-scale(SMI_max_per_yr$Northing_10km)

##remove line 14754 as this is an outlier I think Garden warbler has value of 50
SMI_max_per_yr1<-SMI_max_per_yr[-14754,]

########################################################
###run models max SMI against year and northing#########
#######################################################

#year*northing interaction spring, species with fewest data points removed to allow model convergence
mod<-glmmTMB(SMI~year*Northing_10km_scal+(1|km10)+(year|Eng_Name),data=SMI_max_per_yr1[SMI_max_per_yr1$Season=="spring"&!(SMI_max_per_yr1$Eng_Name %in% c("Melodius warbler")),])
#no year northing interaction spring, species with fewst data points removed to allow model convergence
mod_noint<-glmmTMB(SMI~year+Northing_10km_scal+(1|km10)+(year|Eng_Name),data=SMI_max_per_yr1[SMI_max_per_yr1$Season=="spring"&!(SMI_max_per_yr1$Eng_Name %in% c("Melodius warbler")),])
##year*northing interaction autumn
moda<-glmmTMB(SMI~year*Northing_10km_scal+(1|km10)+(year|Eng_Name), data=SMI_max_per_yr[SMI_max_per_yr$Season=="autumn",])
#no year northing interaction autumn
moda_noint<-glmmTMB(SMI~year+Northing_10km_scal+(1|km10)+(year|Eng_Name), data=SMI_max_per_yr[SMI_max_per_yr$Season=="autumn",])

summary(mod)
summary(mod_noint)
summary(moda)
summary(moda_noint)


###and get varibales needed to back transform northing
unscale_mean<-mean(SMI_max_per_yr$Northing_10km)
northing_unscle_mean_subtracted<-SMI_max_per_yr$Northing_10km-unscale_mean
unscale_std_dv<-sd(northing_unscle_mean_subtracted)


##unscale spring northing effect parameter estimate to days per 100km
##first get upper and lower CI
UCI<-   -0.400491+1.96* 0.079337
LCI<-   -0.400491-1.96* 0.079337
( -0.400491/unscale_std_dv)*100000
UCI/unscale_std_dv*100000
LCI/unscale_std_dv*100000

##unscale autumn northing effect parameter estimate to get days per 100km
UCI<-   -0.477601+1.96*0.035016
LCI<-   -0.477601-1.96*0.035016
( -0.477601/unscale_std_dv)*100000
UCI/unscale_std_dv*100000
LCI/unscale_std_dv*100000

####################################################################
#####models for week of max SMI against northing and year###########
####################################################################

###northing * year interaction
#need to remove 1 species with least data to converge for ineraction model
mod1_int<-glmmTMB(week~year*Northing_10km_scal+(1|km10)+(year|Eng_Name), data=SMI_max_per_yr[SMI_max_per_yr$Season=="spring"&!SMI_max_per_yr$Eng_Name%in%c("Melodius warbler"),])
###northing and year no interaction
mod1<-glmmTMB(week~year+Northing_10km_scal+(1|km10)+(year|Eng_Name), data=SMI_max_per_yr[SMI_max_per_yr$Season=="spring",])

##northing * year interaction
##had to remove one species to allow convergence
mod1a<-glmmTMB(week~year*Northing_10km_scal+(1|km10)+(year|Eng_Name), data=SMI_max_per_yr[SMI_max_per_yr$Season=="autumn"&!SMI_max_per_yr$Eng_Name %in% c("Collared Flycatcher"),])
## no northing  year interaction
mod1a_noint<-glmmTMB(week~year+Northing_10km_scal+(1|km10)+(year|Eng_Name), data=SMI_max_per_yr[SMI_max_per_yr$Season=="autumn"&!SMI_max_per_yr$Eng_Name %in% c("Collared Flycatcher"),])

summary(mod1)
summary(mod1_int)
summary(mod1a)
summary(mod1a_noint)


##unscale spring northing effect parameter estimate
##first get upper and lower CI
UCI<-  1.101809 +1.96*  0.111077
LCI<-  1.101809 -1.96*  0.111077
( 1.101809 /unscale_std_dv)*100000*7
UCI/unscale_std_dv*100000*7
LCI/unscale_std_dv*100000*7

###unscale autumn northign effect parameter estimate
UCI<-  -0.509145 +1.96*  0.101788
LCI<-  -0.509145 -1.96*  0.101788
( -0.509145 /unscale_std_dv)*100000*7
UCI/unscale_std_dv*100000*7
LCI/unscale_std_dv*100000*7


#################################################################
#######plot for spring max SMI against year and northing#########
#################################################################

##produce predicted plots of the effects
##create a new precition data frame
##set a species and site dummy variable this won't be used in the prediction as re.form=NA but is needed for the prediction to run
newdat<-data.frame(year=1980:2020, Northing_10km_scal=min(SMI_max_per_yr$Northing_10km_scal[SMI_max_per_yr$Season=="spring"]), km10="2940000_1590000", Eng_Name="Sedge Warbler")
newdat1<-data.frame(year=1980:2020, Northing_10km_scal=max(SMI_max_per_yr$Northing_10km_scal[SMI_max_per_yr$Season=="spring"]), km10="3660000_4100000", Eng_Name="Sedge Warbler")

###predict from model exclude random efefcts with re.form=NA
pred_yr<-predict(mod, newdata=newdat, re.form=NA, se=T)
pred_yr1<-predict(mod, newdata=newdat1, re.form=NA, se=T)

##make dataframe from prediction
yr_df<-data.frame(newdat, SMI=pred_yr$fit, LCI=pred_yr$fit-1.96*pred_yr$se.fit, UCI=pred_yr$fit+1.96*pred_yr$se.fit)
yr_df$latitude<-"south"
yr_df1<-data.frame(newdat1,SMI=pred_yr1$fit, LCI=pred_yr1$fit-1.96*pred_yr1$se.fit, UCI=pred_yr1$fit+1.96*pred_yr1$se.fit)
yr_df1$latitude<-"north"
yr_df<-rbind(yr_df, yr_df1)

##get summarised raw data for plot
##set median value to divide raw data into north and south
mid_north<-2845000 
lower<-SMI_max_per_yr[SMI_max_per_yr1$Season=="spring"&SMI_max_per_yr1$Northing_10km<=mid_north& !SMI_max_per_yr1$Eng_Name %in%  
                        c("Melodius warbler"), ]

upper<-SMI_max_per_yr[SMI_max_per_yr1$Season=="spring"&SMI_max_per_yr1$Northing_10km>=mid_north & !SMI_max_per_yr1$Eng_Name %in%  
                        c("Melodius warbler"), ]

##then get the species medians each year for upper(north) and lower(south)
##get medians
lower_median<-aggregate(data=lower, SMI~year+Eng_Name, function(x)quantile(x,0.5))
upper_median<-aggregate(data=upper, SMI~year+Eng_Name, function(x)quantile(x,0.5))
##get volume of data for each
lower_number<-aggregate(data=lower, SMI~year+Eng_Name, function(x)length(x))
upper_number<-aggregate(data=upper, SMI~year+Eng_Name, function(x)length(x))
##rename columns for clarity
names(lower_number)<-c("year", "Eng_Name", "N")
names(upper_number)<-c("year", "Eng_Name", "N")
##join tog
lower_median<-left_join(lower_median, lower_number)
upper_median<-left_join(upper_median, upper_number)

##plot
p<-ggplot()+geom_point(data= lower_median,aes(y=SMI, x=year), colour="#E69F00", alpha=lower_median$N/10, position=position_jitter(width=0.25))+
  geom_point(data= upper_median,aes(y=SMI, x=year), colour="#0072B2", alpha=upper_median$N/10,position=position_jitter(width=0.25))

p+geom_line(data=yr_df, aes(x=year, y=SMI, colour=latitude))+geom_line(data=yr_df, aes(x=year, y=LCI, colour=latitude), linetype=2)+
  geom_line(data=yr_df, aes(x=year, y=UCI, colour=latitude), linetype=2)+theme_classic()+ylab("Max SMI")+
  scale_colour_manual("Latitude", values=c("north"="#0072B2", "south"="#E69F00"))+ylim(c(6,40))


########################################################################
##############max SMI against year and northing plot autumn#############
#########################################################################

##produce predicted plots of the effects
##this is the same process as the plots above

##create prediction data.frame
newdat<-data.frame(year=1980:2020, Northing_10km_scal=min(SMI_max_per_yr$Northing_10km_scal[SMI_max_per_yr$Season=="autumn"]), km10="2940000_1590000", Eng_Name="Sedge Warbler")
newdat1<-data.frame(year=1980:2020, Northing_10km_scal=max(SMI_max_per_yr$Northing_10km_scal[SMI_max_per_yr$Season=="autumn"]), km10="3660000_4100000", Eng_Name="Sedge Warbler")

##predict for model
pred_yr<-predict(moda, newdata=newdat, re.form=NA, se=T)
pred_yr1<-predict(moda, newdata=newdat1, re.form=NA, se=T)

yr_df<-data.frame(newdat, SMI=pred_yr$fit, LCI=pred_yr$fit-1.96*pred_yr$se.fit, UCI=pred_yr$fit+1.96*pred_yr$se.fit)
yr_df$latitude<-"south"
yr_df1<-data.frame(newdat1, SMI=pred_yr1$fit, LCI=pred_yr1$fit-1.96*pred_yr1$se.fit, UCI=pred_yr1$fit+1.96*pred_yr1$se.fit)
yr_df1$latitude<-"north"
yr_df<-rbind(yr_df, yr_df1)

##divide raw data into north and south
##the mid point is different for autumn and spring as we have different sides included in each
mid_north<-3245000 
lower<-SMI_max_per_yr[SMI_max_per_yr$Season=="autumn"&SMI_max_per_yr$Northing_10km<=mid_north,]
upper<-SMI_max_per_yr[SMI_max_per_yr$Season=="autumn"&SMI_max_per_yr$Northing_10km>=mid_north,]

##then get the species medians each year for upper and lower from the raw data
lower_median<-aggregate(data=lower,SMI~year+Eng_Name, function(x)quantile(x,0.5))
upper_median<-aggregate(data=upper, SMI~year+Eng_Name, function(x)quantile(x,0.5))

lower_number<-aggregate(data=lower, SMI~year+Eng_Name, function(x)length(x))
upper_number<-aggregate(data=upper, SMI~year+Eng_Name, function(x)length(x))
names(lower_number)<-c("year", "Eng_Name", "N")
names(upper_number)<-c("year", "Eng_Name", "N")

lower_median<-left_join(lower_median, lower_number)
upper_median<-left_join(upper_median, upper_number)

p<-ggplot()+geom_point(data= lower_median,aes(y=SMI, x=year), colour="#E69F00", alpha=lower_median$N/100, position=position_jitter(width=0.25))+
  geom_point(data= upper_median,aes(y=SMI, x=year), colour="#0072B2", alpha=upper_median$N/100,position=position_jitter(width=0.25))

p+geom_line(data=yr_df, aes(x=year, y=SMI, colour=latitude))+geom_line(data=yr_df, aes(x=year, y=LCI, colour=latitude), linetype=2)+
  geom_line(data=yr_df, aes(x=year, y=UCI, colour=latitude), linetype=2)+theme_classic()+ylab("Max SMI")+
  scale_colour_manual("Latitude", values=c("north"="#0072B2", "south"="#E69F00"))+
  ylim(6,25)


######################################################################
########WEEK OF max SMI AUTUMN plot against year and northing########
####################################################################

##produce predicted plots of the effects

##create prediction data.frame
newdat<-data.frame(year=1980:2020, Northing_10km_scal=min(SMI_max_per_yr$Northing_10km_scal[SMI_max_per_yr$Season=="autumn"]), km10="2940000_1590000", Eng_Name="Sedge Warbler")
newdat1<-data.frame(year=1980:2020, Northing_10km_scal=max(SMI_max_per_yr$Northing_10km_scal[SMI_max_per_yr$Season=="autumn"]), km10="3660000_4100000", Eng_Name="Sedge Warbler")
##predict for model
pred_yr<-predict(mod1a, newdata=newdat, re.form=NA, se=T)
pred_yr1<-predict(mod1a, newdata=newdat1, re.form=NA, se=T)

yr_df<-data.frame(newdat, week=pred_yr$fit, LCI=pred_yr$fit-1.96*pred_yr$se.fit, UCI=pred_yr$fit+1.96*pred_yr$se.fit)
yr_df$latitude<-"south"
yr_df1<-data.frame(newdat1,week=pred_yr1$fit, LCI=pred_yr1$fit-1.96*pred_yr1$se.fit, UCI=pred_yr1$fit+1.96*pred_yr1$se.fit)
yr_df1$latitude<-"north"
yr_df<-rbind(yr_df, yr_df1)

##divide raw data into north and south
mid_north<-3245000
lower<-SMI_max_per_yr[SMI_max_per_yr$Season=="autumn"&SMI_max_per_yr$Northing_10km<=mid_north& !SMI_max_per_yr$Eng_Name %in%  
                        c("Collared Flycatcher"), ]
upper<-SMI_max_per_yr[SMI_max_per_yr$Season=="autumn"&SMI_max_per_yr$Northing_10km>=mid_north& !SMI_max_per_yr$Eng_Name %in%  
                        c("Collared Flycatcher"), ]

##then get the species medians each year for upper (north) and lower (south) from raw data
lower_median<-aggregate(data=lower, week~year+Eng_Name, function(x)quantile(x,0.5))
upper_median<-aggregate(data=upper, week~year+Eng_Name, function(x)quantile(x,0.5))

lower_number<-aggregate(data=lower, week~year+Eng_Name, function(x)length(x))
upper_number<-aggregate(data=upper, week~year+Eng_Name, function(x)length(x))
names(lower_number)<-c("year", "Eng_Name", "N")
names(upper_number)<-c("year", "Eng_Name", "N")

lower_median<-left_join(lower_median, lower_number)
upper_median<-left_join(upper_median, upper_number)


p<-ggplot()+geom_point(data= lower_median,aes(y=week, x=year), colour="#E69F00", alpha=lower_median$N/100, position=position_jitter(width=0.25))+
  geom_point(data= upper_median,aes(y=week, x=year), colour="#0072B2", alpha=upper_median$N/100,position=position_jitter(width=0.25))

p+geom_line(data=yr_df, aes(x=year, y=week, colour=latitude))+geom_line(data=yr_df, aes(x=year, y=LCI, colour=latitude), linetype=2)+
  geom_line(data=yr_df, aes(x=year, y=UCI, colour=latitude), linetype=2)+theme_classic()+ylab("Week of max SMI")+
  scale_colour_manual("Latitude", values=c("north"="#0072B2", "south"="#E69F00"))+
  ylim(31,40)


################################################################################################################################################
########################################LOOK AT CLIMATE EFFECTS###################################################################################
##################################################################################################################################################

##read in temperature anomly data
tmp_anomally_aut<-readRDS(file=paste0(path, "data/tmp_anomally_per_yr_3rd_primary.rds"))

##change to long format
tmp_anomally_aut1<-gather(data=tmp_anomally_aut, key="year", value="temp_anomally", -Easting_10km, -Northing_10km, -km10, -over_all_yrs_mean)
head(tmp_anomally_aut1)
tmp_anomally_aut1$year<-sub("Y", "", tmp_anomally_aut1$year)
tmp_anomally_aut1$year<-as.numeric(tmp_anomally_aut1$year)


##look at relationship between year and temp anomly
mod_temp<-glmmTMB(temp_anomally~year+scale(Northing_10km)+scale(Easting_10km), data=tmp_anomally_aut1)
summary(mod_temp) 
                     
                      
#subset max SMI data just to autumn
SMI_max_per_yr_aut<-subset(SMI_max_per_yr, Season=="autumn")

##join temp anomaly data to max SMI data
SMI_max_per_yr_aut<-left_join(SMI_max_per_yr_aut, tmp_anomally_aut1)


##and save
#saveRDS(SMI_max_per_yr_aut, file=paste0(path, "data/SMI_max_per_year_withmeantemp_andtempanomally_3rdprimary_westernpath.rds"))

##read in autumn max SMI data with temperature anomly data extracted
##code to extract temperature anomaly data available on request as it uses a rawer version of the ringing data which we do not have permission to make publically available
#SMI_max_per_yr_aut<-readRDS(file=paste0(path, "data/SMI_max_per_year_withmeantemp_andtempanomally_3rdprimary_westernpath.rds"))


##add factor for pos or negative temp anomaly as preliminary analysis shows a different relationship between max SMI and pos and neg temp anomlies
SMI_max_per_yr_aut$type_anomally[SMI_max_per_yr_aut$temp_anomally<0]<-"neg"
SMI_max_per_yr_aut$type_anomally[SMI_max_per_yr_aut$temp_anomally>0]<-"pos"

################################################################
#####model max SMI against temperature anomaly##################
################################################################
##with northing interaction
mod2<-glmmTMB(SMI~temp_anomally*type_anomally*Northing_10km_scal+(1|km10)+(1|Eng_Name), data=SMI_max_per_yr_aut)
summary(mod2)
##not northing interaction
mod2a<-glmmTMB(SMI~temp_anomally*type_anomally+Northing_10km_scal+(1|km10)+(1|Eng_Name), data=SMI_max_per_yr_aut)
summary(mod2a)


##create prediction data.frame for both positive and negative temperature anomalies for southern europe
temp_anomally_neg<-seq(from=min(SMI_max_per_yr_aut$temp_anomally, na.rm=T), 0, by=0.1)
temp_anomally_pos<-seq(0, to=max(SMI_max_per_yr_aut$temp_anomally, na.rm=T), by=0.1)
newdat_neg<-data.frame(temp_anomally=temp_anomally_neg, type_anomally="neg", Northing_10km_scal=min(SMI_max_per_yr_aut$Northing_10km_scal), Eng_Name="Sedge Warbler", km10="3620000_3290000")
newdat_pos<-data.frame(temp_anomally=temp_anomally_pos, type_anomally="pos", Northing_10km_scal=min(SMI_max_per_yr_aut$Northing_10km_scal), Eng_Name="Sedge Warbler", km10="3620000_3290000")
newdat_df<-rbind(newdat_neg, newdat_pos)

##predict for model
p<-predict(mod2, newdata=newdat_df, re.form=NA, se=T)
LCI<-p$fit-1.96*p$se.fit
UCI<-p$fit+1.96*p$se.fit
low_north<-data.frame(newdat_df, SMI=p$fit,LCI=LCI, UCI=UCI, latitude="south")

##create prediction data.frame for both positive and negative temperature anomalies for northern europe
newdat_neg1<-data.frame(temp_anomally=temp_anomally_neg, type_anomally="neg", Northing_10km_scal=max(SMI_max_per_yr_aut$Northing_10km_scal), Eng_Name="Sedge Warbler", km10="4720000_3690000")
newdat_pos1<-data.frame(temp_anomally=temp_anomally_pos, type_anomally="pos", Northing_10km_scal=max(SMI_max_per_yr_aut$Northing_10km_scal), Eng_Name="Sedge Warblert", km10="4720000_3690000")
newdat_df1<-rbind(newdat_neg1, newdat_pos1)

##predict for model
p<-predict(mod2, newdata=newdat_df1, re.form=NA, se=T)
LCI<-p$fit-1.96*p$se.fit
UCI<-p$fit+1.96*p$se.fit
high_north<-data.frame(newdat_df1, SMI=p$fit,LCI=LCI, UCI=UCI, latitude="north")
df<-rbind(low_north, high_north)

##split raw data into north and south to plot
mid_north<-3245000 
lower<-SMI_max_per_yr_aut[SMI_max_per_yr_aut$Northing_10km<=mid_north,]
upper<-SMI_max_per_yr_aut[SMI_max_per_yr_aut$Northing_10km>=mid_north,]

##split into bins of temp anomaly values
breaks_lower<-seq(from=-2.7786, to=2.6665, by=0.1)
names_lower<-seq(from=-2.7286,to= 2.6165, by=0.1)

temp_anomally_bins_lower<-cut(lower$temp_anomally, breaks=breaks_lower, labels=names_lower)

breaks_upper<-seq(from=-3.6174, to=2.6352, by=0.1)
names_upper<-seq(from= -3.5674,to=2.5852 ,by=0.1)

temp_anomally_bins_upper<-cut(upper$temp_anomally, breaks=breaks_upper, labels=names_upper)

lower$temp_anomally_bins_lower<-temp_anomally_bins_lower
upper$temp_anomally_bins_upper<-temp_anomally_bins_upper

##then get the species medians each temp_anomally for upper and lower
lower_median<-aggregate(data=lower,SMI~temp_anomally_bins_lower+Eng_Name, function(x)quantile(x,0.5))
upper_median<-aggregate(data=upper, SMI~temp_anomally_bins_upper+Eng_Name, function(x)quantile(x,0.5))

names(lower_median)[1]<-"temp_anomally"
names(upper_median)[1]<-"temp_anomally"

##get number of data points for each bin
lower_number<-aggregate(data=lower, SMI~temp_anomally_bins_lower+Eng_Name, function(x)length(x))
upper_number<-aggregate(data=upper, SMI~temp_anomally_bins_upper+Eng_Name, function(x)length(x))
names(lower_number)<-c("temp_anomally", "Eng_Name", "N")
names(upper_number)<-c("temp_anomally", "Eng_Name", "N")

lower_median<-left_join(lower_median, lower_number)
upper_median<-left_join(upper_median, upper_number)

##plot
q<-ggplot()+geom_point(data=lower_median,aes(y=SMI, x=as.numeric(as.character(temp_anomally))), colour="#E69F00", fill=NA, alpha=lower_median$N/50,position=position_jitter(width=0.25))+
  geom_point(data= upper_median,aes(y=SMI, x=as.numeric(as.character(temp_anomally))), colour="#0072B2", alpha=upper_median$N/50,position=position_jitter(width=0.25))+
  xlab("Temperature anomally")+theme_classic()

q+geom_line(data=df, aes(x=temp_anomally, y=SMI, colour=latitude))+ylab("Max SMI")+geom_line(data=df,aes(x=temp_anomally, y=LCI, colour=latitude), linetype=2)+
  geom_line(data=df,aes(x=temp_anomally, y=UCI, colour=latitude), linetype=2)+
  scale_colour_manual("Latitude", values=c("north"="#0072B2", "south"="#E69F00"))+
  xlab("Temperature anomaly ˚C")+
  theme_classic()+ylim(6, 24)

###################################################################
###week of max SMI against temp anomaly and northing models#########
###################################################################

##interaction with northing and temp anomaly
mod3<-glmmTMB(week~temp_anomally*Northing_10km_scal+ (1|km10)+(1|Eng_Name), data=SMI_max_per_yr_aut)
summary(mod3)
##no interaction between northign and temp anomaly
mod3a<-glmmTMB(week~temp_anomally+Northing_10km_scal+ (1|km10)+(1|Eng_Name), data=SMI_max_per_yr_aut)
summary(mod3a)

##create prediction data.frame for south
temp_anomally<-seq(from=min(SMI_max_per_yr_aut$temp_anomally), to=max(SMI_max_per_yr_aut$temp_anomally), by=0.1)
newdat<-data.frame(temp_anomally=temp_anomally, Northing_10km_scal=min(SMI_max_per_yr_aut$Northing_10km_scal), Eng_Name="Sedge Warbler", km10="3620000_3290000")

p<-predict(mod3, newdata=newdat, re.form=NA, se=T)
LCI<-p$fit-1.96*p$se.fit
UCI<-p$fit+1.96*p$se.fit
south<-data.frame(newdat, week=p$fit,LCI=LCI, UCI=UCI, latitude="south")

##create prediction data.frame for north
newdat1<-data.frame(temp_anomally=temp_anomally, Northing_10km_scal=max(SMI_max_per_yr_aut$Northing_10km_scal), Eng_Name="Sedge Warbler", km10="4720000_3690000")

p<-predict(mod3, newdata=newdat1, re.form=NA, se=T)
LCI<-p$fit-1.96*p$se.fit
UCI<-p$fit+1.96*p$se.fit
north<-data.frame(newdat1, week=p$fit,LCI=LCI, UCI=UCI, latitude="north")
df<-rbind(south, north)

##split raw data into north and south
mid_north<-3245000 
lower<-SMI_max_per_yr_aut[SMI_max_per_yr_aut$Northing_10km<=mid_north,]
upper<-SMI_max_per_yr_aut[SMI_max_per_yr_aut$Northing_10km>=mid_north,]

##split into temp anomly bins
breaks_lower<-seq(from=-2.7786, to=2.6665, by=0.1)
names_lower<-seq(from=-2.7286,to= 2.6165, by=0.1)

temp_anomally_bins_lower<-cut(lower$temp_anomally, breaks=breaks_lower, labels=names_lower)

breaks_upper<-seq(from=-3.6174, to=2.6352, by=0.1)
names_upper<-seq(from= -3.5674,to=2.5852 ,by=0.1)

temp_anomally_bins_upper<-cut(upper$temp_anomally, breaks=breaks_upper, labels=names_upper)

lower$temp_anomally_bins_lower<-temp_anomally_bins_lower
upper$temp_anomally_bins_upper<-temp_anomally_bins_upper

##then get the species medians each temp_anomaly bin for upper and lower
lower_median<-aggregate(data=lower,week~temp_anomally_bins_lower+Eng_Name, function(x)quantile(x,0.5))
upper_median<-aggregate(data=upper, week~temp_anomally_bins_upper+Eng_Name, function(x)quantile(x,0.5))

names(lower_median)[1]<-"temp_anomally"
names(upper_median)[1]<-"temp_anomally"

lower_number<-aggregate(data=lower, week~temp_anomally_bins_lower+Eng_Name, function(x)length(x))
upper_number<-aggregate(data=upper, week~temp_anomally_bins_upper+Eng_Name, function(x)length(x))
names(lower_number)<-c("temp_anomally", "Eng_Name", "N")
names(upper_number)<-c("temp_anomally", "Eng_Name", "N")

lower_median<-left_join(lower_median, lower_number)
upper_median<-left_join(upper_median, upper_number)

##plot
q<-ggplot()+geom_point(data=lower_median,aes(y=week, x=as.numeric(as.character(temp_anomally))), colour="#E69F00", fill=NA, alpha=lower_median$N/50, position=position_jitter(width=0.25))+
  geom_point(data= upper_median,aes(y=week, x=as.numeric(as.character(temp_anomally))), colour="#0072B2", alpha=upper_median$N/50,position=position_jitter(width=0.25))+
  xlab("Temperature anomally")+theme_classic()

q+geom_line(data=df, aes(x=temp_anomally, y=week, colour=latitude))+ylab("Week of Max SMI")+geom_line(data=df,aes(x=temp_anomally, y=LCI, colour=latitude), linetype=2)+
  geom_line(data=df,aes(x=temp_anomally, y=UCI, colour=latitude), linetype=2)+
  scale_colour_manual("Latitude", values=c("north"="#0072B2", "south"="#E69F00"))+
  xlab("Temperature anomaly ˚C")+
  theme_classic()+ylim(31, 40)

#######################################################################################
#######and looking at species traits with max SMI and temp anomaly#####################
#######################################################################################

species_traits<-read.csv(paste0(path, "data/species_traits_weight_paper1.csv"))

SMI_max_per_yr_aut1<-left_join(SMI_max_per_yr_aut, species_traits)

#####foraging with neg temp anomaly
mod_sp_traits_SMI2<-glmmTMB(SMI~temp_anomally*Northing_10km_scal*diet+(1|km10)+(1|Eng_Name), data=SMI_max_per_yr_aut1[SMI_max_per_yr_aut1$type_anomally=="neg",])
mod_sp_SMI<-glmmTMB(SMI~temp_anomally*Northing_10km_scal+diet*Northing_10km_scal+(1|km10)+(1|Eng_Name), data=SMI_max_per_yr_aut1[SMI_max_per_yr_aut1$type_anomally=="neg",])
anova(mod_sp_traits_SMI2, mod_sp_SMI)
##diet interaction sig improves model

##look at emtrends to better understand effect of diet
emtrends(mod_sp_traits_SMI2,  specs= pairwise ~diet |Northing_10km_scal, var = "temp_anomally", at=list(Northing_10km_scal=c(min(SMI_max_per_yr_aut1$Northing_10km_scal), max(SMI_max_per_yr_aut1$Northing_10km_scal))))

temp_anomally<-seq(from=min(SMI_max_per_yr_aut1$temp_anomally), to=0, by=0.1)

newdat_foraging_south<-data.frame(diet=rep(c("InsectivoreOnly", "Fruitivore", "AerialInsectivore" ), each=37), temp_anomally=rep(temp_anomally,3),
                                  Northing_10km_scal=min(SMI_max_per_yr_aut1$Northing_10km_scal[SMI_max_per_yr_aut1$temp_anomally<0]), km10="3620000_3290000", Eng_Name="Sedge Warbler")
pred_south<-predict(mod_sp_traits_SMI2, newdata=newdat_foraging_south, se=T, re.form=NA)
pred_south<-data.frame(pred=pred_south$fit, LCI=pred_south$fit-1.96*pred_south$se.fit, UCI=pred_south$fit+1.96*pred_south$se.fit,newdat_foraging_south)
pred_south$Northing_10km_scal<-"south"

newdat_foraging_north<-data.frame(diet=rep(c("InsectivoreOnly", "Fruitivore", "AerialInsectivore" ), each=37), temp_anomally=rep(temp_anomally,3),
                                  Northing_10km_scal=max(SMI_max_per_yr_aut$Northing_10km_scal[SMI_max_per_yr_aut1$temp_anomally<0]), km10="4720000_3690000", Eng_Name="Sedge Warbler")
pred_north<-predict(mod_sp_traits_SMI2, newdata=newdat_foraging_north, se=T, re.form=NA)
pred_north<-data.frame(pred=pred_north$fit, LCI=pred_north$fit-1.96*pred_north$se.fit, UCI=pred_north$fit+1.96*pred_north$se.fit,newdat_foraging_north)
pred_north$Northing_10km_scal<-"north"
pred_foraging<-rbind(pred_south, pred_north)

ggplot(data=pred_foraging)+geom_line(aes(x=temp_anomally, y=pred, colour=diet))+
  geom_ribbon(aes(x = temp_anomally, ymin = LCI, ymax = UCI, fill = diet), alpha = 0.2)+
  geom_line(aes(x=temp_anomally, y=UCI, colour=diet), linetype=2)+
  geom_line(aes(x=temp_anomally, y=LCI, colour=diet), linetype=2)+
  facet_wrap(~Northing_10km_scal)+
  theme_bw()+  scale_colour_manual("Diet", values=c("AerialInsectivore"="#E69F00", "Fruitivore"="#56B4E9", "InsectivoreOnly"="#009E73"))+
  scale_fill_manual("Diet", values=c("AerialInsectivore"="#E69F00", "Fruitivore"="#56B4E9", "InsectivoreOnly"="#009E73"))+
  xlab("Temperature anomaly ˚C")+ylab("max SMI")

##diet with pos temp anomaly
mod_sp_traits_SMI2a<-glmmTMB(SMI~temp_anomally*Northing_10km_scal*diet+(1|km10)+(1|Eng_Name), data=SMI_max_per_yr_aut1[SMI_max_per_yr_aut1$type_anomally=="pos",])
mod_sp_SMI<-glmmTMB(SMI~temp_anomally*Northing_10km_scal+diet*Northing_10km_scal+(1|km10)+(1|Eng_Name), data=SMI_max_per_yr_aut1[SMI_max_per_yr_aut1$type_anomally=="pos",])
anova(mod_sp_traits_SMI2a, mod_sp_SMI)
##diet interaction sig improves model
emtrends(mod_sp_traits_SMI2a, specs= pairwise ~ diet |Northing_10km_scal, var = "temp_anomally", at=list(Northing_10km_scal=c(min(SMI_max_per_yr_aut1$Northing_10km_scal), max(SMI_max_per_yr_aut1$Northing_10km_scal))))

##plot effect
temp_anomally<-seq(from=0,to=max(SMI_max_per_yr_aut1$temp_anomally), by=0.1)

newdat_foraging_south<-data.frame(diet=rep(c( "InsectivoreOnly", "Fruitivore", "AerialInsectivore" ), each=27), temp_anomally=rep(temp_anomally,3),
                                  Northing_10km_scal=min(SMI_max_per_yr_aut1$Northing_10km_scal), km10="3620000_3290000", Eng_Name="Sedge Warbler")
pred_south<-predict(mod_sp_traits_SMI2a, newdata=newdat_foraging_south, se=T, re.form=NA)
pred_south<-data.frame(pred=pred_south$fit, LCI=pred_south$fit-1.96*pred_south$se.fit, UCI=pred_south$fit+1.96*pred_south$se.fit,newdat_foraging_south)
pred_south$Northing_10km_scal<-"south"

newdat_foraging_north<-data.frame(diet=rep(c("InsectivoreOnly", "Fruitivore", "AerialInsectivore" ), each=27), temp_anomally=rep(temp_anomally,3),
                                  Northing_10km_scal=max(SMI_max_per_yr_aut$Northing_10km_scal), km10="4720000_3690000", Eng_Name="Sedge Warbler")
pred_north<-predict(mod_sp_traits_SMI2a, newdata=newdat_foraging_north, se=T, re.form=NA)
pred_north<-data.frame(pred=pred_north$fit, LCI=pred_north$fit-1.96*pred_north$se.fit, UCI=pred_north$fit+1.96*pred_north$se.fit,newdat_foraging_north)
pred_north$Northing_10km_scal<-"north"
pred_foraging<-rbind(pred_south, pred_north)

ggplot(data=pred_foraging)+geom_line(aes(x=temp_anomally, y=pred, colour=diet))+
  geom_ribbon(aes(x = temp_anomally, ymin = LCI, ymax = UCI, fill = diet), alpha = 0.2)+
  geom_line(aes(x=temp_anomally, y=UCI, colour=diet), linetype=2)+
  geom_line(aes(x=temp_anomally, y=LCI, colour=diet), linetype=2)+
  facet_wrap(~Northing_10km_scal)+
  theme_bw()+scale_colour_manual("Diet", values=c("AerialInsectivore"="#E69F00", "Fruitivore"="#56B4E9", "InsectivoreOnly"="#009E73"))+
  scale_fill_manual("Diet", values=c("AerialInsectivore"="#E69F00", "Fruitivore"="#56B4E9", "InsectivoreOnly"="#009E73"))+
  xlab("Temperature anomaly ˚C")+ylab("max SMI")

##breeding cycle with neg temp anomaly
mod_sp_traits_SMI4<-glmmTMB(SMI~temp_anomally*Northing_10km_scal*breeding_cycle+(1|km10)+(1|Eng_Name), data=SMI_max_per_yr_aut1[SMI_max_per_yr_aut1$type_anomally=="neg",])
mod_sp_SMI<-glmmTMB(SMI~temp_anomally*Northing_10km_scal+breeding_cycle*Northing_10km_scal+(1|km10)+(1|Eng_Name), data=SMI_max_per_yr_aut1[SMI_max_per_yr_aut1$type_anomally=="neg",])
anova(mod_sp_traits_SMI4, mod_sp_SMI)

##breeding cycle with pos temp anomaly
mod_sp_traits_SMI4a<-glmmTMB(SMI~temp_anomally*Northing_10km_scal*breeding_cycle+(1|km10)+(1|Eng_Name), data=SMI_max_per_yr_aut1[SMI_max_per_yr_aut1$type_anomally=="pos",])
mod_sp_SMI<-glmmTMB(SMI~temp_anomally*Northing_10km_scal+breeding_cycle*Northing_10km_scal+(1|km10)+(1|Eng_Name), data=SMI_max_per_yr_aut1[SMI_max_per_yr_aut1$type_anomally=="pos",])
anova(mod_sp_traits_SMI4a, mod_sp_SMI)

##number of broods with neg temp anomaly
mod_sp_traits_SMI4a<-glmmTMB(SMI~temp_anomally*Northing_10km_scal*brood_num+(1|km10)+(1|Eng_Name), data=SMI_max_per_yr_aut1[SMI_max_per_yr_aut1$type_anomally=="neg",])
mod_sp_SMI<-glmmTMB(SMI~temp_anomally*Northing_10km_scal+brood_num*Northing_10km_scal+(1|km10)+(1|Eng_Name), data=SMI_max_per_yr_aut1[SMI_max_per_yr_aut1$type_anomally=="neg",])
anova(mod_sp_traits_SMI4a, mod_sp_SMI)

##number of broods with pos temp anomaly
mod_sp_traits_SMI4a<-glmmTMB(SMI~temp_anomally*Northing_10km_scal*brood_num+(1|km10)+(1|Eng_Name), data=SMI_max_per_yr_aut1[SMI_max_per_yr_aut1$type_anomally=="pos",])
mod_sp_SMI<-glmmTMB(SMI~temp_anomally*Northing_10km_scal+brood_num*Northing_10km_scal+(1|km10)+(1|Eng_Name), data=SMI_max_per_yr_aut1[SMI_max_per_yr_aut1$type_anomally=="pos",])
anova(mod_sp_traits_SMI4a, mod_sp_SMI)


###############################################################################################
#######and looking at species traits with week of max SMI and temp anomaly#####################
###############################################################################################

##diet with temp anomaly
mod_sp_traits_week2<-glmmTMB(week~temp_anomally*Northing_10km_scal*diet+(1|km10)+(1|Eng_Name), data=SMI_max_per_yr_aut1)
mod_sp_week<-glmmTMB(week~temp_anomally*Northing_10km_scal+diet*Northing_10km_scal+(1|km10)+(1|Eng_Name), data=SMI_max_per_yr_aut1)
anova(mod_sp_traits_week2, mod_sp_week)
##diet interaction sig improves model

##emtends to understand pattern better
emtrends(mod_sp_traits_week2, spec=pairwise~ diet |Northing_10km_scal, var = "temp_anomally", at=list(Northing_10km_scal=c(min(SMI_max_per_yr_aut1$Northing_10km_scal), max(SMI_max_per_yr_aut1$Northing_10km_scal))))

##plot graph
temp_anomally<-seq(from=min(SMI_max_per_yr_aut1$temp_anomally), to=max(SMI_max_per_yr_aut1$temp_anomally), by=0.1)

newdat_foraging_south<-data.frame(diet=rep(c("AerialInsectivore", "Fruitivore", "InsectivoreOnly"), each=63), temp_anomally=rep(temp_anomally,3),
                                  Northing_10km_scal=min(SMI_max_per_yr_aut1$Northing_10km_scal), km10="3620000_3290000", Eng_Name="Sedge Warbler")
pred_south<-predict(mod_sp_traits_week2, newdata=newdat_foraging_south, se=T, re.form=NA)
pred_south<-data.frame(pred=pred_south$fit, LCI=pred_south$fit-1.96*pred_south$se.fit, UCI=pred_south$fit+1.96*pred_south$se.fit,newdat_foraging_south)
pred_south$Northing_10km_scal<-"south"

newdat_foraging_north<-data.frame(diet=rep(c("AerialInsectivore", "Fruitivore", "InsectivoreOnly"), each=63), temp_anomally=rep(temp_anomally,3),
                                  Northing_10km_scal=max(SMI_max_per_yr_aut$Northing_10km_scal), km10="4720000_3690000", Eng_Name="Sedge Warbler")
pred_north<-predict(mod_sp_traits_week2, newdata=newdat_foraging_north, se=T, re.form=NA)
pred_north<-data.frame(pred=pred_north$fit, LCI=pred_north$fit-1.96*pred_north$se.fit, UCI=pred_north$fit+1.96*pred_north$se.fit,newdat_foraging_north)
pred_north$Northing_10km_scal<-"north"
pred_foraging<-rbind(pred_south, pred_north)

ggplot(data=pred_foraging)+geom_line(aes(x=temp_anomally, y=pred, colour=diet))+
  geom_ribbon(aes(x = temp_anomally, ymin = LCI, ymax = UCI, fill = diet), alpha = 0.2)+
  geom_line(aes(x=temp_anomally, y=UCI, colour=diet), linetype=2)+
  geom_line(aes(x=temp_anomally, y=LCI, colour=diet), linetype=2)+
  facet_wrap(~Northing_10km_scal)+
  theme_bw()+scale_colour_manual("Diet", values=c("AerialInsectivore"="#E69F00", "Fruitivore"="#56B4E9", "InsectivoreOnly"="#009E73"))+
  scale_fill_manual("Diet", values=c("AerialInsectivore"="#E69F00", "Fruitivore"="#56B4E9", "InsectivoreOnly"="#009E73"))+
  xlab("Temperature anomaly ˚C")+ylab("Week of max SMI")


##number of broods with temp anomaly
mod_sp_traits_week5<-glmmTMB(week~temp_anomally*Northing_10km_scal*brood_num+(1|km10)+(1|Eng_Name), data=SMI_max_per_yr_aut1)
mod_sp_week<-glmmTMB(week~temp_anomally*Northing_10km_scal+brood_num*Northing_10km_scal+(1|km10)+(1|Eng_Name), data=SMI_max_per_yr_aut1)
anova(mod_sp_traits_week5, mod_sp_week)

###breeding cycle with temp anomaly
mod_sp_traits_week5<-glmmTMB(week~temp_anomally*Northing_10km_scal*breeding_cycle+(1|km10)+(1|Eng_Name), data=SMI_max_per_yr_aut1)
mod_sp_week<-glmmTMB(week~temp_anomally*Northing_10km_scal+breeding_cycle*Northing_10km_scal+(1|km10)+(1|Eng_Name), data=SMI_max_per_yr_aut1)
anova(mod_sp_traits_week5, mod_sp_week)

##emtrends to understand pattern better
emtrends(mod_sp_traits_week5, specs = pairwise ~ breeding_cycle | Northing_10km_scal, 
         var = "temp_anomally",
         at = list(breeding_cycle = c(22,32,42),
                   Northing_10km_scal = c(min(SMI_max_per_yr_aut1$Northing_10km_scal), max(SMI_max_per_yr_aut1$Northing_10km_scal))))



##plot graph
temp_anomally<-seq(from=min(SMI_max_per_yr_aut1$temp_anomally), to=max(SMI_max_per_yr_aut1$temp_anomally), by=0.1)

newdat_breeding_south<-data.frame(breeding_cycle=rep(c(22,32,42), each=63), temp_anomally=rep(temp_anomally,3),
                                  Northing_10km_scal=min(SMI_max_per_yr_aut1$Northing_10km_scal), km10="3620000_3290000", Eng_Name="Sedge Warbler")
pred_south<-predict(mod_sp_traits_week5, newdata=newdat_breeding_south, se=T, re.form=NA)
pred_south<-data.frame(pred=pred_south$fit, LCI=pred_south$fit-1.96*pred_south$se.fit, UCI=pred_south$fit+1.96*pred_south$se.fit,newdat_breeding_south)
pred_south$Northing_10km_scal<-"south"

newdat_breeding_north<-data.frame(breeding_cycle=rep(c(22,32,42), each=63), temp_anomally=rep(temp_anomally,3),
                                  Northing_10km_scal=max(SMI_max_per_yr_aut$Northing_10km_scal), km10="4720000_3690000", Eng_Name="Sedge Warbler")
pred_north<-predict(mod_sp_traits_week5, newdata=newdat_breeding_north, se=T, re.form=NA)
pred_north<-data.frame(pred=pred_north$fit, LCI=pred_north$fit-1.96*pred_north$se.fit, UCI=pred_north$fit+1.96*pred_north$se.fit,newdat_breeding_north)
pred_north$Northing_10km_scal<-"north"
pred_breeding<-rbind(pred_south, pred_north)
pred_breeding$breeding_cycle<-factor(as.character(pred_breeding$breeding_cycle))

ggplot(data = pred_breeding) +
  geom_line(aes(x = temp_anomally, y = pred, colour = breeding_cycle)) +
  geom_ribbon(aes(x = temp_anomally, ymin = LCI, ymax = UCI, fill = breeding_cycle), alpha = 0.2) +
  geom_line(aes(x = temp_anomally, y = UCI, colour = breeding_cycle), linetype = "dashed") +
  geom_line(aes(x = temp_anomally, y = LCI, colour = breeding_cycle), linetype = "dashed") +
  facet_wrap(~Northing_10km_scal) +
  theme_bw() +
  scale_colour_manual(
    "Breeding Cycle (days)", 
    values = c("22" = "#E69F00",  
               "32" = "#56B4E9",  
               "42" = "#009E73")) +
  scale_fill_manual(
    "Breeding Cycle (days)", 
    values = c("22" = "#E69F00",  
               "32" = "#56B4E9",  
               "42" = "#009E73"))+
  xlab("Temperature anomaly ˚C") +
  ylab("Week of max SMI")
