############################################################################################
#                                                                                          #
#                                          PACKAGES                                        #
#                                                                                          #
############################################################################################

### packages used in analyses

library(ggplot2)
library(ggpubr)
library(rcompanion)
library(MuMIn)
library(car)
library(tidyverse)
library(broom)
library(betareg)
library(rcompanion)

############################################################################################
#                                                                                          #
#                                 LANDSCAPE, NESTS AND BROODS                              #
#                                                                                          #
############################################################################################

### plotting nest locations

nests<-subset(read.csv("XYZ.csv"),stage=="nest")

str(nests)
levels(as.factor(nests$sp))

### for 2018 and 2019

ggplot(nests,aes(X,Y,fill=sp,label=year))+
  geom_hline(yintercept=8440000)+
  geom_hline(yintercept=8452000)+
  geom_vline(xintercept=480000)+
  geom_vline(xintercept=495000)+
  theme_classic()+
  scale_y_continuous(breaks=seq(8440000,8452000,by=1000))+
  scale_x_continuous(breaks=seq(480000,495000,by=1000))+
  geom_point(size=6,shape=21)+
  geom_text(aes(label=year),size=3,color="white")

### for 2018

ggplot(subset(nests,year==18),aes(X,Y,fill=sp,label=year))+
  geom_hline(yintercept=8440000)+
  geom_hline(yintercept=8452000)+
  geom_vline(xintercept=480000)+
  geom_vline(xintercept=495000)+
  theme_classic()+
  scale_y_continuous(breaks=seq(8440000,8452000,by=1000))+
  scale_x_continuous(breaks=seq(480000,495000,by=1000))+
  geom_point(size=6,shape=21)+
  geom_text(aes(label=year),size=3,color="white")

### for 2019

ggplot(subset(nests,year==19),aes(X,Y,fill=sp,label=year))+
  geom_hline(yintercept=8440000)+
  geom_hline(yintercept=8452000)+
  geom_vline(xintercept=480000)+
  geom_vline(xintercept=495000)+
  theme_classic()+
  scale_y_continuous(breaks=seq(8440000,8452000,by=1000))+
  scale_x_continuous(breaks=seq(480000,495000,by=1000))+
  geom_point(size=6,shape=21)+
  geom_text(aes(label=year),size=3,color="white")

### plotting proportion of elevations in the study area

landscape<-read.csv("landscape.csv")

ggplot(landscape,aes(x=elevation,y=part))+
  theme_classic()+
  geom_col()+
  coord_flip()+
  scale_x_continuous(limits=c(0,240),breaks=seq(0,220,by=20))

### plotting nests and broods vs elevation

LS<-subset(read.csv("XYZ.csv"),sp=="LS")

ggplot(LS,aes(y=Z,x=stage,shape=stage,fill=as.factor(year)))+
  geom_boxplot(outlier.shape=NA)+
  geom_point(size=5,position=position_jitterdodge())+
  scale_shape_manual(values=c(24,21))+
  scale_y_continuous(limits=c(0,240),breaks=seq(0,220,by=20))+
  theme_classic()

RK<-subset(read.csv("XYZ.csv"),sp=="RK")

ggplot(RK,aes(y=Z,x=stage,shape=stage,fill=as.factor(year)))+
  geom_boxplot(outlier.shape=NA)+
  geom_point(size=5,position=position_jitterdodge())+
  scale_shape_manual(values=c(24,21))+
  scale_y_continuous(limits=c(0,240),breaks=seq(0,220,by=20))+
  theme_classic()

SQ<-nests[nests$sp %in% c("SPM","STE","STP"),]

ggplot(SQ,aes(y=Z,x=as.factor(year),fill=as.factor(year)))+
  geom_boxplot(outlier.shape=NA)+
  geom_point(size=5,position=position_jitterdodge())+
  scale_y_continuous(limits=c(0,240),breaks=seq(0,220,by=20))+
  theme_classic()

ggplot(SQ,aes(y=Z, x=as.factor(year),shape=stage,fill=as.factor(sp)))+
  geom_boxplot(outlier.shape=NA)+
  geom_point(size=5,position=position_jitterdodge())+
  scale_shape_manual(values=c(21))+
  scale_y_continuous(limits=c(0,240),breaks=seq(0,220,by=20))+
  theme_classic()

