rm(list=ls(all=TRUE))
options(stringsAsFactors=F)

library (ggplot2)
library(DHARMa)
library(lme4)
library(ggcorrplot)
library(ggplot2)
library(GGally)
library(patchwork)
library(RColorBrewer)
library(tidyverse)
library(car)
library(MuMIn)
df <- read.csv("data/NF_data_submitted.csv")

#############################
#model selection output table 
#############################
## global 
NF_global<- lmer(CK_ann_cumul_dens ~ scale(LC) +scale(salinity_top_avg) + (1|year) + scale(June_temp_mean) + scale(Apr30_max) + scale(Ln_fry) + 
                           scale(depth_avg) + scale(pathways) +(1|site), 
                         data = df, na.action=na.fail, REML = F) 
summary(NF_global)
vif_table <- data.frame(vif(NF_global)) # all well under 2
write.csv(vif_table, "global model vif table.csv")

## removal of year specific variables (temp and flow) that do not appear influential
NF.temp.flow<- lmer(CK_ann_cumul_dens ~ scale(LC) + scale(salinity_top_avg) + (1|year) + scale(Ln_fry) + 
                              scale(depth_avg) + scale(pathways) +(1|site), 
                            data = df, na.action=na.fail, REML = F) 
summary(NF.temp.flow)

## eliminating site specific variables (salnity and pathway) that do not appear influential
NF.sal.path<- lmer(CK_ann_cumul_dens ~ scale(LC) + (1|year) + scale(Apr30_max) + scale(Ln_fry) + 
                             scale(depth_avg) + scale(June_temp_mean)+(1|site), 
                           data = df, na.action=na.fail, REML = F) 
summary(NF.sal.path)

###### individual removal of rest of non-influential variables
# flow 
NF.temp.sal.path.flow<- lmer(CK_ann_cumul_dens ~ scale(LC) + (1|year)  + scale(Ln_fry) + 
                                       scale(depth_avg) +(1|site), 
                                     data = df, na.action=na.fail, REML = F) 
summary(NF.temp.sal.path.flow)

# null fit
null_model <- lmer(CK_ann_cumul_dens ~ (1|year)+(1|site), 
                   data = df, na.action=na.fail, REML = F)
summary(null_model)

# removing depth
NF_best_fit <- lmer(CK_ann_cumul_dens ~ scale(LC) + (1|year)  + scale(Ln_fry) + 
                      +(1|site), 
                    data = df, na.action=na.fail, REML = F) 
summary(NF_best_fit)


# only fry
only_fry_fit <- lmer(CK_ann_cumul_dens ~ (1|year)  + scale(Ln_fry) + 
                       (1|site), 
                     data = df, na.action=na.fail, REML = F) 
summary(only_fry_fit)


selection_output <- model.sel(NF_global, NF.temp.flow,
                              NF.sal.path,NF_best_fit,
                              NF.temp.sal.path.flow,
                              null_model, only_fry_fit
)

colnames(selection_output)[1] <- "intercept"
write.csv(selection_output, "final_model_select.csv")






##### ####################model from previous

model_LC_fry<- lmer(CK_ann_cumul_dens ~ scale(LC) + scale(Ln_fry)+ (1|year) +
                      (1|site), 
                    data = df, na.action=na.fail, REML = F)
#summary(model_LC_fry)
df$predict.LC.fry <-  predict(model_LC_fry,
                                 data.frame(
                                   year = df$year,
                                   Ln_fry = df$Ln_fry,
                                   depth_avg = df$depth_avg,
                                   LC = df$LC,
                                   site = df$site), type="response")

summary(lm(predict.LC.fry~CK_ann_cumul_dens, data =df)) #0.7564


###########
# residuals and qq plots
plot(model_LC_fry)
plotQQunif(model_LC_fry)



