#infile2<-"dates_interpolation_03012012.txt"  # list of 10 dates for the regression

infile2<-"dates_interpolation_03052012.txt"

prop<-0.3                                                                            #Proportion of testing retained for validation   

infile3<-"LST_dates_var_names.txt"

#######START OF THE SCRIPT #############

  mod <-ghcn.subsets[[i]][,match(LST_dates[i], names(ghcn.subsets[[i]]))]

  ghcn.subsets[[i]] = transform(ghcn.subsets[[i]],LST = mod)

  n<-nrow(ghcn.subsets[[i]])

  ns<-n-round(n*prop)  #Create a sample from the data frame with 70% of the rows

  nv<-n-ns             #create a sample for validation with prop of the rows

  mod4<- gam(tmax~ s(lat) + s (lon) + s(ELEV_SRTM) + s(Northness) + s (Eastness) + s(DISTOC) + s(LST), data=data_s)

  mod5<- gam(tmax~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST), data=data_s)

  mod6<- gam(tmax~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST,LC1), data=data_s)

  ####Regression part 3: Calculating and storing diagnostic measures

  results_AIC[i,1]<- dates[i]  #storing the interpolation dates in the first column

  results_AIC[i,6]<- AIC (mod4)

  results_AIC[i,7]<- AIC (mod5)

  results_AIC[i,8]<- AIC (mod6)

  results_GCV[i,1]<- dates[i]  #storing the interpolation dates in the first column

  results_GCV[i,2]<- ns        #number of stations used in the training stage

  results_GCV[i,6]<- mod4$gcv.ubre

  results_GCV[i,7]<- mod5$gcv.ubre

  results_GCV[i,8]<- mod6$gcv.ubre

  results_DEV[i,1]<- dates[i]  #storing the interpolation dates in the first column

  results_DEV[i,2]<- ns        #number of stations used in the training stage

  results_DEV[i,6]<- mod4$deviance

  results_DEV[i,7]<- mod5$deviance

  results_DEV[i,8]<- mod6$deviance

  #####VALIDATION: Prediction checking the results using the testing data########

  y_mod4<- predict(mod4, newdata=data_v, se.fit = TRUE) 

  y_mod5<- predict(mod5, newdata=data_v, se.fit = TRUE) 

  y_mod6<- predict(mod6, newdata=data_v, se.fit = TRUE)

  res_mod1<- data_v$tmax - y_mod1$fit #Residuals for GMA model that resembles the ANUSPLIN interpolation

  res_mod2<- data_v$tmax - y_mod2$fit   #Residuals for GAM model that resembles the PRISM interpolation                               

  res_mod4<- data_v$tmax - y_mod4$fit

  res_mod5<- data_v$tmax - y_mod5$fit

  res_mod6<- data_v$tmax - y_mod6$fit

  RMSE_mod1 <- sqrt(sum(res_mod1^2)/nv)          

  RMSE_mod2 <- sqrt(sum(res_mod2^2)/nv)

  RMSE_mod4 <- sqrt(sum(res_mod4^2)/nv)

  RMSE_mod5 <- sqrt(sum(res_mod5^2)/nv)

  RMSE_mod6 <- sqrt(sum(res_mod6^2)/nv)

  results_RMSE[i,1]<- dates[i]  #storing the interpolation dates in the first column

  results_RMSE[i,2]<- ns        #number of stations used in the training stage

  results_RMSE[i,6]<- RMSE_mod4

  results_RMSE[i,7]<- RMSE_mod5

  results_RMSE[i,8]<- RMSE_mod6

  #Saving dataset in dataframes

  data_name<-paste("ghcn_v_",dates[[i]],sep="")

  assign(data_name,data_v)

## Plotting and saving diagnostic measures

results_RMSEnum <-results_RMSE

mode(results_RMSEnum)<- "numeric"

# Make it numeric first

# Now turn it into a data.frame...

results_table_RMSE<-as.data.frame(results_RMSEnum)

#results_table_RMSE

# End of script##########