### calculating medians and IQR for birds

XYZ<-read.csv("XYZ.csv")
str(XYZ)


quantile(subset(LS,stage=="nest"&year==18)$Z)
quantile(subset(LS,stage=="nest"&year==19)$Z)
quantile(subset(LS,stage=="brood"&year==18)$Z)
quantile(subset(LS,stage=="brood"&year==19)$Z)
quantile(subset(RK,stage=="nest"&year==18)$Z)
quantile(subset(RK,stage=="nest"&year==19)$Z)
quantile(subset(RK,stage=="brood"&year==18)$Z)
quantile(subset(RK,stage=="brood"&year==19)$Z)
quantile(subset(SQ,year==18)$Z)
quantile(subset(SQ,year==19)$Z)

### The median (interquartile range) elevation of the nests observed 
### was 41.0 (37.1 - 53.5) m a.s.l. for Little Stints and 103.0 (95.0 - 117.75) m a.s.l. 
### for Red Knots. The median (interquartile range) elevation of broods caught 
### was 46.1 (39.3 - 55.0) m a.s.l. for Little Stints and 123.0 (103.0 - 139.9) m a.s.l.
### for Red Knots.

### Mann-Witney test: nests vs broods

LS$year<-as.factor(LS$year)
wilcox.test(subset(LS,stage=="nest"&year=="19")$Z,subset(LS,stage=="brood"&year=="19")$Z,paired=FALSE)
wilcoxonZ(subset(LS,stage=="nest"&year=="19")$Z,subset(LS,stage=="brood"&year=="19")$Z)

RK$year<-as.factor(RK$year)
wilcox.test(subset(RK,stage=="nest"&year=="18")$Z,subset(RK,stage=="brood"&year=="18")$Z,paired=FALSE)
wilcoxonZ(subset(RK,stage=="nest"&year=="18")$Z,subset(RK,stage=="brood"&year=="18")$Z)

wilcox.test(subset(RK,stage=="nest"&year=="19")$Z,subset(RK,stage=="brood"&year=="19")$Z,paired=FALSE)
wilcoxonZ(subset(RK,stage=="nest"&year=="19")$Z,subset(RK,stage=="brood"&year=="19")$Z)

### calculating per elevation class

per_elevation_class<-read.csv("XYZ.csv")
per_elevation_class<-uncount(per_elevation_class,N)
per_elevation_class$Z_cat<-as.factor(per_elevation_class$Z_cat)
table(subset(per_elevation_class,animal=="Lemming"&year==18)$Z_cat)
table(subset(per_elevation_class,animal=="Lemming"&year==19)$Z_cat)
table(subset(per_elevation_class,animal=="Squa"&year==18)$Z_cat)
table(subset(per_elevation_class,animal=="Squa"&year==19)$Z_cat)
table(subset(per_elevation_class,animal=="LS"&year==18&stage=="nest")$Z_cat)
table(subset(per_elevation_class,animal=="LS"&year==19&stage=="nest")$Z_cat)
table(subset(per_elevation_class,animal=="LS"&year==18&stage=="brood")$Z_cat)
table(subset(per_elevation_class,animal=="LS"&year==19&stage=="brood")$Z_cat)
table(subset(per_elevation_class,animal=="RK"&year==18&stage=="nest")$Z_cat)
table(subset(per_elevation_class,animal=="RK"&year==19&stage=="nest")$Z_cat)
table(subset(per_elevation_class,animal=="RK"&year==18&stage=="brood")$Z_cat)
table(subset(per_elevation_class,animal=="RK"&year==19&stage=="brood")$Z_cat)

############################################################################################
#                                                                                          #
#                                       PREFERENCE SCORE                                   #
#                                                                                          #
############################################################################################

### plotting preference score (Ivlev Index) for red knots, little stints, lemmings and squas

Ivlev<-read.csv("Ivlev_calculation.csv")
str(Ivlev)
Ivlev$Ivlev_index<-as.numeric(Ivlev$Ivlev_index)

