---
title: "Untitled"
output: html_document
date: '2022-09-14'
---

---
title: "NPN_Gams_Analysis"

output: html_notebook
---

```{r}

knitr::opts_chunk$set(echo = TRUE)
rm(list = ls())
library(data.table) # Parallel, much faster version of data.frames
#library(mgcViz)
library(plyr)
library(MCMCpack)
library(quantreg)
library(DataCombine)
#library(tidyverse)
set.seed(2)


```


```{r}
#raw_data <- fread("E:/Dropbox/Macro_2021/Papers/Paper1_assessments/Data/SimHerb_Data_v2/RQ_Sim_Distribution_slope4_skewDistribution.csv")

raw_data <- fread("E:/Dropbox/Macro_2021/Papers/Paper1_assessments/Data/SimHerb_Data_v2/RQ_Sim_Distribution_slope4_skewDistribution.csv")

raw_data$Latitude <- raw_data$Y
raw_data$Longitude <- raw_data$X
raw_data$First_Yes_Year <- raw_data$year
print(nrow(raw_data))
#raw_data <- raw_data[raw_data$genus_species == "Syringa vulgaris"]
```

```{r}

#formula_locked <- paste("Local_sample ~ Ann_tmean_13")

#DOY_column = "Local_sample"
#year_column = 'First_Yes_Year'
#taxon_column = "genus_species"


#iter_fit <- rq(as.formula(formula_locked),
#                         data = raw_data, tau = 0.1)

#iter_fit$coefficients

```




```{r}
Count_dt <-  function(raw_data,
                      group_column = "genus_species", 
                      threshold = 10) { 
  
  #param raw_data = input data table
  #param group_column = colum name to be grouped on (genus_species)
  #param threshold = minimum number of samples per species
  
print(group_column)
count_data <- count(raw_data,  group_column)#get number of rows per species
count_data$sample_count <- count_data$freq
setDT(count_data)
count_data[, c("freq") := NULL]   #remove N column

print(nrow(raw_data))
print('sneps')
raw_data <- merge(raw_data, count_data, by = c(group_column), no.dups = TRUE) 
print('fleps')

raw_data <- raw_data[raw_data$sample_count >= threshold] #filter out species with <1000 specimens
print(nrow(raw_data))

return(raw_data)

}

```


