# SCRIPT NR. 9

gam.variable.select <- function(y,a,z, rownames, taxon.name, family, corr.plot = F, criterion = "significance", cross.criterion = FALSE, simp.threshold = 0.01, silent = T) {
  
  #### variable definition and dependencies ####
  
  #y = dataframe with input data
  
  #a = explanatory columns
  
  #z = response column
  
  #rownames = character list with the names of the explanatory variables in the correct order corresponding to explanatory columns
  
  #taxon.name = name of the species (needed to name the output data frames)
  
  library(mgcv)                                                                                                   #load mgcv package for gerneralized addattive models
  
  library(car)
  
  x = 0                                                                                                           #count variable to name the models from 1 to length (a)
  
  n = length(a)                                                                                                   #needed to set proggress bar + end point for loop
  
  Models <- c()                                                                                                   #create list called models (currently not needed)
  
  Names <- c()                                                                                                    #create list called names (currently not needed)
  
  NullDeviance <- c()                                                                                             #create list to write in each null deviance
  
  Deviance <- c()                                                                                                 #create list to write in each deviance
  
  Data <- c()                                                                                                     #create dataframe to write all results in
  
  
  
  message(paste("####----#####----####----#### Variable selection process ####----#####----####----####"))
  
  message("\n")
  
  
  
  pb <- txtProgressBar(min = 0, max = n, style = 3)                                                               #progress bar, style = 3 (text based)
  
  
  
  #### calculate single variable pseudo r^2 #### 
  
  for (i in a){
    
    x <- x + 1
    
    gammodel <- gam(y[,z] ~ s(y[,i], k=5), family = family, se.fit = T, method = "REML")
    
    Models <- append(Models,paste0("gam_bio", x))
    
    Names <- append(Names,paste0("bio", x))
    
    ov <- deviance(gammodel) / df.residual(gammodel)
    
    NullDeviance <- append(NullDeviance, gammodel$null.deviance)
    
    Deviance<- append(Deviance, gammodel$deviance)
    
    Sys.sleep(0.001) 
    
    if (silent == T) {
      
      setTxtProgressBar(pb, x)
      
      gc()
      
    }
    
    if (silent == F) {
      
      message("\n")
      
      message(paste0("gam_bio", x,"  Overdispersion: ", ov))
      
      if (abs(ov) > 1) {
        
        message("Your model is overdispersed! Check your model family again.")
        
        message("\n")
        
      }
      
      else {
        
        message("Your model is not overdispersed! Congratulations!")
        
        message("\n")
        
      }
      
      setTxtProgressBar(pb, x)
      
      gc()
      
    }
    
  }
  
  
  
  Data <- as.data.frame(cbind(NullDeviance,Deviance))
  
  Data$r2 <- 1-(Data$Deviance/Data$NullDeviance)
  
  Data$r2 <- Data$r2*100
  
  rownames(Data) <- rownames
  
  #assign(paste0("gam_all_var_",taxon.name), Data, envir = .GlobalEnv)
  
  
  
  
  
  #### correlation plot #### 
  
  
  
  if (corr.plot == T) {
    
    library(corrplot)
    
    cor_matrix <- cor(y[c(a)], use = "pairwise.complete.obs")
    
    cor_matrix <- abs(cor_matrix)
    
    plot <- corrplot(cor_matrix, 
                     
                     method = "circle",
                     
                     tl.cex = 0.8, 
                     
                     p.mat = cor_matrix, 
                     
                     sig.level = 0.7)
    
  }
  
  
  
  #### select variables based on highest r^2 and correlation coefficient #### 
  
  vgl_df <- Data
  
  b <- a
  
  final_var_df <- c()
  
  final_var_df <- as.data.frame(final_var_df)
  
  l <- c(1:n)
  
  
  
  options(warn=-1)
  
  for (i in c(1:nrow(Data))) {
    
    var_selected_value <- max(Data$r2)
    
    var_selected_num <- which(Data$r2 == var_selected_value)
    
    col_num_selected_var <- b[var_selected_num]
    
    final_var_df <- rbind(final_var_df, Data[var_selected_num,])
    
    
    
    list_corr_coef <- c()
    
    
    
    for (v in b) {
      
      cor_coef <- stats::cor(y[,col_num_selected_var],y[,v], use = "na.or.complete")
      
      list_corr_coef <- append(list_corr_coef, paste0(abs(cor_coef))) 
      
    }
    
    
    
    remaining_var_count <- which(list_corr_coef < 0.7)
    
    as.numeric(remaining_var_count)
    
    l <- remaining_var_count
    
    
    
    Data <- Data[remaining_var_count,]
    
    b <- b[remaining_var_count]
    
  }
  
  options(warn=0)
  
  assign(paste0("gam_selected_var_",taxon.name), final_var_df, envir = .GlobalEnv)
  
  message("\n")
  
  message(paste("####----#####----####----#### Global Envirnoment Message ####----#####----####----#### "))
  
  message("\n")
  
  message(paste0("'","gam_all_var_",taxon.name,"'", " Full Variable list was sent to your Global Evironment"))
  
  message(paste0("'","gam_selected_var_",taxon.name,"'", " Final Variable list was sent to your Global Evironment"))
  
  
  
  
  
  r2_selected <- intersect(vgl_df$r2, final_var_df$r2)
  
  model_building_index <- which(vgl_df$r2 %in% r2_selected)
  
  c <- a[model_building_index]
  
  colnames <- colnames(y)
  
  
  
  
  
  
  
  #### build final model #### 
  
  gam.formula <- paste0(colnames[z], " ~ ",paste0(" s(",colnames[c],", k =3) ",collapse="+"))
  
  gam.formula <- as.formula(gam.formula)
  
  final.gam.model <- gam(formula = gam.formula,data = y, family = family, se.fit = T, method = "REML")
  
  message("\n")
  
  message(paste("####----#####----####----#### Model stats ####----#####----####----#### "))
  
  message("\n")
  
  message(paste0("Pseudo R^2: ",1-(final.gam.model$deviance/final.gam.model$null.deviance)))
  
  #assign(paste0("final_gam_model",taxon.name), final.gam.model, envir = .GlobalEnv)
  
  message("\n")
  
  
  
  
  
  
  
  #### model simplification ####
  
  message(paste("####----#####----####----#### Model simplification process ####----#####----####----#### "))
  
  message("\n")
  
  message(paste0("Simplification process starts: "))
  
  message("\n")
  
  
  
  if (criterion  == "significance" | cross.criterion == T) {
    
    message(paste0("Simplification criterion: significance"))
    
    pb_2 <- txtProgressBar(min = 0, max = 1, style = 3)
    
    model_building_index_simplification <- c  
    
    gam.running.model <- final.gam.model
    
    running_r_2 <- 1-(final.gam.model$deviance/final.gam.model$null.deviance)
    
    simp_r_2 <- 1-(final.gam.model$deviance/final.gam.model$null.deviance)
    
    simplification_df <- final_var_df
    
    
    
    anova.gam.model <- anova.gam(gam.running.model)
    
    simplification_df$significance <- anova.gam.model$s.pv
    
    simplification_df$significance <- as.numeric(simplification_df$significance)
    
    
    
    while (running_r_2 == simp_r_2 && sum(simplification_df$significance) > 0) {
      
      anova.gam.model <- anova.gam(gam.running.model)
      
      simplification_df$significance <- anova.gam.model$s.pv
      
      simplification_df$significance <- as.numeric(simplification_df$significance)
      
      simp_count <- which(simplification_df$significance == max(simplification_df$significance))
      
      simp_index_significance <- model_building_index_simplification
      
      model_building_index_simplification <- model_building_index_simplification[-simp_count]
      
      
      
      if(length(model_building_index_simplification) > 0) {
        
        
        
        #### build simplified model #### 
        
        gam.formula.simp <- paste0(colnames[z], " ~ ",paste0(" s(",colnames[model_building_index_simplification],", k =3) ",collapse="+"))
        
        gam.formula.simp <- as.formula(gam.formula.simp)
        
        gam.simp.model <- gam(formula = gam.formula.simp,data = y, family = family, se.fit = T, method = "REML")
        
        simp_r_2 <- 1-(gam.simp.model$deviance/gam.simp.model$null.deviance)
        
        
        
        if (running_r_2-simp_r_2 <= simp.threshold && running_r_2 != simp_r_2) {
          
          running_r_2 <- simp_r_2
          
          gam.running.model <- gam.simp.model
          
          simplification_df <- simplification_df[-simp_count,]
          
        } else {
          
          running_r_2 <- running_r_2
          
          gam.running.model <- gam.running.model
          
        }
        
        
        
      } else {
        
        running_r_2 <- 1
        
        simp_r_2 <- 0
        
      }
      
    }
    
    
    
    # check if there was any signigicance simplification. if not write unsimplified formula in gam.formula.simp
    
    if(exists("gam.formula.simp") == F) {
      
      gam.formula.simp <- gam.formula
      
      simp_index_significance <- c
      
      
      
    }
    
    
    
    gam.formula.simp.significance <- gam.formula.simp
    
    
    
    setTxtProgressBar(pb_2, 1)
    
    message("\n")
    
    message(paste0("Simplification criterion: significance (finished)"))
    
    message("\n")
    
  }
  
  
  
  if (criterion == "square" | cross.criterion == T) {
    
    message(paste0("Simplification criterion: square"))
    
    model_building_index_simplification <- c  
    
    gam.running.model <- final.gam.model
    
    running_r_2 <- 1-(final.gam.model$deviance/final.gam.model$null.deviance)
    
    simp_r_2 <- 1-(final.gam.model$deviance/final.gam.model$null.deviance)
    
    simplification_df <- final_var_df
    
    
    
    
    
    simp_condition <- 1
    
    
    
    while_count <- 0
    
    while (simp_condition == 1) {
      
      while_count <- while_count + 1
      
      simp_square <- c()
      
      message(paste0("Run ", while_count))
      
      pb_3 <- txtProgressBar(min = 0, max = length(model_building_index_simplification), style = 3) 
      
      if (length(model_building_index_simplification) > 1) {                                                     ### adjusted
        
        for (i in c(1:length(model_building_index_simplification))) {
          
          model_building_index_simplification_i <- model_building_index_simplification[-i]
          
          gam.formula.simp <- paste0(colnames[z], " ~ ",paste0(" s(",colnames[model_building_index_simplification_i],", k =3) ",collapse="+"))
          
          gam.formula.simp <- as.formula(gam.formula.simp)
          
          gam.simp.model <- gam(formula = gam.formula.simp,data = y, family = family, se.fit = T, method = "REML")
          
          simp_r_2 <- 1-(gam.simp.model$deviance/gam.simp.model$null.deviance)
          
          simp_square <- append(simp_square, simp_r_2)
          
          setTxtProgressBar(pb_3, i)
          
        }
        
        if (running_r_2-max(simp_square) <= simp.threshold) {
          
          discard_val <- colnames[model_building_index_simplification[which(simp_square == max(simp_square))]]
          
          model_building_index_simplification <- model_building_index_simplification[-which(simp_square == max(simp_square))]
          
          simp_condition <- 1
          
          if (silent == F) {
            
            message("\n")
            
            message(paste0("Discarded: ", discard_val))
            
          }
          
        } else {
          
          simp_condition <- 0
          
          gam.formula.simp.square <- paste0(colnames[z], " ~ ",paste0(" s(",colnames[model_building_index_simplification],", k =3) ",collapse="+"))
          
          gam.formula.simp.square <- as.formula(gam.formula.simp.square)
          
          simp_index_square <- model_building_index_simplification
          
        }
        
      } else {                                                                        ### adjusted
        
        simp_condition <- 0 
        
        simp_index_square <- model_building_index_simplification    
        
        gam.formula.simp.square <- paste0(colnames[z], " ~ ",paste0(" s(",colnames[model_building_index_simplification],", k =3) ",collapse="+"))
        
        gam.formula.simp.square <- as.formula(gam.formula.simp.square)                ### adjusted
        
      }                                                                               ### adjusted
      
      message("\n")
      
    }
    
    message(paste0("Simplification criterion: significance (finished)"))
    
    message("\n")
    
  }
  
  
  
  ###### adjusted
  
  assign(paste0("var_list"), colnames[simp_index_square], envir = .GlobalEnv)
  
  ######
  
  
  
  #### write output of model simplification #### 
  
  message(paste("####----#####----####----#### Global Envirnoment Message ####----#####----####----#### "))
  
  message("\n")
  
  
  
  if (criterion == "significance" | cross.criterion == T) {
    
    gam.final.simp.model.significance <- gam(formula = gam.formula.simp.significance, data = y, family = family, se.fit = T, method = "REML")
    
    #assign(paste0("final_gam_model_simplified_significance_",taxon.name), gam.final.simp.model.significance, envir = .GlobalEnv)
    
    message(paste0("'","final_gam_model_simplified_significance_",taxon.name,"'", " Simplified model was sent to your Global Evironment"))
    
    message("\n")
    
    
    
    if (cross.criterion == F) {
      
      message(paste("####----#####----####----#### Model simplification stats ####----#####----####----#### "))
      
      message("\n")
      
      message(paste0("Pseudo R^2: ",1-(gam.final.simp.model.significance$deviance/gam.final.simp.model.significance$null.deviance)))
      
      message(paste("simplified variables siginificance : "))
      
      message(paste(" ",colnames[simp_index_significance], " "))
      
    }
    
  }
  
  
  
  if (criterion == "square" | cross.criterion == T) {
    
    gam.final.simp.model.square <- gam(formula = gam.formula.simp.square, data = y, family = family, se.fit = T, method = "REML")
    
    #assign(paste0("final_gam_model_simplified_square_",taxon.name), gam.final.simp.model.square, envir = .GlobalEnv)
    
    message(paste0("'","final_gam_model_simplified_square_",taxon.name,"'", " Simplified model was sent to your Global Evironment"))
    
    message("\n")
    
    
    
    if (cross.criterion == F) {
      
      message(paste("####----#####----####----#### Model simplification stats ####----#####----####----#### "))
      
      message("\n")
      
      message(paste0("Pseudo R^2: ",1-(gam.final.simp.model.square$deviance/gam.final.simp.model.square$null.deviance)))
      
      message(paste("simplified variables quare         : "))
      
      message(paste(" ",colnames[simp_index_square], " "))
      
    }
    
  }
  
  
  
  if (cross.criterion == T){
    
    message(paste("####----#####----####----#### Model simplification stats ####----#####----####----#### "))
    
    message("\n")
    
    message(paste0("Pseudo R^2: ",1-(gam.final.simp.model.significance$deviance/gam.final.simp.model.significance$null.deviance)))
    
    message(paste("simplified variables siginificance : "))
    
    message(paste(" ",colnames[simp_index_significance], " "))
    
    message("\n")
    
    message(paste0("Pseudo R^2: ",1-(gam.final.simp.model.square$deviance/gam.final.simp.model.square$null.deviance)))
    
    message(paste("simplified variables quare         : "))
    
    message(paste(" ",colnames[simp_index_square], " "))
    
  }
  
}