Select_BestMod_Zbins_ObsAnalysis_PDSI_higherPeakSWE.R

#load in libraries
library(splines)
library(foreach)
library(nlme)
library(mgcv)
library(Metrics)
library(MASS)
library(sm) 
library(mgcv.helper)
library(tibble)
library(dplyr)
library(fmsb)
library(randomForest)
library(KRLS)
library(matrixStats)
#remove variables from prior run
rm(list = ls())
#define model concurvity threshold:
concurvity_threshold = 0.4

##=========================pick best model based on BurnArea_GAM_v0*.R
#read in data

#0-1100m:
Climate_Fire_Data = read.table("/Volumes/Pruina_External_Elements/DroughtFireSnow/Data/AnlysisData/Annual_Fire_Climate_Obs_forecast_Ziter1_higherPeakSWE.csv",sep=",")
#define variables:
WYs = Climate_Fire_Data[1:37,1]
BAs = Climate_Fire_Data[1:37,2]
BAs_Z1 = (BAs*0.000247105) /(10^6);
SWEIs_Z1 = Climate_Fire_Data[1:37,3]
WinPRCP_Z1 = Climate_Fire_Data[1:37,4]
WinTMP_Z1 = Climate_Fire_Data[1:37,5]
WinterPDSI_Z1 = Climate_Fire_Data[1:37,6]
Spring_PRCP_Z1 = Climate_Fire_Data[1:37,7] 
Spring_TMP_Z1 = Climate_Fire_Data[1:37,8] 
Spring_PDSI_Z1 = Climate_Fire_Data[1:37,9] 
Spring_SWEI_Z1 = Climate_Fire_Data[1:37,10] 
Spring_VPD_Z1 = Climate_Fire_Data[1:37,11] 
Winter_VPD_Z1 = Climate_Fire_Data[1:37,12] 
MODIS_BA_Z1 = Climate_Fire_Data[1:37,13] 
Spring_ET_Z1 = Climate_Fire_Data[1:37,14] 
Spring_PET_Z1 = Climate_Fire_Data[1:37,15] 
Spring_PETminusET_Z1 = Climate_Fire_Data[1:37,16] 
Winter_ET_Z1 = Climate_Fire_Data[1:37,17] 
Winter_PET_Z1 = Climate_Fire_Data[1:37,18] 
Winter_PETminusET_Z1 = Climate_Fire_Data[1:37,19] 
Winter_DroughtArea_Z1 = Climate_Fire_Data[1:37,31] 
Spring_DroughtArea_Z1 = Climate_Fire_Data[1:37,32] 

Z1_ID = rep(1,37)

#relate MODIS and MTBS BA linearly to allow MODIS to gap fill MTBS 2020 BA:
BA_MTBS <- BAs_Z1[18:36]
BA_MODIS <- MODIS_BA_Z1[18:36]
BA_DF = data.frame(BA_MTBS=BA_MTBS,BA_MODIS=BA_MODIS)
BA_lm <- lm(data=BA_DF,BA_MTBS~BA_MODIS)
data_2020 = data.frame(BA_MODIS = MODIS_BA_Z1[37])
BA_2020_MTBS=predict(BA_lm,newdata=data_2020)
BAs_Z1[37] = BA_2020_MTBS
if (BAs_Z1[37] < 0){
  BAs_Z1[37] = 0
}

#1100-2200:
Climate_Fire_Data = read.table("/Volumes/Pruina_External_Elements/DroughtFireSnow/Data/AnlysisData/Annual_Fire_Climate_Obs_forecast_Ziter2_higherPeakSWE.csv",sep=",")
#define variables:
#define variables:
WYs = Climate_Fire_Data[1:37,1]
BAs = Climate_Fire_Data[1:37,2]
BAs_Z2 = (BAs*0.000247105) /(10^6);
SWEIs_Z2 = Climate_Fire_Data[1:37,3]
WinPRCP_Z2 = Climate_Fire_Data[1:37,4]
WinTMP_Z2 = Climate_Fire_Data[1:37,5]
WinterPDSI_Z2 = Climate_Fire_Data[1:37,6]
Spring_PRCP_Z2 = Climate_Fire_Data[1:37,7] 
Spring_TMP_Z2 = Climate_Fire_Data[1:37,8] 
Spring_PDSI_Z2 = Climate_Fire_Data[1:37,9] 
Spring_SWEI_Z2 = Climate_Fire_Data[1:37,10] 
Spring_VPD_Z2 = Climate_Fire_Data[1:37,11] 
Winter_VPD_Z2 = Climate_Fire_Data[1:37,12] 
MODIS_BA_Z2 = Climate_Fire_Data[1:37,13] 
Spring_ET_Z2 = Climate_Fire_Data[1:37,14] 
Spring_PET_Z2 = Climate_Fire_Data[1:37,15] 
Spring_PETminusET_Z2 = Climate_Fire_Data[1:37,16] 
Winter_ET_Z2 = Climate_Fire_Data[1:37,17] 
Winter_PET_Z2 = Climate_Fire_Data[1:37,18] 
Winter_PETminusET_Z2 = Climate_Fire_Data[1:37,19] 
Winter_DroughtArea_Z2 = Climate_Fire_Data[1:37,31] 
Spring_DroughtArea_Z2 = Climate_Fire_Data[1:37,32]  
Z2_ID = rep(2,37)

#relate MODIS and MTBS BA linearly to allow MODIS to gap fill MTBS 2020 BA:
BA_MTBS <- BAs_Z2[18:36]
BA_MODIS <- MODIS_BA_Z2[18:36]
BA_DF = data.frame(BA_MTBS=BA_MTBS,BA_MODIS=BA_MODIS)
BA_lm <- lm(data=BA_DF,BA_MTBS~BA_MODIS)
data_2020 = data.frame(BA_MODIS = MODIS_BA_Z2[37])
BA_2020_MTBS=predict(BA_lm,newdata=data_2020)
BAs_Z2[37] = BA_2020_MTBS
if (BAs_Z2[37] < 0){
  BAs_Z2[37] = 0
}

