## library("MASS")
## library("magrittr")
## library("QuantPsyc")
## lm.beta(nb)


library("MASS")
library("magrittr")
library("ggplot2")
install.packages('TMB', type = 'source')

#Generalized linear models of responses expresed as proportions with with binomial family and logit link functions
data <- Norway_Sweden_GR_plant_data_1_

#Total plant cover GLM
glm.TPC <- glm(data$`prop TPC` ~ data$`T mean`+data$Precipitation+data$LDP+data$FTC +data$`WHL+Sub`+data$Slope, family = binomial(link="logit"))
summary(glm.TPC)

res.tpc <- residuals.glm(glm.TPC)


plot(res.tpc)   
#Vascular plant cover GLM
glm.Vasc <- glm(data$`prop vasc` ~ data$`T mean`+data$Precipitation+data$LDP+data$FTC +data$`WHL+Sub`+data$Slope, family = binomial(link="logit"))
summary(glm.Vasc)

res.Vasc <- residuals.glm(glm.Vasc)

plot(res.Vasc)

#Intended vegetation cover GLM
glm.intended <- glm(data$`prop Intended cover` ~ data$`T mean`+data$Precipitation+data$LDP+data$FTC +data$`WHL+Sub`+data$Slope, family = binomial(link="logit"))
summary(glm.sedum)

res.intended <- residuals(glm.intended)

plot(res.intended)

#Unintended cover GLM
glm.unintended <- glm(data$`prop unintended` ~  data$`T mean`+data$Precipitation+data$LDP+data$FTC +data$`WHL+Sub`+data$Slope, family = binomial(link="logit"))
summary(glm.unintended)

plot(data$Slope, data$`prop bare`)

#Moss cover GLM
glm.Moss <- glm(data$`prop moss` ~ data$`T mean`+data$Precipitation+data$LDP+data$FTC +data$`WHL+Sub`+data$Slope, family = binomial(link="logit"))
summary(glm.Moss)

res.Moss <- residuals(glm.Moss)
plot(res.Moss)

#Bare substrate GLM
glm.bare <- glm(data$`prop bare` ~  data$`T mean`+data$Precipitation+data$LDP+data$FTC +data$`WHL+Sub`+data$Slope, family = binomial(link="logit"))
summary(glm.bare)

res.bare <- residuals(glm.bare)
plot(res.bare)


#Proportion "refound" species GLM
glm.refound <- glm(data$`prop refound` ~  data$`T mean`+data$Precipitation+data$LDP+data$FTC +data$`WHL+Sub`+data$Slope, family = binomial(link="logit"))
summary(glm.Succ)

res.Succ <- residuals(glm.Succ)
plot(res.Succ)

library(dplyr) #basic data management & %>%
library(MASS) #for negative bin regression
library(magrittr)
library(stargazer) #for combining output
library(ggplot2) #for graphing
library(pscl) #for zero inflated models & predprob
library(reshape2) #convert wide to tall
library(summarytools) #freq
library(jtools)
library(ggplot2)
library(forcats)


#Visualisations
#Figure 3 box plots of vegetations and bare substrate
data <- Norway_Sweden_GR_plant_data_1_
cbind(data$Location, data$year)
cbind(data$`Species mix`, data$year)
data$`Species mix`<- as.factor(data$`Species mix`)
data$`Species mix` <- as.numeric(data$`Species mix`)

summary(data)

year<-factor(data$year)

#Intended cover by location boxplot
Intended_bp <- ggplot(data=Norway_Sweden_GR_plant_data_1_, aes(x=fct_reorder(data$Location, data$Latitude), y=data$`Intended cover`, fill=factor(year)))+ geom_boxplot(lwd=0.5, position=position_dodge2(preserve = "single"))+ ylim(0,100)+theme(legend.position = 'bottom')+scale_fill_brewer()+ylab("% intended cover")+xlab("Location")+labs(fill="year")+
  theme(legend.position="none", axis.title.x = element_blank(), axis.text.x = element_text(color="black", size=6), axis.title.y = element_text(color="black", size=10), axis.text.y = element_text(color="black", size=6))

Intended_bp

#Unintended cover by location boxplot
Unintended_bp <- ggplot(data=Norway_Sweden_GR_plant_data_1_, aes(x=fct_reorder(data$Location, data$Latitude), y=data$`Unintended cover`, fill=factor(year)))+ geom_boxplot(lwd=0.5, position=position_dodge2(preserve = "single"))+ ylim(0,100)+theme(legend.position = 'bottom')+scale_fill_brewer()+ylab("% unintedned cover")+xlab("Location")+labs(fill="year")+
  theme(legend.position="none", axis.title.x = element_blank(), axis.text.x = element_text(color="black", size=6), axis.title.y = element_text(color="black", size=10), axis.text.y = element_text(color="black", size=6))

Unintended_bp