######## best model LC vs predicted
# single fit 
# actual vs predicted
predict_year <- ggplot(df, aes(x=CK_ann_cumul_dens, y= predict.LC.fry))+
  geom_point(aes(color = year))+
  scale_color_gradientn(colours = rainbow(2))+
  theme_bw()+
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
  #theme(legend.position="none")+
  #ggtitle("measured vs predicted annual CK cumulative density\n CK~fry+LC+depth")+
  geom_smooth(method='lm')+
  ylab("Predicted Cumulative CK density ln(fish/ha)")+
  xlab("Observed Cumulative CK density ln(fish/ha)")+
  annotate("text", x = 1, y = 12, label = "B", size =5)+
  expand_limits(x = 0, y = 0) +
  coord_cartesian(expand = FALSE)+
  scale_x_continuous(breaks = seq(0, 13, by = 1), limits = c(0,13))+ 
  scale_y_continuous(breaks = seq(0, 13, by = 1), limits = c(0,13))

predict_site <- ggplot(df, aes(x=CK_ann_cumul_dens, y= predict.LC.fry))+
  geom_point(aes(color = site))+
  theme_bw()+
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
  #theme(legend.position="none")+
  #ggtitle("measured vs predicted annual CK cumulative density\n CK~fry+LC+depth")+
  geom_smooth(method='lm')+
  ylab("Predicted Cumulative CK density ln(fish/ha)")+
  xlab("Observed CK density ln(fish/ha)")+
  annotate("text", x = 1, y = 12, label = "A", size =5)+
  expand_limits(x = 0, y = 0) +
  coord_cartesian(expand = FALSE)+
  scale_x_continuous(breaks = seq(0, 13, by = 1), limits = c(0,13))+ 
  scale_y_continuous(breaks = seq(0, 13, by = 1), limits = c(0,13))

predict_out <- predict_site + predict_year







###############################################################
############### Re-fitting models - one before, one after
df_before <- df[df$year %in% c(2000:2005),]
df_after <-df[df$year %in% c(2016:2021),]

fit_before <- lmer(CK_ann_cumul_dens ~ scale(LC) + (1|year)  + scale(Ln_fry) +
                     (1|site),
                   data = df_before, na.action=na.fail, REML = F)
summary(fit_before)

fit_after <- lmer(CK_ann_cumul_dens ~ scale(LC) + (1|year)  + scale(Ln_fry) +
                    (1|site),
                  data = df_after, na.action=na.fail, REML = F)
summary(fit_after)


#### actual LC measurments

df_before$predict <- predict(fit_before,
                             data.frame(
                               year = df_before$year,
                               Ln_fry = df_before$Ln_fry,
                               LC = df_before$LC,
                               site = df_before$site), type="response")

df_after$predict <- predict(fit_after,
                            data.frame(
                              year = df_after$year,
                              Ln_fry = df_after$Ln_fry,
                              LC = df_after$LC,
                              site = df_after$site), type="response")


df_before$bef_aft <- "before: 2000-2005"
df_after$bef_aft <-  "after: 2016-2021"
# rbind to one dataframe fro plotting
df_two_fits <- rbind(df_before,df_after)

LC_two_fits_bef_aft <- ggplot(df_two_fits, aes(x=LC, y= predict))+
  geom_point(aes(color = bef_aft))+
  theme_bw()+
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
  #theme(legend.position="none")+
  #ggtitle("Preditions from two separate models: before and after\n actual LC vs fitted annual CK cumulative density\n CK~fry+LC+depth")+
  ggtitle("")+
  geom_smooth(aes(group = bef_aft, color =  bef_aft),method='lm')+
  ylab("Predicted Cumulative CK density ln(fish/ha)")+
  xlab("Observed LC")+
  #annotate("text", x = 0.060, y = 5, label ="dens = 28.4906*LC + 7.8721", size = 3, color ="#00BFC4")+
  #annotate("text", x = 0.060, y = 4.5, label ="dens = 29.0569*LC + 7.4941", size = 3, color ="#F8766D")+
  expand_limits(x = 0, y = 0) +
  coord_cartesian(expand = FALSE)+
  scale_y_continuous(breaks = seq(0, 12, by = 1), limits = c(0,12))+
  scale_x_continuous(breaks = seq(0, 0.1, by = 0.02), limits = c(0,0.1))+
  theme(legend.title=element_blank())

# linear fits
summary(lm(predict ~ LC, data = df_two_fits[df_two_fits$bef_aft == "before: 2000-2005",]))
summary(lm(predict ~ LC, data = df_two_fits[df_two_fits$bef_aft == "after: 2016-2021",]))