#2200-3300:
Climate_Fire_Data = read.table("/Volumes/Pruina_External_Elements/DroughtFireSnow/Data/AnlysisData/Annual_Fire_Climate_Obs_forecast_Ziter3_higherPeakSWE.csv",sep=",")
#define variables:
#define variables:
WYs = Climate_Fire_Data[1:37,1]
BAs = Climate_Fire_Data[1:37,2]
BAs_Z3 = (BAs*0.000247105) /(10^6);
SWEIs_Z3 = Climate_Fire_Data[1:37,3]
WinPRCP_Z3 = Climate_Fire_Data[1:37,4]
WinTMP_Z3 = Climate_Fire_Data[1:37,5]
WinterPDSI_Z3 = Climate_Fire_Data[1:37,6]
Spring_PRCP_Z3 = Climate_Fire_Data[1:37,7] 
Spring_TMP_Z3 = Climate_Fire_Data[1:37,8] 
Spring_PDSI_Z3 = Climate_Fire_Data[1:37,9] 
Spring_SWEI_Z3 = Climate_Fire_Data[1:37,10] 
Spring_VPD_Z3 = Climate_Fire_Data[1:37,11] 
Winter_VPD_Z3 = Climate_Fire_Data[1:37,12] 
MODIS_BA_Z3 = Climate_Fire_Data[1:37,13] 
Spring_ET_Z3 = Climate_Fire_Data[1:37,14] 
Spring_PET_Z3 = Climate_Fire_Data[1:37,15] 
Spring_PETminusET_Z3 = Climate_Fire_Data[1:37,16] 
Winter_ET_Z3 = Climate_Fire_Data[1:37,17] 
Winter_PET_Z3 = Climate_Fire_Data[1:37,18] 
Winter_PETminusET_Z3 = Climate_Fire_Data[1:37,19] 
Winter_DroughtArea_Z3 = Climate_Fire_Data[1:37,31] 
Spring_DroughtArea_Z3 = Climate_Fire_Data[1:37,32]  
Z3_ID = rep(3,37)

#relate MODIS and MTBS BA linearly to allow MODIS to gap fill MTBS 2020 BA:
BA_MTBS <- BAs_Z3[18:36]
BA_MODIS <- MODIS_BA_Z3[18:36]
BA_DF = data.frame(BA_MTBS=BA_MTBS,BA_MODIS=BA_MODIS)
BA_lm <- lm(data=BA_DF,BA_MTBS~BA_MODIS)
data_2020 = data.frame(BA_MODIS = MODIS_BA_Z3[37])
BA_2020_MTBS=predict(BA_lm,newdata=data_2020)
BAs_Z3[37] = BA_2020_MTBS
if (BAs_Z3[37] < 0){
  BAs_Z3[37] = 0
}

idx<-c(1:37)
WYs = rep(WYs[idx],3)
BAs = cbind(t(BAs_Z1[idx]),t(BAs_Z2[idx]),t(BAs_Z3[idx]))
SWEIs = cbind(t(SWEIs_Z1[idx]),t(SWEIs_Z2[idx]),t(SWEIs_Z3[idx]))
WinPRCP = cbind(t(WinPRCP_Z1[idx]),t(WinPRCP_Z2[idx]),t(WinPRCP_Z3[idx]))
WinTMP = cbind(t(WinTMP_Z1[idx]),t(WinTMP_Z2[idx]),t(WinTMP_Z3[idx]))
WinterPDSI = cbind(t(WinterPDSI_Z1[idx]),t(WinterPDSI_Z2[idx]),t(WinterPDSI_Z3[idx]))
Spring_PRCP = cbind(t(Spring_PRCP_Z1[idx]),t(Spring_PRCP_Z2[idx]),t(Spring_PRCP_Z3[idx]))
Spring_TMP = cbind(t(Spring_TMP_Z1[idx]),t(Spring_TMP_Z2[idx]),t(Spring_TMP_Z3[idx]))
Spring_PDSI = cbind(t(Spring_PDSI_Z1[idx]),t(Spring_PDSI_Z2[idx]),t(Spring_PDSI_Z3[idx]))
Spring_SWEI = cbind(t(Spring_SWEI_Z1[idx]),t(Spring_SWEI_Z2[idx]),t(Spring_SWEI_Z3[idx]))
Spring_VPD = cbind(t(Spring_VPD_Z1[idx]),t(Spring_VPD_Z2[idx]),t(Spring_VPD_Z3[idx]))
Winter_VPD = cbind(t(Winter_VPD_Z1[idx]),t(Winter_VPD_Z2[idx]),t(Winter_VPD_Z3[idx]))
Spring_ET = cbind(t(Spring_ET_Z1[idx]),t(Spring_ET_Z2[idx]),t(Spring_ET_Z3[idx]))
Spring_PET = cbind(t(Spring_PET_Z1[idx]),t(Spring_PET_Z2[idx]),t(Spring_PET_Z3[idx]))
Spring_PETminusET = cbind(t(Spring_PETminusET_Z1[idx]),t(Spring_PETminusET_Z2[idx]),t(Spring_PETminusET_Z3[idx]))
Winter_ET = cbind(t(Winter_ET_Z1[idx]),t(Winter_ET_Z2[idx]),t(Winter_ET_Z3[idx]))
Winter_PET = cbind(t(Winter_PET_Z1[idx]),t(Winter_PET_Z2[idx]),t(Winter_PET_Z3[idx]))
Winter_PETminusET = cbind(t(Winter_PETminusET_Z1[idx]),t(Winter_PETminusET_Z2[idx]),t(Winter_PETminusET_Z3[idx]))
Winter_DroughtArea = cbind(t(Winter_DroughtArea_Z1[idx]),t(Winter_DroughtArea_Z2[idx]),t(Winter_DroughtArea_Z3[idx]))
Spring_DroughtArea = cbind(t(Spring_DroughtArea_Z1[idx]),t(Spring_DroughtArea_Z2[idx]),t(Spring_DroughtArea_Z3[idx]))
Winter_Spring_Temp = (Spring_TMP+WinTMP)/2
Winter_Spring_Precip = (Spring_PRCP+WinPRCP)/2
Winter_Spring_VPD = (Winter_VPD+Spring_VPD)/2
Winter_Spring_PET = (Winter_PET+Spring_PET)/2
Winter_Spring_ET = (Winter_ET+Spring_ET)/2

Zs = cbind(t(Z1_ID[idx]),t(Z2_ID[idx]),t(Z3_ID[idx]))