#Moss cover by location boxplot
Moss_bp <- ggplot(data=Norway_Sweden_GR_plant_data_1_, aes(x=fct_reorder(Location, data$Latitude), y=data$Moss, fill=factor(year)))+ geom_boxplot(lwd=0.5, position=position_dodge2(preserve = "single"))+ ylim(0,100)+theme(legend.position = 'bottom')+scale_fill_brewer()+ylab("% moss cover") +xlab("Location")+labs(fill="year")+
  theme(legend.position="none", axis.title.x = element_blank(), axis.text.x = element_text(color="black", size=6), axis.title.y = element_text(color="black", size=10), axis.text.y = element_text(color="black", size=6))

Moss_bp

#Bare substrate by location cover boxplot
Bare_bp <- ggplot(data=Norway_Sweden_GR_plant_data_1_, aes(x=fct_reorder(Location, data$Latitude), y=`Bare substr`, fill=factor(year)))+ geom_boxplot(lwd=0.5, position=position_dodge2(preserve = "single"))+ ylim(0,100)+theme(legend.position = 'bottom')+scale_fill_brewer()+ylab("% bare substrate") + xlab("Location")+labs(fill="year")+
  theme(legend.position="none", axis.title.x = element_blank(), axis.text.x = element_text(color="black", size=6), axis.title.y = element_text(color="black", size=10), axis.text.y = element_text(color="black", size=6))

Bare_bp

#Intended cover by species mixture
Intended_by_mix_bp <- ggplot(data=Norway_Sweden_GR_plant_data_1_, aes(x=data$`Species mix`, y=data$`Intended cover`, fill=factor(year)))+ geom_boxplot(lwd=0.5, position=position_dodge2(preserve = "single"))+ ylim(0,100)+theme(legend.position = 'bottom')+scale_fill_brewer()+ylab("% intended cover") + xlab("Species mixture")+labs(fill="year")+
  theme(legend.position="bottom", axis.title.x = element_text(color="black", size=10), axis.text.x = element_text(color="black", size=6), axis.title.y = element_text(color="black", size=10), axis.text.y = element_text(color="black", size=10))

Intended_by_mix_bp

#Species richness by location
Richness_plt <- ggplot(data, aes(x=fct_reorder(Location, Latitude), y=`Species richness`, fill=factor(year)))+ geom_boxplot(lwd=0.5, position=position_dodge2(preserve = "single"))+ ylim(0,8)+theme(legend.position = 'bottom')+scale_fill_brewer()+ylab("Species richness") + xlab("Location")+labs(fill="year")
Richness_plt

#All covers merged in one plot
percent_covers <- ggarrange(Intended_bp, Unintended_bp, Moss_bp,Bare_bp, 
                            labels = c("A", "B", "C", "D"),
                            ncol = 2, nrow = 2, legend="bottom", common.legend=TRUE)

percent_covers



#Shannon_bp <- ggplot(data, aes(x=fct_reorder(Location, data$Latitude), y=data$`Shannon div index H`, fill=factor(year)))+ geom_boxplot(lwd=0.5, position=position_dodge2(preserve = "single"))+ ylim(0,3)+theme(legend.position = 'bottom')+scale_fill_viridis_d()+ylab("Shannon diversity index")+xlab("Location")+labs(fill="year")
Shannon_bp

data=Norway_Sweden_GR_plant_data_1_

data.f<- data.frame(data$Location,data$Latitude,data$`Originally planted`,data$Refound)

library(reshape2)

dat.m <-melt(data.f,idvars=cbind("Location", Latitude,`Originally planted`,Refound))
dat.m


my_comparisons <- list( c("`Originally intended`", "2016"), c("2016", "2017"), c("`Originally intended`", "2017") )

Planted_and_refound <- ggplot(b_a, aes(x=fct_reorder(Location, Latitude), y=Species, fill=fct_reorder(factor(Time),Dummy)))+ geom_boxplot(lwd=0.5, position=position_dodge2(preserve = "single"))+ ylim(0,12.5)+theme(legend.position = 'bottom')+scale_fill_brewer()+ylab("Number of species") + xlab("Location")+labs(fill="")+theme(legend.position="bottom", axis.title.x = element_blank(), axis.text.x = element_text(color="black", size=8), axis.title.y = element_text(color="black", size=10), axis.text.y = element_text(color="black", size=10))


Planted_and_refound

levels(b_a$Dummy)
b_a$Dummy <- ordered(b_a$Dummy,
                     levels = c("1", "2", "3"))

Location <- data.frame(cbind(b_a$Location))
Location

b_a$Time=factor(b_a$Time, levels=c("`Originally planted`", "2016", "2017"))

#KRUSKAL TEST
compare_means(Species~Time, b_a)
Sandnes_krusk <- kruskal.test(b_a$Species ~ b_a$Dummy, subset =b_a$Location=="Sandnes")
Sandnes_krusk

summary(Sandnes_krusk)

wilcox.paired

dt_ba <- dunn.test(b_a$Species, b_a$Time, method="bonferroni", kw=TRUE, label=TRUE, wrap=TRUE, table=TRUE, list=TRUE)