#### simulated across all possible LC measurements

pos_LC <- seq(0,1, 0.001)

## Prepping before dataframe for graphing
df_before$predict <- NULL
df_before$station_year <- paste(df_before$site, df_before$year, sep="_")
station_year_before <- df_before$station_year
pos_LC_before <- data.frame(crossing(var1 = pos_LC, var2 = station_year_before))
colnames(pos_LC_before) <- c("pos_LC","station_year")

df_before2 <- merge(pos_LC_before, df_before, all=T)

df_before2$predict2 <- predict(fit_before,
                              data.frame(
                                year = df_before2$year,
                                Ln_fry = df_before2$Ln_fry,
                                LC = df_before2$pos_LC,
                                site = df_before2$site), type="response")

#aggregate for means by LC
agg_before <- aggregate(df_before2$predict2, by = list(df_before2$pos_LC, df_before2$bef_aft), FUN = mean)
              colnames(agg_before) <- c("pos_LC", "bef_aft", "mean_CK_dens")

## Prepping after dataframe for graphing
df_after$predict <- NULL
df_after$station_year <- paste(df_after$site, df_before$after, sep="_")
station_year_after <- df_after$station_year
pos_LC_after <- data.frame(crossing(var1 = pos_LC, var2 = station_year_after))
colnames(pos_LC_after) <- c("pos_LC","station_year")

df_after2 <- merge(pos_LC_after, df_after, all=T)

df_after2$predict2 <- predict(fit_after,
                               data.frame(
                                 year = df_after2$year,
                                 Ln_fry = df_after2$Ln_fry,
                                 LC = df_after2$pos_LC,
                                 site = df_after2$site), type="response")

# aggreagate for means by LC
agg_after <- aggregate(df_after2$predict2, by = list(df_after2$pos_LC, df_after2$bef_aft), FUN = mean)
              colnames(agg_after) <- c("pos_LC", "bef_aft", "mean_CK_dens")


# combine before graphging
df_2_mod_combined <- rbind(agg_before, agg_after)



LC_two_fits_bef_aft_short_range <- ggplot(df_2_mod_combined, aes(x=pos_LC, y= mean_CK_dens))+
  geom_line(aes(color = bef_aft), size = 1.5)+
  theme_bw()+
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
  #theme(legend.position="none")+
  #ggtitle("Preditions from two separate models: before and after\npossible LC values vs mean fitted annual CK cumulative density\n CK~fry+LC+depth")+
  ggtitle("")+
  #geom_smooth(aes(group = bef_aft, color =  bef_aft),method='lm')+
  xlim(0, 0.1)+
  ylim(0, 15)+
  ylab("Mean Predicted Cumulative CK density ln(fish/ha)")+
  xlab("Possible LC Values")+
  #annotate("text", x = 0.03, y = 14, label ="dens = 25.53*LC + 7.870", size = 3, color ="#00BFC4")+
  #annotate("text", x = 0.03, y = 13, label ="dens = 36.55*LC + 7.171", size = 3, color ="#F8766D")+  expand_limits(x = 0, y = 0) +
  expand_limits(x = 0, y = 0) +
  coord_cartesian(expand = FALSE)+
  theme(legend.title=element_blank())

# linear fits
summary(lm(mean_CK_dens~pos_LC, data = df_2_mod_combined[df_2_mod_combined$bef_aft == "before: 2000-2005",]))
summary(lm(mean_CK_dens~pos_LC, data = df_2_mod_combined[df_2_mod_combined$bef_aft == "after: 2016-2021",]))

LC_out1 <- LC_two_fits_bef_aft / LC_two_fits_bef_aft_short_range
########################################################



###################
# Summary Tables
###################
sdmean <- rep(c("sd","mean"), length(unique(df$site)))

## means and sd by site
ck_mean <- aggregate(df$CK_ann_cumul_dens, by = list(df$site), FUN = mean)
ck_sd <- aggregate(df$CK_ann_cumul_dens, by = list(df$site), FUN = sd)
ck_site <- merge(ck_sd,ck_mean, all =T)  
colnames(ck_site)<- c("site", "Ann_ck_dens")
rownames(ck_site)<- paste(ck_site$site, sdmean, sep="_")