Climate_Fire_DF = data.frame(Spring_PDSI=t(Spring_PDSI),WinPRCP = t(WinPRCP), WinTMP = t(WinTMP),Spring_PRCP = t(Spring_PRCP),Spring_TMP=t(Spring_TMP),Spring_VPD=t(Spring_VPD),Winter_VPD=t(Winter_VPD),Spring_ET=t(Spring_ET),Spring_PET=t(Spring_PET),Winter_ET=t(Winter_ET),Winter_PET=t(Winter_PET),Spring_DroughtArea=t(Spring_DroughtArea),Winter_Spring_Temp=t(Winter_Spring_Temp),Winter_Spring_Precip=t(Winter_Spring_Precip),Winter_Spring_VPD=t(Winter_Spring_VPD),Winter_Spring_PET=t(Winter_Spring_PET),Winter_Spring_ET=t(Winter_Spring_ET))

#loop through all combinations of covariates:
Y = as.vector(BAs)
nvars <- 17


#for  combos of 2 variables:
# n_predictors <- 2
# ncombos = factorial(nvars)/(factorial(n_predictors)*factorial(nvars-n_predictors))
# x=1:nvars
# combo_IDs<-combn(x,n_predictors)
# 
# store_AIC_2vars = c()
# store_R_2vars = c()
# store_i_2vars=c()
# for (i in 1:ncombos){
#   col1 <- combo_IDs[1,i]
#   col2 <- combo_IDs[2,i]
#   
#   predictor1 <-Climate_Fire_DF[,col1]
#   predictor2 <-Climate_Fire_DF[,col2]
#   #PCA combo model:
#   Covariates_DF = data.frame(predictor1=predictor1,predictor2=predictor2)
#   df.pca <- prcomp(Covariates_DF, center = TRUE,scale. = TRUE)
#   pcs = df.pca$x
#   pc1 = pcs[,1]
#   pc2 = pcs[,2]
#   current_DF_pcs <- data.frame(Y=Y,predictor1=pc1,predictor2=pc2,predictor3=as.vector(Zs))
#   mod <- gam(data=current_DF_pcs,Y~s(predictor1)+s(predictor2)+predictor3,family="gaussian")
#   
#   X = concurvity(mod,full=FALSE)
#   X=X$estimate
#   idx <- which(X==1)
#   X[idx] = 0
#   Max_Concuvrity <- max(X)
#   if (Max_Concuvrity <= concurvity_threshold){
#     AIC = mod$aic
#     R = cor(predict(mod),Y)
#     
#     store_R_2vars = c(store_R_2vars,R)
#     store_AIC_2vars = c(store_AIC_2vars,AIC)
#     store_i_2vars = c(store_i_2vars,i)
#   }
# }
# 
# idx_bestmod_2vars <- which(store_AIC_2vars==min(store_AIC_2vars))
# store_R_2vars[idx_bestmod_2vars]
# idx_bestmod_2vars = store_i_2vars[idx_bestmod_2vars]
# 
# #for  combos of 3 variables:
# n_predictors <- 3
# ncombos = factorial(nvars)/(factorial(n_predictors)*factorial(nvars-n_predictors))
# x=1:nvars
# combo_IDs<-combn(x,n_predictors)
# 
# store_AIC_3vars = c()
# store_R_3vars = c()
# store_i_3vars = c()
# for (i in 1:ncombos){
#   col1 <- combo_IDs[1,i]
#   col2 <- combo_IDs[2,i]
#   col3 <- combo_IDs[3,i]
#   
#   predictor1 <-Climate_Fire_DF[,col1]
#   predictor2 <-Climate_Fire_DF[,col2]
#   predictor3 <-Climate_Fire_DF[,col3]
#   #PCA combo model:
#   Covariates_DF = data.frame(predictor1=predictor1,predictor2=predictor2,predictor3=predictor3)
#   df.pca <- prcomp(Covariates_DF, center = TRUE,scale. = TRUE)
#   pcs = df.pca$x
#   pc1 = pcs[,1]
#   pc2 = pcs[,2]
#   pc3 = pcs[,3]
#   current_DF_pcs <- data.frame(Y=Y,predictor1=pc1,predictor2=pc2,predictor3=pc3,predictor4=as.vector(Zs))
#   mod <- gam(data=current_DF_pcs,Y~s(predictor1)+s(predictor2)+s(predictor3)+predictor4,family="gaussian")
#   
#   X = concurvity(mod,full=FALSE)
#   X=X$estimate
#   idx <- which(X==1)
#   X[idx] = 0
#   Max_Concuvrity <- max(X)
#   if (Max_Concuvrity <= concurvity_threshold){
#     AIC = mod$aic
#     R = cor(predict(mod),Y)
#     
#     store_R_3vars = c(store_R_3vars,R)
#     store_AIC_3vars = c(store_AIC_3vars,AIC)
#     store_i_3vars = c(store_i_3vars,i)
#   }
# }
# 
# idx_bestmod_3vars <- which(store_AIC_3vars==min(store_AIC_3vars))
# store_R_3vars[idx_bestmod_3vars]
# idx_bestmod_3vars = store_i_3vars[idx_bestmod_3vars]
# 
# 
# #for  combos of 4 variables:
# n_predictors <- 4
# ncombos = factorial(nvars)/(factorial(n_predictors)*factorial(nvars-n_predictors))
# x=1:nvars
# combo_IDs<-combn(x,n_predictors)
# 
# store_AIC_4vars = c()
# store_R_4vars = c()
# store_i_4vars=c()
# for (i in 1:ncombos){
#   col1 <- combo_IDs[1,i]
#   col2 <- combo_IDs[2,i]
#   col3 <- combo_IDs[3,i]
#   col4 <- combo_IDs[4,i]
#   
#   predictor1 <-Climate_Fire_DF[,col1]
#   predictor2 <-Climate_Fire_DF[,col2]
#   predictor3 <-Climate_Fire_DF[,col3]
#   predictor4 <-Climate_Fire_DF[,col4]
#   #PCA combo model:
#   Covariates_DF = data.frame(predictor1=predictor1,predictor2=predictor2,predictor3=predictor3,predictor4=predictor4)
#   df.pca <- prcomp(Covariates_DF, center = TRUE,scale. = TRUE)
#   pcs = df.pca$x
#   pc1 = pcs[,1]
#   pc2 = pcs[,2]
#   pc3 = pcs[,3]
#   pc4 = pcs[,4]
#   current_DF_pcs <- data.frame(Y=Y,predictor1=pc1,predictor2=pc2,predictor3=pc3,predictor4=pc4,predictor5=as.vector(Zs))
#   mod <- gam(data=current_DF_pcs,Y~s(predictor1)+s(predictor2)+s(predictor3)+s(predictor4)+predictor5,family="gaussian")
#   
#   X = concurvity(mod,full=FALSE)
#   X=X$estimate
#   idx <- which(X==1)
#   X[idx] = 0
#   Max_Concuvrity <- max(X)
#   if (Max_Concuvrity <= concurvity_threshold){
#     AIC = mod$aic
#     R = cor(predict(mod),Y)
#     
#     store_R_4vars = c(store_R_4vars,R)
#     store_AIC_4vars = c(store_AIC_4vars,AIC)
#     store_i_4vars = c(store_i_4vars,i)
#   }
# }
# 
# idx_bestmod_4vars <- which(store_AIC_4vars==min(store_AIC_4vars))
# store_R_4vars[idx_bestmod_4vars]
# idx_bestmod_4vars = store_i_3vars[idx_bestmod_4vars]
# 
# 
# #for  combos of 5 variables:
# n_predictors <- 5
# ncombos = factorial(nvars)/(factorial(n_predictors)*factorial(nvars-n_predictors))
# x=1:nvars
# combo_IDs<-combn(x,n_predictors)
# 
# store_AIC_5vars = c()
# store_R_5vars = c()
# store_i_5vars=c()
# for (i in 1:ncombos){
#   col1 <- combo_IDs[1,i]
#   col2 <- combo_IDs[2,i]
#   col3 <- combo_IDs[3,i]
#   col4 <- combo_IDs[4,i]
#   col5 <- combo_IDs[5,i]
# 
#   predictor1 <-Climate_Fire_DF[,col1]
#   predictor2 <-Climate_Fire_DF[,col2]
#   predictor3 <-Climate_Fire_DF[,col3]
#   predictor4 <-Climate_Fire_DF[,col4]
#   predictor5 <-Climate_Fire_DF[,col5]
#   #PCA combo model:
#   Covariates_DF = data.frame(predictor1=predictor1,predictor2=predictor2,predictor3=predictor3,predictor4=predictor4,predictor5=predictor5)
#   df.pca <- prcomp(Covariates_DF, center = TRUE,scale. = TRUE)
#   pcs = df.pca$x
#   pc1 = pcs[,1]
#   pc2 = pcs[,2]
#   pc3 = pcs[,3]
#   pc4 = pcs[,4]
#   pc5 = pcs[,5]
#   current_DF_pcs <- data.frame(Y=Y,predictor1=pc1,predictor2=pc2,predictor3=pc3,predictor4=pc4,predictor5=pc5,predictor6=as.vector(Zs))
#   mod <- gam(data=current_DF_pcs,Y~s(predictor1)+s(predictor2)+s(predictor3)+s(predictor4)+s(predictor5)+predictor6,family="gaussian")
# 
#   X = concurvity(mod,full=FALSE)
#   X=X$estimate
#   idx <- which(X==1)
#   X[idx] = 0
#   Max_Concuvrity <- max(X)
#   if (Max_Concuvrity <= concurvity_threshold){
#     AIC = mod$aic
#     R = cor(predict(mod),Y)
# 
#     store_R_5vars = c(store_R_5vars,R)
#     store_AIC_5vars = c(store_AIC_5vars,AIC)
#     store_i_5vars = c(store_i_5vars,i)
#   }
# }
# 
# idx_bestmod_5vars <- which(store_AIC_5vars==min(store_AIC_5vars))
# store_R_5vars[idx_bestmod_5vars]
# idx_bestmod_5vars = store_i_5vars[idx_bestmod_5vars]