ggplot(subset(Ivlev,type=="LS_brood"|type=="LS_nest"),aes(x=elevation,y=Ivlev_index,shape=type,fill=as.factor(year)),label=year)+
  geom_point(size=8)+
  scale_shape_manual(values=c(24,21))+
  scale_x_continuous(limits=c(0,220),breaks=seq(0,220,by=20))+
  scale_y_continuous(limits=c(-1,1),breaks=seq(-1,1,by=0.5))+
  geom_hline(yintercept=0,linetype="solid",size=.1)+
  geom_text(aes(label=year),size=3,color="white")+
  theme_classic()

ggplot(subset(Ivlev,type=="RK_brood"|type=="RK_nest"),aes(x=elevation,y=Ivlev_index,shape=type,fill=as.factor(year)),label=year)+
  geom_point(size=8)+
  scale_shape_manual(values=c(24,21))+
  scale_x_continuous(limits=c(0,220),breaks=seq(0,220,by=20))+
  scale_y_continuous(limits=c(-1,1),breaks=seq(-1,1,by=0.5))+
  geom_hline(yintercept=0,linetype="solid",size=.1)+
  geom_text(aes(label=year),size=3,color="white")+
  theme_classic()

ggplot(subset(Ivlev,type=="Lemming"),aes(x=elevation,y=Ivlev_index,fill=as.factor(year)),label=year)+
  geom_point(size=8,shape=21)+
  scale_x_continuous(limits=c(0,220),breaks=seq(0,220,by=20))+
  scale_y_continuous(limits=c(-1,1),breaks=seq(-1,1,by=0.5))+
  geom_hline(yintercept=0,linetype="solid",size=.1)+
  geom_text(aes(label=year),size=3,color="white")+
  theme_classic()

ggplot(subset(Ivlev,type=="Squa"),aes(x=elevation,y=Ivlev_index,fill=as.factor(year)),label=year)+
  geom_point(size=8,shape=21)+
  scale_x_continuous(limits=c(0,220),breaks=seq(0,220,by=20))+
  scale_y_continuous(limits=c(-1,1),breaks=seq(-1,1,by=0.5))+
  geom_hline(yintercept=0,linetype="solid",size=.1)+
  geom_text(aes(label=year),size=3,color="white")+
  theme_classic()

### Little Stints were attracted to lower elevation for both building nests and rearing chicks, as shown by positive preference score (Ivlev Index > 0) for elevations between 20 and 60 m a.s.l.
### The maximum score for nesting was observed at elevation between 20 and 40 m a.s.l. (Ivlev Index = 0.54), while for broods it was between 40 and 60 m a.s.l. (Ivlev Index = 0.55).
### No nests or broods were encountered above 140 m a.s.l. during any of the years (Ivlev Index = -1.00).
### The elevations between 80 and 160 m a.s.l. were prefered by Red Knots based on locations of nests and broods during the season (Ivlev Index > 0).
### The highest preference was observed between 100 and 120 m a.s.l. for nests (Ivlev Index = 0.52) and between 140 and 160 m a.s.l. for broods (Ivlev Index = 0.62).
### Elevations lower than 80 m a.s.l. were averted by Red Knots during the entire breeding period (Ivlev Index < 0).

############################################################################################
#                                                                                          #
#                                         ARTHROPODS                                       #
#                                                                                          #
############################################################################################

### arthropods number per trap per day vs elevation

chi_and_tip<-read.csv('Tip_and_Chi_n.csv')
str(chi_and_tip)

ggplot(subset(chi_and_tip,year==2018),aes(x=elevation,y=(n+1)/n_sampling_days,fill=family))+
  geom_smooth(method="lm",alpha=.4,color="black")+
  geom_point(size=5,shape=21)+
  scale_fill_manual(values=c("grey25","grey75"))+
  theme_classic()+
  scale_y_log10(breaks=c(0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,2,3,4,5,6,7,8,9,10,20,30,40,50,60,70,80,90,100,200,300),limits=c(0.02,8))+
  xlim(40,125)

ggplot(subset(chi_and_tip,year==2019),aes(x=elevation,y=(n+1)/n_sampling_days,fill=family))+
  geom_smooth(method="lm",alpha=.4,color="black")+
  geom_point(size=5,shape=21)+
  scale_fill_manual(values=c("grey25","grey75"))+
  theme_classic()+
  scale_y_log10(breaks=c(0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,2,3,4,5,6,7,8,9,10,20,30,40,50,60,70,80,90,100,200,300),limits=c(0.02,8))+
  xlim(40,125)