```{r}


  
qreg_analysis <-  function(raw_data, #data table for input data
                           predict_data, #data table to use for predictions (usually the same as raw_data)
                           out_path, #string path for where to de[posit output files]
                           monthly_params, #vector of montly parameter types to evaluate
                           DOY_column, #string of the name of the DOY column
                           year_column, #string of the year column name
                           taxon_column, #string of the taxon column name
                           sample_count_min = 50,#minimum number of samples in a species to evaluate it
                           quantiles = c(0.1, 0.5, 0.9)) 
                            {  
  
  predict_data_trimmed <- FALSE
  
  raw_data <- Count_dt(raw_data, taxon_column, sample_count_min)

  
  rq_eval <- data.table(genus_species=character(), parameter=character(), month=numeric(),
                       quantiles = numeric(), AIC = numeric(), Slope = numeric())
  species_count = 0
  
  for (iter_species in unique(raw_data$genus_species)){
    species_count <- species_count + 1
    try({
      
    
    
    
    print(nrow(raw_data[raw_data[, get(taxon_column)] == iter_species]))
    
    
    if (nrow(raw_data[raw_data[, get(taxon_column)] == iter_species]) >= sample_count_min){
      iter_data <- raw_data[raw_data[,get(taxon_column)] == iter_species]
      predict_data_raw_iter <-  raw_data[raw_data[,get(taxon_column)] == iter_species]
      predict_data_iter <-  predict_data[predict_data[,get(taxon_column)] == iter_species]
    
      for (iter_quantile in quantiles){
      print(paste('QUANTILE:   ', iter_quantile))
       prev_params <- '' # set previously selected parameters to none
       
    
      for (param in monthly_params){
        rq_eval_species <- rq_eval[rq_eval[, get(taxon_column)] == iter_species & rq_eval[, get("quantiles")] == iter_quantile]
        
        formula_locked <- paste(DOY_column, " ~ ")
        
        
        
        if (nrow(rq_eval_species)>0) { #if previous parameters have been selected for this species, add them to formula   
                                          #before testing additional parameters
        
          for (i in nrow(rq_eval_species)) {
            
              if (rq_eval_species[i,month]<10){
                selected_month <- paste('0', as.character(rq_eval_species[i,month]), sep = '')
                  
              }
              else{
                selected_month <- as.character(rq_eval_species[i,month])
              }
                
              formula_locked = paste(formula_locked, rq_eval_species[i,parameter], selected_month, ' + ', sep = '')
            
                
          }   
        }
      
       #create blank data.table to document selectwed parameters  
       rq_eval_monthly <- data.table(genus_species=character(), parameter=character(), month=numeric(), 
                                     quantiles = numeric(),
                                      AIC = numeric(), Slope = numeric()) 
       
       
       #iterate over months for each parameter
       for (iter_month in 13:13){
        if (iter_month <10){
          iter_month_str <- paste("0", as.character(iter_month), sep = '')
             }
        else  {
          iter_month_str <- as.character(iter_month)
          }  
        iter_param <- paste( param, iter_month_str, sep = '')
        print(iter_param)
        iter_fit <- rq(as.formula(paste(formula_locked, as.character(iter_param), sep = '')),
                         data = iter_data, tau = iter_quantile)
        print(paste("AIC:",as.character(AIC(iter_fit))))
        
        #calculate MSE
        fit_summ <-summary(iter_fit)
        #print(paste("MSE:",iter_MSE))
        #print(plot(iter_fit))
        #print(summary(iter_fit))
        #print(gam.check(iter_fit))
        
        iter_AIC <- AIC.rq(iter_fit)
        
        
        if (dim(summary(iter_fit)$coefficients)[1] == 2){
          iter_Slope <- iter_fit$coefficients[2] 
          }
        
        else{
          iter_Slope <- 0.0
          }
        
        
        #print(iter_species)
        #print(param)
        #print(iter_month)
        #print(iter_dev_expl)
        #print(iter_AIC)
        #print(plot(iter_fit), pages = 1)
        
        
        iter_eval <- data.table(genus_species=c(iter_species), parameter=c(param), month=c(iter_month), 
                                quantiles = c(iter_quantile), AIC = c(iter_AIC), Slope = iter_Slope)
        
        rq_eval_monthly <- rbind(rq_eval_monthly, iter_eval)
      }
        
      
      #identify month for a given parameter that minimizes AIC
      rq_eval_grouped <- rq_eval_monthly[rq_eval_monthly[,.(AIC == min(AIC)), by = list(genus_species, quantiles)]$V1 == TRUE]
      
      rq_eval <- rbind(rq_eval, rq_eval_grouped)  #add that row to the permanent selected parameters file
       
       
       } 
      
    
    
    #Once model selected, fit actual model
    
     
      #write formula from selected params  
      rq_eval_iter <- rq_eval[rq_eval[, get(taxon_column)] == iter_species & rq_eval[, get("quantiles")] == iter_quantile]
    
    
    #create formula
    formula_final <- paste(DOY_column, " ~ ")
    
    print("start)")
    for (i in 1:nrow(rq_eval_iter)) {
      #print(i)
            if (rq_eval_iter[i,month]<10){
              selected_month <- paste('0', as.character(rq_eval_iter[i,month]), sep = '')
                
            }
            else{
              selected_month <- as.character(rq_eval_iter[i,month])
            }
        
            if (i < nrow(rq_eval_iter)){
              formula_final <- paste(formula_final, rq_eval_iter[i,parameter], selected_month, '+ ', sep = '')   
              
            }
            else{
              formula_final <- paste(formula_final, rq_eval_iter[i,parameter], selected_month, sep = '')
              print('final_formula')
              print(formula_final)
            }
    
      #ensurre that rows with NA data are removed from data used for predictions
      predict_data_iter <- predict_data_iter[!is.na(predict_data_iter[[paste(rq_eval_iter[i,parameter], selected_month, sep = '')]]),]
      predict_data_raw_iter <- predict_data_raw_iter[!is.na(predict_data_raw_iter[[paste(rq_eval_iter[i,parameter], selected_month, sep = '')]]),]
      
      } 
    
    
    print(paste("predict length: ", nrow(predict_data_iter)))
    
    
    
    1
    final_fit <- rq(as.formula(as.character(formula_final)), data = iter_data, tau = iter_quantile)                       
              
    
    if (predict_data_trimmed == FALSE){
      
      #isolate key columns to trim fat for data to be output with predictions
        #predict_data_trimmed <- predict_data_iter[, list(unique_ID, genus_species, slope, param, sigma, duration, skew, Local_indiv_mean_DOY, Local_indiv_onset_DOY, Local_indiv_termination_DOY, Latitude, Longitude, First_Yes_Year)]
        
        predict_data_trimmed <- predict_data_iter[, list(unique_ID, genus_species, slope, param, sigma, duration, Local_indiv_mean_DOY, Local_indiv_onset_DOY, Local_indiv_termination_DOY, Latitude, Longitude, First_Yes_Year)]
        
        predict_data_raw_trimmed <- predict_data_raw_iter[, list(unique_ID, genus_species, slope, param, sigma, duration, skew, Local_10pct_mean_DOY, Local_50pct_mean_DOY, Local_90pct_mean_DOY,  Local_10pct_onset_DOY, Local_90pct_mean_DOY, Local_90pct_termination_DOY,  Local_indiv_mean_DOY, Local_indiv_onset_DOY, Local_indiv_termination_DOY, Latitude, Longitude, First_Yes_Year, get(DOY_column))] 
        
        #predict_data_raw_trimmed <- predict_data_raw_iter[, list(unique_ID, genus_species, slope, param, sigma, duration, Local_10pct_mean_DOY, Local_50pct_mean_DOY, Local_90pct_mean_DOY,  Local_10pct_onset_DOY, Local_90pct_mean_DOY, Local_90pct_termination_DOY,  Local_indiv_mean_DOY, Local_indiv_onset_DOY, Local_indiv_termination_DOY, Latitude, Longitude, First_Yes_Year, get(DOY_column))] 
        
        setnames(predict_data_raw_trimmed, "V20", DOY_column)
    }
    
      
      
      colname_test <- paste("RQ_Pred", as.character(iter_quantile), as.character(DOY_column), sep = "_")
      print(colname_test)
      
            print(paste("predict length2: ", nrow(predict_data_iter)))
            print(paste("predict length trimmed: ", nrow(predict_data_trimmed)))
            
            
            print(paste('PREDICTED_LEN: ', nrow(predict(final_fit, predict_data_iter, type = "none", interval = c("none")))))
            #print(predict(final_fit, predict_data_iter, type = "none", interval = c("none")))
            #return(predict_data_iter)
            #return(predict(final_fit, predict_data_iter, type = "none", interval = c("none")))
      
            predict_data_trimmed[[paste("RQ_Pred", as.character(iter_quantile), as.character(DOY_column), sep = "_")]]  <- predict(final_fit, predict_data_iter, type = "none", interval = c("none"))
            predict_data_raw_trimmed[[paste("RQ_Pred", as.character(iter_quantile), as.character(DOY_column), sep = "_")]]  <- predict(final_fit, predict_data_raw_iter, type = "none", interval = c("none"))
            #predict_data_iter[[paste("GAM_Pred", as.character(quantile), as.character(DOY_column), sep = "_")]] <-
            
            print(paste("RQ_Pred", as.character(iter_quantile), as.character(DOY_column), sep = "_"))
      
      }#end of Quantiles
      
    
    if(exists("predictions_concat")== FALSE){
      print('ZAKDORN')
      predictions_concat <- predict_data_trimmed
      predictions_concat_raw <- predict_data_raw_trimmed
      }
    else{
      print('FELPNAR')
      predictions_concat <- rbindlist(list(predictions_concat, predict_data_trimmed))
      predictions_concat_raw <- rbindlist(list(predictions_concat_raw, predict_data_raw_trimmed))
      
    }
    predict_data_trimmed <- FALSE
    predict_data_raw_trimmed <- FALSE

    
    }
  })     
  }
  #####end of species iteration
  
  
#once all models selected, write selected params to csv
write.csv(rq_eval, paste(out_path, "RQ_SelectedParams_Simulated_Distribution_SlopeComp_", DOY_column, ".csv", sep = ''))
write.csv(predictions_concat, paste(out_path, "RQ_Predict_Test_rawData_Simulated_Distribution_SlopeComp_", DOY_column, ".csv", sep = ''))
write.csv(predictions_concat_raw, paste(out_path, "RQ_Predict_Test_simData_Simulated_Distribution_SlopeComp_", DOY_column, ".csv", sep = ''))

print(species_count)
print(nrow(predictions_concat))
print(nrow(predictions_concat_raw))



return(predictions_concat)  
}






```


```{r}


sim_data_1000 <- raw_data[,.SD[sample(.N, min(1000,.N))],by = genus_species]

q <- qreg_analysis(sim_data_1000, sim_data_1000, out_path= "E:/Dropbox/Macro_2021/Papers/Paper1_assessments/Data/Simulated_Distribution_Models/RQ_SimHerb_1000_Simulated_Distribution_Slope4_30sigma_normal_skew_test20_",
                   monthly_params = c("Ann_tmean_"),
                   DOY_column = "Local_sample",
                   year_column = 'First_Yes_Year',
                   taxon_column = "genus_species",
                   quantile = c(0.1, 0.5, 0.9))

```