# #for  combos of 6 variables:
# n_predictors <- 6
# ncombos = factorial(nvars)/(factorial(n_predictors)*factorial(nvars-n_predictors))
# x=1:nvars
# combo_IDs<-combn(x,n_predictors)
# 
# store_AIC_6vars = c()
# store_R_6vars = c()
# store_i_6vars=c()
# for (i in 1:ncombos){
#   col1 <- combo_IDs[1,i]
#   col2 <- combo_IDs[2,i]
#   col3 <- combo_IDs[3,i]
#   col4 <- combo_IDs[4,i]
#   col5 <- combo_IDs[5,i]
#   col6 <- combo_IDs[6,i]
#   
#   predictor1 <-Climate_Fire_DF[,col1]
#   predictor2 <-Climate_Fire_DF[,col2]
#   predictor3 <-Climate_Fire_DF[,col3]
#   predictor4 <-Climate_Fire_DF[,col4]
#   predictor5 <-Climate_Fire_DF[,col5]
#   predictor6 <-Climate_Fire_DF[,col6]
#   
#   #PCA combo model:
#   Covariates_DF = data.frame(predictor1=predictor1,predictor2=predictor2,predictor3=predictor3,predictor4=predictor4,predictor5=predictor5,predictor6=predictor6)
#   df.pca <- prcomp(Covariates_DF, center = TRUE,scale. = TRUE)
#   pcs = df.pca$x
#   pc1 = pcs[,1]
#   pc2 = pcs[,2]
#   pc3 = pcs[,3]
#   pc4 = pcs[,4]
#   pc5 = pcs[,5]
#   pc6 = pcs[,6]
#   current_DF_pcs <- data.frame(Y=Y,predictor1=pc1,predictor2=pc2,predictor3=pc3,predictor4=pc4,predictor5=pc5,predictor6=pc6,predictor7=as.vector(Zs))
#   mod <- gam(data=current_DF_pcs,Y~s(predictor1)+s(predictor2)+s(predictor3)+s(predictor4)+s(predictor5)+s(predictor6)+predictor7,family="gaussian")
#   
#   X = concurvity(mod,full=FALSE)
#   X=X$estimate
#   idx <- which(X==1)
#   X[idx] = 0
#   Max_Concuvrity <- max(X)
#   if (Max_Concuvrity <= concurvity_threshold){
#     AIC = mod$aic
#     R = cor(predict(mod),Y)
#     
#     store_R_6vars = c(store_R_6vars,R)
#     store_AIC_6vars = c(store_AIC_6vars,AIC)
#     store_i_6vars = c(store_i_6vars,i)
#   }
# }
# 
# idx_bestmod_6vars <- which(store_AIC_6vars==min(store_AIC_6vars))
# store_R_6vars[idx_bestmod_6vars]
# idx_bestmod_6vars = store_i_6vars[idx_bestmod_6vars]