###model: Chironomidae ~ elvation X year 
###vif

vif(lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=log10((n+1)/n_sampling_days)~elevation*as.factor(year)))
vif(lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=log10((n+1)/n_sampling_days)~elevation+as.factor(year)))

###checking assumptions

par(mfrow = c(2, 2))
plot(lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=log10((n+1)/n_sampling_days)~elevation*as.factor(year)))

par(mfrow = c(2, 2))
plot(lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=n_sampling_days~elevation*as.factor(year)))

par(mfrow = c(2, 2))
plot(lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=log10((n+1)/n_sampling_days)~elevation*as.factor(year)))

par(mfrow = c(2, 2))
plot(lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=n_sampling_days~elevation*as.factor(year)))

summary(lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=log10((n+1)/n_sampling_days)~elevation*as.factor(year)))
summary(lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=log10((n+1)/n_sampling_days)~elevation+as.factor(year)))
summary(lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=log10((n+1)/n_sampling_days)~elevation))
summary(lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=log10((n+1)/n_sampling_days)~as.factor(year)))
summary(lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=log10((n+1)/n_sampling_days)~1))

AICc(lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=log10((n+1)/n_sampling_days)~elevation*as.factor(year)))
AICc(lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=log10((n+1)/n_sampling_days)~elevation+as.factor(year)))
AICc(lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=log10((n+1)/n_sampling_days)~elevation))
AICc(lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=log10((n+1)/n_sampling_days)~as.factor(year)))
AICc(lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=log10((n+1)/n_sampling_days)~1))

50.96575-39.74495
79.14883-39.74495
68.48257-39.74495
87.29563-39.74495
Weights(AICc(lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=log10((n+1)/n_sampling_days)~elevation*as.factor(year)),lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=log10((n+1)/n_sampling_days)~elevation+as.factor(year)),lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=log10((n+1)/n_sampling_days)~elevation),lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=log10((n+1)/n_sampling_days)~as.factor(year)),lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=log10((n+1)/n_sampling_days)~1)))

logLik(lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=log10((n+1)/n_sampling_days)~elevation*as.factor(year)))
logLik(lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=log10((n+1)/n_sampling_days)~elevation+as.factor(year)))
logLik(lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=log10((n+1)/n_sampling_days)~elevation))
logLik(lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=log10((n+1)/n_sampling_days)~as.factor(year)))
logLik(lm(data=subset(chi_and_tip,family=="Chironomidae"),formula=log10((n+1)/n_sampling_days)~1))

##Tiipulidae

##vif

vif(lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=log10((n+1)/n_sampling_days)~elevation*as.factor(year)))
vif(lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=log10((n+1)/n_sampling_days)~elevation+as.factor(year)))

summary(lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=log10((n+1)/n_sampling_days)~elevation*as.factor(year)))
summary(lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=log10((n+1)/n_sampling_days)~elevation+as.factor(year)))
summary(lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=log10((n+1)/n_sampling_days)~elevation))
summary(lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=log10((n+1)/n_sampling_days)~as.factor(year)))
summary(lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=log10((n+1)/n_sampling_days)~1))

AICc(lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=log10((n+1)/n_sampling_days)~elevation*as.factor(year)))
AICc(lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=log10((n+1)/n_sampling_days)~elevation+as.factor(year)))
AICc(lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=log10((n+1)/n_sampling_days)~elevation))
AICc(lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=log10((n+1)/n_sampling_days)~as.factor(year)))
AICc(lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=log10((n+1)/n_sampling_days)~1))

100.1398-76.21432
100.6706-76.21432
78.48956-76.21432
77.73353-76.21432

Weights(AICc(lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=log10((n+1)/n_sampling_days)~elevation*as.factor(year)),lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=log10((n+1)/n_sampling_days)~elevation+as.factor(year)),lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=log10((n+1)/n_sampling_days)~elevation),lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=log10((n+1)/n_sampling_days)~as.factor(year)),lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=log10((n+1)/n_sampling_days)~1)))