sal_mean <- aggregate(df$salinity_top_avg, by = list(df$site), FUN = mean)
sal_sd <-aggregate(df$salinity_top_avg, by = list(df$site), FUN = sd)
sal_site <- merge(sal_sd,sal_mean, all =T)  
colnames(sal_site)<- c("site", "salinity")
rownames(sal_site)<- paste(sal_site$site, sdmean, sep="_")
  
depth_mean <- aggregate(df$depth_avg, by = list(df$site), FUN = mean)
depth_sd<- aggregate(df$depth_avg, by = list(df$site), FUN = sd)
depth_site <- merge(depth_sd,depth_mean, all =T)  
colnames(depth_site)<- c("site", "depth")
rownames(depth_site)<- paste(depth_site$site, sdmean, sep="_")
  
temp_mean<-aggregate(df$June_temp_mean, by = list(df$site), FUN = mean)
temp_sd<- aggregate(df$June_temp_mean, by = list(df$site), FUN = sd)
temp_site <- merge(temp_sd,temp_mean, all =T)  
colnames(temp_site)<- c("site", "temp")
rownames(temp_site)<- paste(temp_site$site, sdmean, sep="_")
  
LC_mean <- aggregate(df$LC, by = list(df$site), FUN = mean)
LC_sd <- aggregate(df$LC, by = list(df$site), FUN = sd)
LC_site <- merge(LC_sd,LC_mean, all =T)  
colnames(LC_site)<- c("site", "LC")
rownames(LC_site)<- paste(LC_site$site, sdmean, sep="_")
  
flow_mean <- aggregate(df$Apr30_max, by = list(df$site), FUN = mean)
flow_sd <- aggregate(df$Apr30_max, by = list(df$site), FUN = sd)
flow_site <- merge(flow_sd,flow_mean, all =T)  
colnames(flow_site)<- c("site", "flow")
rownames(flow_site)<- paste(flow_site$site, sdmean, sep="_")

path_mean <- aggregate(df$pathways, by = list(df$site), FUN = mean)
path_sd <- aggregate(df$pathways, by = list(df$site), FUN = sd)
path_site <- merge(path_sd,path_mean, all =T)  
colnames(path_site)<- c("site", "pathways")
rownames(path_site)<- paste(path_site$site, sdmean, sep="_")

fry_mean <- aggregate(df$Ln_fry, by = list(df$site), FUN = mean)
fry_sd <- aggregate(df$Ln_fry, by = list(df$site), FUN = sd)
fry_site <- merge(fry_sd,fry_mean, all =T)  
colnames(fry_site)<- c("site", "fry")
rownames(fry_site)<- paste(fry_site$site, sdmean, sep="_")

site_table <- do.call(cbind, list(ck_site,sal_site,depth_site,temp_site, LC_site, fry_site, flow_site, path_site))
site_table <- site_table[,c("site","Ann_ck_dens","fry", "LC", "salinity","depth","temp","flow", "pathways")]
site_table <- cbind(site_table,sdmean)
site_table <- site_table[order(site_table$site, site_table$sdmean),]
rownames(site_table)<- NULL
site_table <- site_table[,c("site","sdmean", "Ann_ck_dens", "LC", "salinity", "depth", 
                            "pathways")] # temp flow and fry were not site specific parameters


write.csv(site_table,"Parameter summary by site.csv", row.names = F)


## means and sd by year
sdmean2 <- rep(c("sd","mean"), length(unique(df$year)))

## means and sd by site
ck_mean_y <- aggregate(df$CK_ann_cumul_dens, by = list(df$year), FUN = mean)
ck_sd_y <- aggregate(df$CK_ann_cumul_dens, by = list(df$year), FUN = sd)
ck_year <- merge(ck_sd_y,ck_mean_y, all =T)  
colnames(ck_year)<- c("year", "Ann_ck_dens")
rownames(ck_year)<- paste(ck_year$year, sdmean2, sep="_")