#for  combos of 7 variables:
n_predictors <- 7
ncombos = factorial(nvars)/(factorial(n_predictors)*factorial(nvars-n_predictors))
x=1:nvars
combo_IDs<-combn(x,n_predictors)

store_AIC_7vars = c()
store_R_7vars = c()
store_i_7vars = c()
for (i in 1:ncombos){
  col1 <- combo_IDs[1,i]
  col2 <- combo_IDs[2,i]
  col3 <- combo_IDs[3,i]
  col4 <- combo_IDs[4,i]
  col5 <- combo_IDs[5,i]
  col6 <- combo_IDs[6,i]
  col7 <- combo_IDs[7,i]

  predictor1 <-Climate_Fire_DF[,col1]
  predictor2 <-Climate_Fire_DF[,col2]
  predictor3 <-Climate_Fire_DF[,col3]
  predictor4 <-Climate_Fire_DF[,col4]
  predictor5 <-Climate_Fire_DF[,col5]
  predictor6 <-Climate_Fire_DF[,col6]
  predictor7 <-Climate_Fire_DF[,col7]
  #PCA combo model:
  Covariates_DF = data.frame(predictor1=predictor1,predictor2=predictor2,predictor3=predictor3,predictor4=predictor4,predictor5=predictor5,predictor6=predictor6,predictor7=predictor7)
  df.pca <- prcomp(Covariates_DF, center = TRUE,scale. = TRUE)
  pcs = df.pca$x
  pc1 = pcs[,1]
  pc2 = pcs[,2]
  pc3 = pcs[,3]
  pc4 = pcs[,4]
  pc5 = pcs[,5]
  pc6 = pcs[,6]
  pc7 = pcs[,7]
  current_DF_pcs <- data.frame(Y=Y,predictor1=pc1,predictor2=pc2,predictor3=pc3,predictor4=pc4,predictor5=pc5,predictor6=pc6,predictor7=pc7,predictor8=as.vector(Zs))
  mod <- gam(data=current_DF_pcs,Y~s(predictor1)+s(predictor2)+s(predictor3)+s(predictor4)+s(predictor5)+s(predictor6)+s(predictor7)+predictor8,family="gaussian")

  X = concurvity(mod,full=FALSE)
  X=X$estimate
  idx <- which(X==1)
  X[idx] = 0
  Max_Concuvrity <- max(X)
  if (Max_Concuvrity <= concurvity_threshold){
    AIC = mod$aic
    R = cor(predict(mod),Y)

    store_R_7vars = c(store_R_7vars,R)
    store_AIC_7vars = c(store_AIC_7vars,AIC)
    store_i_7vars = c(store_i_7vars,i)
  }
}

idx_bestmod_7vars <- which(store_AIC_7vars==min(store_AIC_7vars))
store_R_7vars[idx_bestmod_7vars]
idx_bestmod_7vars = store_i_7vars[idx_bestmod_7vars]