#Pairwise Wilcox test
Sandnes_pwxt <- pairwise.wilcox.test(Sandnes$Species, Sandnes$Time, paired=TRUE, p.adjust.method = "bonferroni")
Sandnes_pwxt

Drammen_pwxt <- pairwise.wilcox.test(Drammen$Species, Drammen$Time, paired=TRUE, p.adjust.method = "bonferroni")
Drammen_pwxt

Oslo1_pwxt <- pairwise.wilcox.test(Oslo1$Species, Oslo1$Time, paired=TRUE, p.adjust.method = "bonferroni")
Oslo1_pwxt

Oslo2_pwxt <- pairwise.wilcox.test(Oslo2$Species, Oslo2$Time, paired=TRUE, p.adjust.methods)
Oslo2_pwxt

Bergen_pwxt <- pairwise.wilcox.test(Bergen$Species, Bergen$Time, paired=TRUE, p.adjust.method = "bonferroni")
Bergen_pwxt

Trondheim_pwxt <- pairwise.wilcox.test(Trondheim$Species, Trondheim$Time, paired=TRUE, p.adjust.method = "bonferroni")
Trondheim_pwxt

Umea_pwxt <- pairwise.wilcox.test(Umea$Species, Umea$Time, p.adjust.method = "bonferroni")
Umea_pwxt

Lulea_pwxt <- pairwise.wilcox.test(Lulea$Species, Lulea$Time, p.adjust.method = "bonferroni")
Lulea_pwxt


Kiruna_pwxt <- pairwise.wilcox.test(Kiruna$Species, Kiruna$Time, p.adjust.method = "bonferroni")
Kiruna_pwxt


Tromso_pwxt <- pairwise.wilcox.test(Tromso$Species, Tromso$Time, p.adjust.method = "bonferroni")
Tromso_pwxt



#Weather data compiled

library(tidyverse)

ggplot(data=data, aes(x=Location, y=data$`T mean`, fill=year)) +  geom_bar()

H <- c(data$Location,data$`T mean`, data$year)
M <- c(data$Location)

weather <- cbind(Summary_weather_data_2)
weather
year=c(2017,2016)

df.long<-melt(df,id.vars=c("ID","Type","year"))
df.long$Location<-substring(df.long$Location,2)

Temp_m <- ggplot(weather,aes(x=fct_reorder(weather$Location, weather$Latitude),y=weather$`Mean temperature`,fill=as.factor(year)))+geom_bar(color="black", size=0.5,position=position_dodge2(preserve = "single"),stat="identity")+theme(legend.position = 'bottom')+scale_fill_brewer()+ylab("Temperature mean °C ") + xlab("Location")+labs(fill="")+
  theme(legend.position="none", axis.title.x = element_blank(), axis.text.x = element_text(color="black", size=5), axis.title.y = element_text(color="black", size=10), axis.text.y = element_text(color="black", size=6))

Temp_m

FTC_bc <- ggplot(weather,aes(x=fct_reorder(Location, Latitude),weather$FTC,fill=as.factor(year)))+geom_bar(color="black", size=0.5,position=position_dodge2(preserve = "single"),stat="identity")+theme(legend.position = 'bottom')+scale_fill_brewer()+ylab("Freeze thaw cycles") + xlab("Location")+labs(fill="")+
  theme(legend.position="none", axis.title.x = element_blank(), axis.text.x = element_text(color="black", size=5), axis.title.y = element_text(color="black", size=10), axis.text.y = element_text(color="black", size=6))

FTC_bc


Total_P_bc <- ggplot(weather,aes(x=fct_reorder(Location, Latitude),weather$`total P`,fill=as.factor(year)))+geom_bar(color="black", size=0.5, position=position_dodge2(preserve = "single"),stat="identity")+theme(legend.position = 'bottom')+scale_fill_brewer()+ylab("Total precipitation (mm)") + xlab("Location")+labs(fill="")+
  theme(legend.position="none", axis.title.x = element_blank(), axis.text.x = element_text(color="black", size=5), axis.title.y = element_text(color="black", size=10), axis.text.y = element_text(color="black", size=6))
Total_P_bc

LDP_bc <- ggplot(weather,aes(x=fct_reorder(Location, Latitude),weather$LDP,fill=as.factor(year)))+geom_bar(color="black", size=0.5,position=position_dodge2(preserve = "single"),stat="identity")+theme(legend.position = 'bottom')+scale_fill_brewer()+ylab("Longest dought (days)") + xlab("Location")+labs(fill="")+
  theme(legend.position="none", axis.title.x = element_blank(), axis.text.x = element_text(color="black", size=5), axis.title.y = element_text(color="black", size=10), axis.text.y = element_text(color="black", size=6))
LDP_bc

weather_bp <- ggarrange(Temp_m, FTC_bc,Total_P_bc, LDP_bc, 
                        labels = c("A", "B", "C", "D"),
                        ncol = 2, nrow = 2, legend="bottom", common.legend=TRUE)

weather_bp