sal_mean_y <- aggregate(df$salinity_top_avg, by = list(df$year), FUN = mean)
sal_sd_y <-aggregate(df$salinity_top_avg, by = list(df$year), FUN = sd)
sal_year <- merge(sal_sd_y,sal_mean_y, all =T)  
colnames(sal_year)<- c("year", "salinity")
rownames(sal_year)<- paste(sal_year$year, sdmean2, sep="_")

depth_mean_y <- aggregate(df$depth_avg, by = list(df$year), FUN = mean)
depth_sd_y<- aggregate(df$depth_avg, by = list(df$year), FUN = sd)
depth_year <- merge(depth_sd_y,depth_mean_y, all =T)  
colnames(depth_year)<- c("year", "depth")
rownames(depth_year)<- paste(depth_year$year, sdmean2, sep="_")

temp_mean_y<-aggregate(df$June_temp_mean, by = list(df$year), FUN = mean)
temp_sd_y<- aggregate(df$June_temp_mean, by = list(df$year), FUN = sd)
temp_year <- merge(temp_sd_y,temp_mean_y, all =T)  
colnames(temp_year)<- c("year", "temp")
rownames(temp_year)<- paste(temp_year$year, sdmean2, sep="_")

LC_mean_y <- aggregate(df$LC, by = list(df$year), FUN = mean)
LC_sd_y <- aggregate(df$LC, by = list(df$year), FUN = sd)
LC_year <- merge(LC_sd_y,LC_mean_y, all =T)  
colnames(LC_year)<- c("year", "LC")
rownames(LC_year)<- paste(LC_year$year, sdmean2, sep="_")

flow_mean_y <- aggregate(df$Apr30_max, by = list(df$year), FUN = mean)
flow_sd_y <- aggregate(df$Apr30_max, by = list(df$year), FUN = sd)
flow_year <- merge(flow_sd_y,flow_mean_y, all =T)  
colnames(flow_year)<- c("year", "flow")
rownames(flow_year)<- paste(flow_year$year, sdmean2, sep="_")

path_mean_y <- aggregate(df$pathways, by = list(df$year), FUN = mean)
path_sd_y <- aggregate(df$pathways, by = list(df$year), FUN = sd)
path_year <- merge(path_sd_y,path_mean_y, all =T)  
colnames(path_year)<- c("year", "pathways")
rownames(path_year)<- paste(path_year$year, sdmean2, sep="_")

fry_mean_y <- aggregate(df$Ln_fry, by = list(df$year), FUN = mean)
fry_sd_y <- aggregate(df$Ln_fry, by = list(df$year), FUN = sd)
fry_year <- merge(fry_sd_y,fry_mean_y, all =T)  
colnames(fry_year)<- c("year", "fry")
rownames(fry_year)<- paste(fry_year$year, sdmean2, sep="_")


year_table <- do.call(cbind, list(ck_year,fry_year, sal_year,depth_year,temp_year, LC_year, flow_year, path_year))
year_table <- year_table[,c("year","Ann_ck_dens","fry","LC", "salinity","depth","temp","flow", "pathways")]
year_table <- cbind(year_table,sdmean2)
year_table <- year_table[order(year_table$year, year_table$sdmean2),]
rownames(year_table)<- NULL
year_table <- year_table[,c("year","sdmean2", "Ann_ck_dens","fry", "LC", "salinity", "depth", "temp", "flow", 
                            "pathways")]

write.csv(year_table,"Parameter summary by year.csv", row.names = F)



######################





############################
## final correlation plots for appendix
df_cor_plots <- df[,c("CK_ann_cumul_dens","LC", "pathways", "salinity_top_avg","depth_avg","Ln_fry", "Apr30_max",  
                          "June_temp_mean"
                          )]
colnames(df_cor_plots) <- c("Chinook\nDensity","LC", "Pathways", "Salinity","Depth","Fry\nOutmigrants", "River\nDischarge",  
                            "Temperature")
ggpairs(df_cor_plots)

df_cor_plots2 <- df_cor_plots[,c("LC", "Pathways", "Salinity","Depth","Fry\nOutmigrants", "River\nDischarge",  
                                 "Temperature")]
ggpairs(df_cor_plots2, upper = list(continuous = wrap(ggally_cor, stars = FALSE)))



######################