# #for  combos of 8 variables:
# n_predictors <- 8
# ncombos = factorial(nvars)/(factorial(n_predictors)*factorial(nvars-n_predictors))
# x=1:nvars
# combo_IDs<-combn(x,n_predictors)
# 
# store_AIC_8vars = c()
# store_R_8vars = c()
# store_i_8vars = c()
# for (i in 1:ncombos){
#   col1 <- combo_IDs[1,i]
#   col2 <- combo_IDs[2,i]
#   col3 <- combo_IDs[3,i]
#   col4 <- combo_IDs[4,i]
#   col5 <- combo_IDs[5,i]
#   col6 <- combo_IDs[6,i]
#   col7 <- combo_IDs[7,i]
#   col8 <- combo_IDs[8,i]
# 
#   predictor1 <-Climate_Fire_DF[,col1]
#   predictor2 <-Climate_Fire_DF[,col2]
#   predictor3 <-Climate_Fire_DF[,col3]
#   predictor4 <-Climate_Fire_DF[,col4]
#   predictor5 <-Climate_Fire_DF[,col5]
#   predictor6 <-Climate_Fire_DF[,col6]
#   predictor7 <-Climate_Fire_DF[,col7]
#   predictor8 <-Climate_Fire_DF[,col8]
#   #PCA combo model:
#   Covariates_DF = data.frame(predictor1=predictor1,predictor2=predictor2,predictor3=predictor3,predictor4=predictor4,predictor5=predictor5,predictor6=predictor6,predictor7=predictor7,predictor8=predictor8)
#   df.pca <- prcomp(Covariates_DF, center = TRUE,scale. = TRUE)
#   pcs = df.pca$x
#   pc1 = pcs[,1]
#   pc2 = pcs[,2]
#   pc3 = pcs[,3]
#   pc4 = pcs[,4]
#   pc5 = pcs[,5]
#   pc6 = pcs[,6]
#   pc7 = pcs[,7]
#   pc8 = pcs[,8]
#   current_DF_pcs <- data.frame(Y=Y,predictor1=pc1,predictor2=pc2,predictor3=pc3,predictor4=pc4,predictor5=pc5,predictor6=pc6,predictor7=pc7,predictor8=pc8,predictor9=as.vector(Zs))
#   mod <- gam(data=current_DF_pcs,Y~s(predictor1)+s(predictor2)+s(predictor3)+s(predictor4)+s(predictor5)+s(predictor6)+s(predictor7)+s(predictor8)+predictor9,family="gaussian")
# 
#   X = concurvity(mod,full=FALSE)
#   X=X$estimate
#   idx <- which(X==1)
#   X[idx] = 0
#   Max_Concuvrity <- max(X)
#   if (Max_Concuvrity <= concurvity_threshold){
#     AIC = mod$aic
#     R = cor(predict(mod),Y)
# 
#     store_R_8vars = c(store_R_8vars,R)
#     store_AIC_8vars = c(store_AIC_8vars,AIC)
#     store_i_8vars = c(store_i_8vars,i)
#   }
# }
# 
# idx_bestmod_8vars <- which(store_AIC_8vars==min(store_AIC_8vars))
# store_R_8vars[idx_bestmod_8vars]
# idx_bestmod_8vars = store_i_8vars[idx_bestmod_8vars]
# 
# #for  combos of 9 variables:
# n_predictors <- 9
# ncombos = factorial(nvars)/(factorial(n_predictors)*factorial(nvars-n_predictors))
# x=1:nvars
# combo_IDs<-combn(x,n_predictors)
# 
# store_AIC_9vars = c()
# store_R_9vars = c()
# store_i_9vars = c()
# for (i in 1:ncombos){
#   col1 <- combo_IDs[1,i]
#   col2 <- combo_IDs[2,i]
#   col3 <- combo_IDs[3,i]
#   col4 <- combo_IDs[4,i]
#   col5 <- combo_IDs[5,i]
#   col6 <- combo_IDs[6,i]
#   col7 <- combo_IDs[7,i]
#   col8 <- combo_IDs[8,i]
#   col9 <- combo_IDs[9,i]
# 
#   predictor1 <-Climate_Fire_DF[,col1]
#   predictor2 <-Climate_Fire_DF[,col2]
#   predictor3 <-Climate_Fire_DF[,col3]
#   predictor4 <-Climate_Fire_DF[,col4]
#   predictor5 <-Climate_Fire_DF[,col5]
#   predictor6 <-Climate_Fire_DF[,col6]
#   predictor7 <-Climate_Fire_DF[,col7]
#   predictor8 <-Climate_Fire_DF[,col8]
#   predictor9 <-Climate_Fire_DF[,col9]
# 
#   #PCA combo model:
#   Covariates_DF = data.frame(predictor1=predictor1,predictor2=predictor2,predictor3=predictor3,predictor4=predictor4,predictor5=predictor5,predictor6=predictor6,predictor7=predictor7,predictor8=predictor8,predictor9=predictor9)
#   df.pca <- prcomp(Covariates_DF, center = TRUE,scale. = TRUE)
#   pcs = df.pca$x
#   pc1 = pcs[,1]
#   pc2 = pcs[,2]
#   pc3 = pcs[,3]
#   pc4 = pcs[,4]
#   pc5 = pcs[,5]
#   pc6 = pcs[,6]
#   pc7 = pcs[,7]
#   pc8 = pcs[,8]
#   pc9 = pcs[,9]
#   current_DF_pcs <- data.frame(Y=Y,predictor1=pc1,predictor2=pc2,predictor3=pc3,predictor4=pc4,predictor5=pc5,predictor6=pc6,predictor7=pc7,predictor8=pc8,predictor9=pc9,predictor10=as.vector(Zs))
#   mod <- gam(data=current_DF_pcs,Y~s(predictor1)+s(predictor2)+s(predictor3)+s(predictor4)+s(predictor5)+s(predictor6)+s(predictor7)+s(predictor8)+s(predictor9)+predictor10,family="gaussian")
# 
#   X = concurvity(mod,full=FALSE)
#   X=X$estimate
#   idx <- which(X==1)
#   X[idx] = 0
#   Max_Concuvrity <- max(X)
#   if (Max_Concuvrity <= concurvity_threshold){
#     AIC = mod$aic
#     R = cor(predict(mod),Y)
# 
#     store_R_9vars = c(store_R_9vars,R)
#     store_AIC_9vars = c(store_AIC_9vars,AIC)
#     store_i_9vars = c(store_i_9vars,i)
#   }
# }
# 
# idx_bestmod_9vars <- which(store_AIC_9vars==min(store_AIC_9vars))
# store_R_9vars[idx_bestmod_9vars]
# idx_bestmod_9vars = store_i_9vars[idx_bestmod_9vars]
# 
# #for  combos of 10 variables:
# n_predictors <- 10
# ncombos = factorial(nvars)/(factorial(n_predictors)*factorial(nvars-n_predictors))
# x=1:nvars
# combo_IDs<-combn(x,n_predictors)
# 
# store_AIC_10vars = c()
# store_R_10vars = c()
# store_i_10vars = c()
# for (i in 1:ncombos){
#   col1 <- combo_IDs[1,i]
#   col2 <- combo_IDs[2,i]
#   col3 <- combo_IDs[3,i]
#   col4 <- combo_IDs[4,i]
#   col5 <- combo_IDs[5,i]
#   col6 <- combo_IDs[6,i]
#   col7 <- combo_IDs[7,i]
#   col8 <- combo_IDs[8,i]
#   col9 <- combo_IDs[9,i]
#   col10 <- combo_IDs[10,i]
# 
#   predictor1 <-Climate_Fire_DF[,col1]
#   predictor2 <-Climate_Fire_DF[,col2]
#   predictor3 <-Climate_Fire_DF[,col3]
#   predictor4 <-Climate_Fire_DF[,col4]
#   predictor5 <-Climate_Fire_DF[,col5]
#   predictor6 <-Climate_Fire_DF[,col6]
#   predictor7 <-Climate_Fire_DF[,col7]
#   predictor8 <-Climate_Fire_DF[,col8]
#   predictor9 <-Climate_Fire_DF[,col9]
#   predictor10 <-Climate_Fire_DF[,col10]
# 
#   #PCA combo model:
#   Covariates_DF = data.frame(predictor1=predictor1,predictor2=predictor2,predictor3=predictor3,predictor4=predictor4,predictor5=predictor5,predictor6=predictor6,predictor7=predictor7,predictor8=predictor8,predictor9=predictor9,predictor10=predictor10)
#   df.pca <- prcomp(Covariates_DF, center = TRUE,scale. = TRUE)
#   pcs = df.pca$x
#   pc1 = pcs[,1]
#   pc2 = pcs[,2]
#   pc3 = pcs[,3]
#   pc4 = pcs[,4]
#   pc5 = pcs[,5]
#   pc6 = pcs[,6]
#   pc7 = pcs[,7]
#   pc8 = pcs[,8]
#   pc9 = pcs[,9]
#   pc10 = pcs[,10]
#   current_DF_pcs <- data.frame(Y=Y,predictor1=pc1,predictor2=pc2,predictor3=pc3,predictor4=pc4,predictor5=pc5,predictor6=pc6,predictor7=pc7,predictor8=pc8,predictor9=pc9,predictor10=pc10,predictor11=as.vector(Zs))
#   mod <- gam(data=current_DF_pcs,Y~s(predictor1)+s(predictor2)+s(predictor3)+s(predictor4)+s(predictor5)+s(predictor6)+s(predictor7)+s(predictor8)+s(predictor9)+s(predictor10)+predictor11,family="gaussian")
# 
#   X = concurvity(mod,full=FALSE)
#   X=X$estimate
#   idx <- which(X==1)
#   X[idx] = 0
#   Max_Concuvrity <- max(X)
#   if (Max_Concuvrity <= concurvity_threshold){
#     AIC = mod$aic
#     R = cor(predict(mod),Y)
# 
#     store_R_10vars = c(store_R_10vars,R)
#     store_AIC_10vars = c(store_AIC_10vars,AIC)
#     store_i_10vars = c(store_i_10vars,i)
#   }
# }
# 
# idx_bestmod_10vars <- which(store_AIC_10vars==min(store_AIC_9vars))
# store_R_10vars[idx_bestmod_10vars]
# idx_bestmod_10vars = store_i_10vars[idx_bestmod_10vars]
# 
# #for  combos of 11 variables:
# n_predictors <- 11
# ncombos = factorial(nvars)/(factorial(n_predictors)*factorial(nvars-n_predictors))
# x=1:nvars
# combo_IDs<-combn(x,n_predictors)
# 
# store_AIC_11vars = c()
# store_R_11vars = c()
# store_i_11vars = c()
# for (i in 1:ncombos){
#   col1 <- combo_IDs[1,i]
#   col2 <- combo_IDs[2,i]
#   col3 <- combo_IDs[3,i]
#   col4 <- combo_IDs[4,i]
#   col5 <- combo_IDs[5,i]
#   col6 <- combo_IDs[6,i]
#   col7 <- combo_IDs[7,i]
#   col8 <- combo_IDs[8,i]
#   col9 <- combo_IDs[9,i]
#   col10 <- combo_IDs[10,i]
#   col11 <- combo_IDs[11,i]
#   
#   predictor1 <-Climate_Fire_DF[,col1]
#   predictor2 <-Climate_Fire_DF[,col2]
#   predictor3 <-Climate_Fire_DF[,col3]
#   predictor4 <-Climate_Fire_DF[,col4]
#   predictor5 <-Climate_Fire_DF[,col5]
#   predictor6 <-Climate_Fire_DF[,col6]
#   predictor7 <-Climate_Fire_DF[,col7]
#   predictor8 <-Climate_Fire_DF[,col8]
#   predictor9 <-Climate_Fire_DF[,col9]
#   predictor10 <-Climate_Fire_DF[,col10]
#   predictor11 <-Climate_Fire_DF[,col11]
#   
#   #PCA combo model:
#   Covariates_DF = data.frame(predictor1=predictor1,predictor2=predictor2,predictor3=predictor3,predictor4=predictor4,predictor5=predictor5,predictor6=predictor6,predictor7=predictor7,predictor8=predictor8,predictor9=predictor9,predictor10=predictor10,predictor11=predictor11)
#   df.pca <- prcomp(Covariates_DF, center = TRUE,scale. = TRUE)
#   pcs = df.pca$x
#   pc1 = pcs[,1]
#   pc2 = pcs[,2]
#   pc3 = pcs[,3]
#   pc4 = pcs[,4]
#   pc5 = pcs[,5]
#   pc6 = pcs[,6]
#   pc7 = pcs[,7]
#   pc8 = pcs[,8]
#   pc9 = pcs[,9]
#   pc10 = pcs[,10]
#   pc11 = pcs[,11]
#   current_DF_pcs <- data.frame(Y=Y,predictor1=pc1,predictor2=pc2,predictor3=pc3,predictor4=pc4,predictor5=pc5,predictor6=pc6,predictor7=pc7,predictor8=pc8,predictor9=pc9,predictor10=pc10,predictor11=pc11,predictor12=as.vector(Zs))
#   mod <- gam(data=current_DF_pcs,Y~s(predictor1)+s(predictor2)+s(predictor3)+s(predictor4)+s(predictor5)+s(predictor6)+s(predictor7)+s(predictor8)+s(predictor9)+s(predictor10)+s(predictor11)+predictor12,family="gaussian")
#   
#   X = concurvity(mod,full=FALSE)
#   X=X$estimate
#   idx <- which(X==1)
#   X[idx] = 0
#   Max_Concuvrity <- max(X)
#   if (Max_Concuvrity <= concurvity_threshold){
#     AIC = mod$aic
#     R = cor(predict(mod),Y)
#     
#     store_R_11vars = c(store_R_11vars,R)
#     store_AIC_11vars = c(store_AIC_11vars,AIC)
#     store_i_11vars = c(store_i_11vars,i)
#   }
# }
# 
# idx_bestmod_11vars <- which(store_AIC_11vars==min(store_AIC_9vars))
# store_R_11vars[idx_bestmod_11vars]
# idx_bestmod_11vars = store_i_11vars[idx_bestmod_11vars]
# 
# n_predictors <- 12
# ncombos = factorial(nvars)/(factorial(n_predictors)*factorial(nvars-n_predictors))
# x=1:nvars
# combo_IDs<-combn(x,n_predictors)
# 
# store_AIC_12vars = c()
# store_R_12vars = c()
# store_i_12vars = c()
# for (i in 1:ncombos){
#   col1 <- combo_IDs[1,i]
#   col2 <- combo_IDs[2,i]
#   col3 <- combo_IDs[3,i]
#   col4 <- combo_IDs[4,i]
#   col5 <- combo_IDs[5,i]
#   col6 <- combo_IDs[6,i]
#   col7 <- combo_IDs[7,i]
#   col8 <- combo_IDs[8,i]
#   col9 <- combo_IDs[9,i]
#   col10 <- combo_IDs[10,i]
#   col11 <- combo_IDs[11,i]
#   col12 <- combo_IDs[12,i]
#   
#   predictor1 <-Climate_Fire_DF[,col1]
#   predictor2 <-Climate_Fire_DF[,col2]
#   predictor3 <-Climate_Fire_DF[,col3]
#   predictor4 <-Climate_Fire_DF[,col4]
#   predictor5 <-Climate_Fire_DF[,col5]
#   predictor6 <-Climate_Fire_DF[,col6]
#   predictor7 <-Climate_Fire_DF[,col7]
#   predictor8 <-Climate_Fire_DF[,col8]
#   predictor9 <-Climate_Fire_DF[,col9]
#   predictor10 <-Climate_Fire_DF[,col10]
#   predictor11 <-Climate_Fire_DF[,col11]
#   predictor12 <-Climate_Fire_DF[,col12]
#   
#   #PCA combo model:
#   Covariates_DF = data.frame(predictor1=predictor1,predictor2=predictor2,predictor3=predictor3,predictor4=predictor4,predictor5=predictor5,predictor6=predictor6,predictor7=predictor7,predictor8=predictor8,predictor9=predictor9,predictor10=predictor10,predictor11=predictor11,predictor12=predictor12)
#   df.pca <- prcomp(Covariates_DF, center = TRUE,scale. = TRUE)
#   pcs = df.pca$x
#   pc1 = pcs[,1]
#   pc2 = pcs[,2]
#   pc3 = pcs[,3]
#   pc4 = pcs[,4]
#   pc5 = pcs[,5]
#   pc6 = pcs[,6]
#   pc7 = pcs[,7]
#   pc8 = pcs[,8]
#   pc9 = pcs[,9]
#   pc10 = pcs[,10]
#   pc11 = pcs[,11]
#   pc12 = pcs[,12]
#   current_DF_pcs <- data.frame(Y=Y,predictor1=pc1,predictor2=pc2,predictor3=pc3,predictor4=pc4,predictor5=pc5,predictor6=pc6,predictor7=pc7,predictor8=pc8,predictor9=pc9,predictor10=pc10,predictor11=pc11,predictor12=pc12,predictor13=as.vector(Zs))
#   mod <- gam(data=current_DF_pcs,Y~s(predictor1)+s(predictor2)+s(predictor3)+s(predictor4)+s(predictor5)+s(predictor6)+s(predictor7)+s(predictor8)+s(predictor9)+s(predictor10)+s(predictor11)+s(predictor12)+predictor13,family="gaussian")
#   
#   X = concurvity(mod,full=FALSE)
#   X=X$estimate
#   idx <- which(X==1)
#   X[idx] = 0
#   Max_Concuvrity <- max(X)
#   if (Max_Concuvrity <= concurvity_threshold){
#     AIC = mod$aic
#     R = cor(predict(mod),Y)
#     
#     store_R_12vars = c(store_R_12vars,R)
#     store_AIC_12vars = c(store_AIC_12vars,AIC)
#     store_i_12vars = c(store_i_12vars,i)
#   }
# }
# 
# idx_bestmod_12vars <- which(store_AIC_12vars==min(store_AIC_12vars))
# store_R_12vars[idx_bestmod_12vars]
# idx_bestmod_12vars = store_i_12vars[idx_bestmod_12vars]
# store_best_mod_idx = cbind(idx_bestmod_2vars,idx_bestmod_3vars,idx_bestmod_4vars,idx_bestmod_5vars,idx_bestmod_6vars,idx_bestmod_7vars,idx_bestmod_8vars)
# 
# idx_bestmod_2vars <- which(store_AIC_2vars==min(store_AIC_2vars))
# idx_bestmod_3vars <- which(store_AIC_3vars==min(store_AIC_3vars))
# idx_bestmod_4vars <- which(store_AIC_4vars==min(store_AIC_4vars))
# idx_bestmod_5vars <- which(store_AIC_5vars==min(store_AIC_5vars))
# idx_bestmod_6vars <- which(store_AIC_6vars==min(store_AIC_6vars))
# idx_bestmod_7vars <- which(store_AIC_7vars==min(store_AIC_7vars))
# idx_bestmod_8vars <- which(store_AIC_8vars==min(store_AIC_8vars))
# store_best_R = cbind(store_R_2vars[idx_bestmod_2vars],store_R_3vars[idx_bestmod_3vars],store_R_4vars[idx_bestmod_4vars],store_R_5vars[idx_bestmod_4vars],store_R_6vars[idx_bestmod_6vars],store_R_7vars[idx_bestmod_7vars],store_R_8vars[idx_bestmod_8vars])
# store_best_AIC = cbind(store_AIC_2vars[idx_bestmod_2vars],store_AIC_3vars[idx_bestmod_3vars],store_AIC_4vars[idx_bestmod_4vars],store_AIC_5vars[idx_bestmod_5vars],store_AIC_6vars[idx_bestmod_6vars],store_AIC_7vars[idx_bestmod_7vars],store_AIC_8vars[idx_bestmod_8vars])
# 
# 
# idx_bestmod_2vars = store_i_2vars[idx_bestmod_2vars]
# idx_bestmod_3vars = store_i_3vars[idx_bestmod_3vars]
# idx_bestmod_4vars = store_i_4vars[idx_bestmod_4vars]
# idx_bestmod_5vars = store_i_5vars[idx_bestmod_5vars]
# idx_bestmod_6vars = store_i_6vars[idx_bestmod_6vars]
# idx_bestmod_7vars = store_i_7vars[idx_bestmod_7vars]
# idx_bestmod_8vars = store_i_8vars[idx_bestmod_8vars]

#define the best models by AIC:

# vars=c("SWEI","WinPRCP","WinTMP","Spring-prcp","Spring TMP","Spring VPD","Winter VPD","Spring ET","Spring PET","Winter ET","Winter PET","Spring area","winter+spring temp","winter+spring precip","winter+spring VPD","winter+spring PET")

#record best 100 models:
AIC_7vars=store_AIC_7vars
i_7vars=store_i_7vars
store_i=c()
for (b in 1:100){
  idx_bestmod_7vars <- which(AIC_7vars==min(AIC_7vars))
  i_bestmod_7vars = i_7vars[idx_bestmod_7vars]
  store_i =cbind(store_i,i_bestmod_7vars)
  
  AIC_7vars=AIC_7vars[-idx_bestmod_7vars]
  i_7vars=i_7vars[-idx_bestmod_7vars]
}

best_is=store_i
outfilename = sprintf("/Volumes/Pruina_External_Elements/DroughtFireSnow/Data/AnlysisData/Best_Obs_Forecast_outputs/BAs_1984_2020_Obs_Forecast_SE_Zbins_PDSImod_PCA_bestmods_higherPeakSWE.csv")
write.table(best_is,file=outfilename,sep=",",row.names = FALSE)