logLik(lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=log10((n+1)/n_sampling_days)~elevation*as.factor(year)))
logLik(lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=log10((n+1)/n_sampling_days)~elevation+as.factor(year)))
logLik(lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=log10((n+1)/n_sampling_days)~elevation))
logLik(lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=log10((n+1)/n_sampling_days)~as.factor(year)))
logLik(lm(data=subset(chi_and_tip,family=="Tipulidae"),formula=log10((n+1)/n_sampling_days)~1))

summary(lm(data=subset(chi_and_tip,family=="Chironomidae"&year==2018),formula=log10((n+1)/n_sampling_days)~elevation))
summary(lm(data=subset(chi_and_tip,family=="Chironomidae"&year==2019),formula=log10((n+1)/n_sampling_days)~elevation))

summary(lm(data=subset(chi_and_tip,family=="Tipulidae"&year==2018),formula=log10((n+1)/n_sampling_days)~elevation))
summary(lm(data=subset(chi_and_tip,family=="Tipulidae"&year==2019),formula=log10((n+1)/n_sampling_days)~elevation))

### comparison of 40 and 120 m a.s.l. based on the models outcome

10^(0.641-0.014*40)/10^(0.641-0.014*120)

10^(0.375-0.004*40)/10^(0.375-0.004*120)

1/(10^(-1.046+0.007*40)/10^(-1.046+0.007*120))

1/(10^(-1.495+0.010*40)/10^(-1.495+0.010*120))

############################################################################################
#                                                                                          #
#                                       DIET EXPLANATION                                   #
#                                                                                          #
############################################################################################

### % of arthropod in diet (preliminary estimation)

diet<-read.csv("diet.csv")
diet_LS<-subset(diet,bird=="LS")
diet_RK<-subset(diet,bird=="RK")

ggplot(diet_LS,aes(x=2,y=percentage))+
  geom_col(color="black")+
  geom_text(aes(label=paste0(round(percentage,digits=1),"%")),position=position_stack(vjust=.5))+
  coord_polar("y",start=1)+
  theme_classic()+
  xlim(.01,2.5)

ggplot(diet_RK,aes(x=2,y=percentage))+
  geom_col(color="black")+
  geom_text(aes(label=paste0(round(percentage,digits=1),"%")),position=position_stack(vjust=.5))+
  coord_polar("y",start=1)+
  theme_classic()+
  xlim(0.01,2.5)

### selection of the best model explaining the proportion of Tipulidae in the diet (full model: Tipulidae in the diet ~ elevation X species)

df<-read.csv("faeces_data.csv")
str(df)
df<-subset(df,family=="Tipulidae"&(sp=="RK"|sp=="LS")&year==2019)
nrow(df)
df$reads_ZI<-(53^-1)*(df$reads*(53-1)+0.5)
df<-df[complete.cases(df), ]

df<-aggregate(df$reads_ZI,list(df$sp,df$brood,df$elevation),FUN="mean")
names(df)<-c("sp","brood","elevation","reads_ZI")
str(df)

plot(df)

m_full<-betareg(formula=reads_ZI~elevation*sp,data=df)
summary(m_full)

m_1<-betareg(formula=reads_ZI~elevation+sp,data=df)

vif(m_1)

m_2<-betareg(formula=reads_ZI~elevation,data=df)
m_3<-betareg(formula=reads_ZI~sp,data=df)
m_4<-betareg(formula=reads_ZI~1,data=df)

sum_all<-summary(model.avg(m_full,m_1,m_2,m_3,m_4))
sum_all
write.table(sum_all$msTable,file="models_Tip.csv",sep=",")

summary(m_2)

### plotting the best model: Crane flies ~ elevation

ggplot(data=df,aes(y=reads_ZI,x=elevation))+
  geom_smooth(method = 'glm',  method.args = list(family = "quasibinomial"))+
  geom_point(shape=21,size=3,aes(fill=sp))+
  theme_classic()

############################################################################################
#                                                                                          #
#                                           LEMMINGS                                       #
#                                                                                          #
############################################################################################

### lemmings distribution, preference and abundance

### new

lemmings<-subset(read.csv("XYZ.csv"),animal=="Lemming")
lemmings<-uncount(lemmings,N)
str(lemmings)
nrow(lemmings)

ggplot(lemmings,aes(y=Z,x=as.factor(year),fill=as.factor(year)))+
  geom_hline(yintercept=108.5)+
  geom_hline(yintercept=103)+
  geom_hline(yintercept=39.6)+
  geom_hline(yintercept=44.2)+
  geom_boxplot(outlier.shape=NA)+
  geom_point(size=5,shape=21,position=position_jitterdodge())+
  scale_y_continuous(limits=c(0,240),breaks=seq(0,220,by=20))+
  theme_classic()

quantile(subset(lemmings,year==18)$Z)
quantile(subset(lemmings,year==19)$Z)

lemmings$year<-as.factor(lemmings$year)
wilcox.test(subset(lemmings,year=="18")$Z,subset(lemmings,year=="19")$Z,paired=FALSE)
wilcoxonZ(subset(lemmings,year=="18")$Z,subset(lemmings,year=="19")$Z)

### lemmings per day

lemmings_per_day<-read.csv('Lemmings_per_day.csv')
str(lemmings_per_day)
lemmings_per_day$year<-as.factor(lemmings_per_day$year)

ggplot(lemmings_per_day,aes(y=n,x=year,fill=year))+
  geom_boxplot(outlier.shape=NA)+
  geom_point(size=5,shape=21,position=position_jitterdodge())+
  scale_fill_manual(values=c("orange","lightblue"))+
  theme_classic()+
  scale_y_continuous(limits=c(0,30),breaks=seq(0,30,by=5))

quantile(na.omit(subset(lemmings_per_day,year=="2018")$n))
quantile(na.omit(subset(lemmings_per_day,year=="2019")$n))

mean(na.omit(subset(lemmings_per_day,year=="2018")$n))
max(na.omit(subset(lemmings_per_day,year=="2018")$n))
mean(na.omit(subset(lemmings_per_day,year=="2019")$n))
max(na.omit(subset(lemmings_per_day,year=="2019")$n))

wilcox.test(na.omit(lemmings_per_day$n~lemmings_per_day$year),paired=FALSE)

wilcoxonZ(subset(lemmings_per_day,year=="2018")$n,subset(lemmings_per_day,year=="2019")$n)

### The total number of lemmings observed was 30 during 50 field days in 2018 and 
### 226 during 38 field days in 2019. The abundance of lemmings was lower in 2018 
### with median (interquartile range) encounters per day of 0 (0.0 - 1.0) against 
### 2.0 (0.3 - 8.5) in 2019 (U-test, p < 0.001, Z = -4.47).

############################################################################################
#                                                                                          #
#                                CORRELATION BETWEEN PREFERENCES                           #
#                                                                                          #
############################################################################################

cor_ivlev_both<-read.csv("for_correlation.csv")

### correlation between lemmings and shorebirds habitat preference

ggpar(ylim=c(-1,1),xlim=c(-1,1),
      ggscatter(subset(cor_ivlev_both,year==2018), x = "Lemming", y = "LS_nest", 
                add = "reg.line", conf.int = FALSE, 
                cor.coef = TRUE, cor.method = "spearman"))

ggpar(ylim=c(-1,1),xlim=c(-1,1),
      ggscatter(subset(cor_ivlev_both,year==2018), x = "Lemming", y = "RK_nest", 
                add = "reg.line", conf.int = FALSE, 
                cor.coef = TRUE, cor.method = "spearman"))

ggpar(ylim=c(-1,1),xlim=c(-1,1),
      ggscatter(subset(cor_ivlev_both,year==2019), x = "Lemming", y = "LS_nest", 
                add = "reg.line", conf.int = FALSE, 
                cor.coef = TRUE, cor.method = "spearman"))

ggpar(ylim=c(-1,1),xlim=c(-1,1),
      ggscatter(subset(cor_ivlev_both,year==2019), x = "Lemming", y = "RK_nest", 
                add = "reg.line", conf.int = FALSE, 
                cor.coef = TRUE, cor.method = "spearman"))

### correlation between squas and shorebirds habitat preference

ggpar(ylim=c(-1,1),xlim=c(-1,1),
      ggscatter(subset(cor_ivlev_both,year==2018), x = "LS_nest", y = "Squa", 
                add = "reg.line", conf.int = FALSE, 
                cor.coef = TRUE, cor.method = "spearman"))

ggpar(ylim=c(-1,1),xlim=c(-1,1),
      ggscatter(subset(cor_ivlev_both,year==2019), x = "LS_nest", y = "Squa", 
                add = "reg.line", conf.int = FALSE, 
                cor.coef = TRUE, cor.method = "spearman"))

ggpar(ylim=c(-1,1),xlim=c(-1,1),
      ggscatter(subset(cor_ivlev_both,year==2018), x = "RK_nest", y = "Squa", 
                add = "reg.line", conf.int = FALSE, 
                cor.coef = TRUE, cor.method = "spearman"))

ggpar(ylim=c(-1,1),xlim=c(-1,1),
      ggscatter(subset(cor_ivlev_both,year==2019), x = "RK_nest", y = "Squa", 
                add = "reg.line", conf.int = FALSE, 
                cor.coef = TRUE, cor.method = "spearman"))

## cor between squas and lemmings

ggpar(ylim=c(-1,1),xlim=c(-1,1),
      ggscatter(subset(cor_ivlev_both,year==2018), x = "Lemming", y = "Squa", 
                add = "reg.line", conf.int = FALSE, 
                cor.coef = TRUE, cor.method = "spearman"))

ggpar(ylim=c(-1,1),xlim=c(-1,1),
      ggscatter(subset(cor_ivlev_both,year==2019), x = "Lemming", y = "Squa", 
                add = "reg.line", conf.int = FALSE, 
                cor.coef = TRUE, cor.method = "spearman"))

############################################################################################
#                                                                                          #
#                                       MAP CALLIBRATION                                   #
#                                                                                          #
############################################################################################

### supplementary GPS vs topographic map. Accuracy check

check<-read.csv("Taymir_reserve_map_vs_ArcDEM_accuracy.csv")
str(check)

ggscatter(check,x="Z_cat_ArcDEM",y="Z_cat_map",cor.method = "spearman",add="reg.line",conf.int=TRUE,cor.coef = TRUE, )

ggplot(check,aes(y=Z_cat_map,x=Z_cat_ArcDEM))+
  geom_vline(xintercept=0,colour="lightgrey")+
  geom_vline(xintercept=20,colour="lightgrey")+
  geom_vline(xintercept=40,colour="lightgrey")+
  geom_vline(xintercept=60,colour="lightgrey")+
  geom_vline(xintercept=80,colour="lightgrey")+
  geom_vline(xintercept=100,colour="lightgrey")+
  geom_vline(xintercept=120,colour="lightgrey")+
  geom_vline(xintercept=140,colour="lightgrey")+
  geom_vline(xintercept=160,colour="lightgrey")+
  geom_vline(xintercept=180,colour="lightgrey")+
  geom_vline(xintercept=200,colour="lightgrey")+
  geom_vline(xintercept=220,colour="lightgrey")+
  geom_hline(yintercept=0,colour="lightgrey")+
  geom_hline(yintercept=20,colour="lightgrey")+
  geom_hline(yintercept=40,colour="lightgrey")+
  geom_hline(yintercept=60,colour="lightgrey")+
  geom_hline(yintercept=80,colour="lightgrey")+
  geom_hline(yintercept=100,colour="lightgrey")+
  geom_hline(yintercept=120,colour="lightgrey")+
  geom_hline(yintercept=140,colour="lightgrey")+
  geom_hline(yintercept=160,colour="lightgrey")+
  geom_hline(yintercept=180,colour="lightgrey")+
  geom_hline(yintercept=200,colour="lightgrey")+
  geom_hline(yintercept=220,colour="lightgrey")+
  geom_abline(slope=1,intercept=0)+
  geom_smooth(method="lm",alpha=.4,fill="darkgrey",color="black")+
  geom_point(shape=21,size=4,fill="grey25")+
  theme_classic()+
  scale_x_continuous(limits=c(0,220),breaks=seq(0,220,by=20))+
  scale_y_continuous(limits=c(0,220),breaks=seq(0,220,by=20))

table(check[,c(5,6)])

### citation

citation()