05_MLR_SSN_Model.Rmd

---
title:  "James MLR and SSN Modeling"
author: "Michael McManus & Travis Linscome-Hatfield"
date: "12/03/2024"
output:
  html_document:
    fig_caption: yes
    number_sections: false
    toc: yes
    toc_float:  yes
    code_folding: hide
    self_contained: yes
    theme: lumen
editor_options: 
  chunk_output_type: console
---

# Outline
The purpose of this script is to model the relationship between VSCI and covariates across the 4 geographies of Watershed (Ws), Watershed-Riparian (WsRp), Catchment (Cat), and Catchment-Riparian (CatRp). For each of the 4 geographies, the modeling is done as follows:

* Use bestglm to select covariates 
* Check model diagnostics of the MLR model
* Run a SSN nugget model (SSN0), which is the equivalent to an MLR model, and check diagnostic plots
* Run the SSN1 with covariates and spatial autocovariances. Note the spatial autocovariances are identified by running script 03_Y_Intercept_SSN_Model.
* Compare SSN0 to SSN1 to see if modeling spatial autocorrelation of residuals improved model fit and model prediction
* Plot loocv observed and predicted
* Map residuals
* Test for random effects of trend stations versus status stations. Trends stations are sampled repeated, with station_id_2 being the identifier that distinguishes trend from status stations. Status stations are only sampled one time.
* Compare SSN models across the 4 geographies using ML estimation method as the models have different covariates.

Portions of this analysis are summarized in the Presentation folder, specifically 2024_sfc_philadelphia_bioassessment_session_mcmanus.pdf. That presentation provides an overview of the detailed analysis presented below.

Log Keeping log below for now when reach final version will delete.
11/21/2024 Only retained modeling chunks. Loading SSN that was produced in script 03.

On 10/15/2024 added waterbody data for Travis to evaluate on github.

On 09/19/2024 Ellen provided update waterbody datasets as csv files.

09/11/2024 Ellen and I worked out how to get NHD waterbody associated with 199 monitoring observations. 

08/28/2024 explored STATSGO_Set1 and STATSGO_Set2 variables from StreamCat. Also, brought in waterbodies data created by Ellen D'Amico from Pegasus. Ellen used comid from edges, but I am not sure they all match to the featureid of the points.


On 08/15/2024 added data, figures_sfs, and ssn_object folders and trying to run code pointing to those folders.

On 08/14/2024 pushing to github repository.

On 07/31/2024 bringing in WWTP_VA.vcs to do a left merge with DFobs so I can compare stations with WWTPs to stations without WWTPs in the VSCI response.

On 07/11/2024 bring the new ssn into R.
On 07/10/2024 downloaded James_071024.ssn that had PRISM climate data, including normals and deviations from normals.

On 05/14/2024 running code to get SSN output for SSN models for the 4 geographies and because 05/08/2024 updated R, RTools, and RStudio.

I am going to bring in two SSNs, the 040424 version and the latest 041924 version. I will take the PRISM climate variables out of 041924 and put into 040424.
Downloaded James_041924.ssn_NEWEST.zip and it is stored in James_041924.ssn with file modified on 04/22/2024 at 7:34am. Ellen removed netid and pid columns and reran through SSN compilation. Bestglm still gave error

On 04/19/2024 imported James041924.ssn, still could not get bestglm to run.

On 04/17/2024 imported James041724.ssn that now include annual maximum temperature, annual average temperature, and annual average precipitation. Also, I will simplify coding of wetlands to be absence/presence across all 4 geographies:  watershed, watershed-riparian, catchment, and catchment-riparian. For grass cover at catchment and catchment-riparian I will also use binary code for absence/presence.

On 04/16/2024 imported James041924.ssn, which now has PRISM climate data from StreamCat.

This code will only use vsci as the response variable. Also, I will examine if impervious surface, forest, hay, and grass at 0 values show dramatically different vsci y-intercept, when x = 0, compared to vsci y-intercepts when x > 0.

On 04/09/2024, I had a call with Mike Dumelle, and he suggested 1) use untransformed response variable, 2) include binary absence/presence along with continuous landscape percentages, and 3) try including year.

On 04/04/2024 now using James040424.ssn as Ellen D'Amico from Pegasus updated Preds_2021_2022 points. I entered WQ data as described in section 1.10 Missing Water Chemistry Data below.

On 01/09/2023 read in Pegasus' SSN object from 01/09/2024, which has additional StreamCat variables. I am interested in the shrub and wetland metrics. In previous SSN, Note that VDEQ switched fonts between Wadeable_ProbMon_2001-2018_Final_Final.xlsx
and Wadeable_ProbMon_2001-2020.xlsx. I used janitor package to put all variable names in lowercase.


On 12/26/2023 read in Pegasus' SSN object from 11/16/2023, which as fields of COMID, AreaSqkM, TotDASqKM, h2oAreaKm2, and rcaAreaKM2. Using this new SSN object means I re-entered missing water chemistry data for two stations in ArcGIS.

Now using SSN2 from CRAN on 11/16/2023.

Note on SSNs I receive from Pegasus, I have to hand enter some of the water chemistry values for stations 2-WLS023.10 and 2AMLC000.84. This is because Emma Jones at VDEQ sent FinalFieldAnalytes_MissingJamesSites.csv by email on 09/29/2023. That csv file has water chemistry entries for those two stations, which in the original data file were NAs. I edit those entries in ArcGIS before I do import into R. Once the SSN is in R I can't export back to ArcGIS. 

On 11/09/2023 unzipped. James.ssn_110923 and saved to James.ssn. This James.ssn_11092023 has additonal covariates produced by Ellen D'Amico, GIS analyst with Pegasus.
The James.ssn_102323 had stream order added to edges.
On 11/08/2023 unzipped James.ssn_102323 into ssn_objects folder and renamed it as James.ssn so I would not need to change R code.
Can I use SSN2 package for some of ESDA of James River? Yes.
# Library and Functions Load
```{r setup, collapse=TRUE}

library(tidyverse)
library(janitor)
library(sf)
library(dummy)
library(bestglm)
library(performance)
library(car)
library(SSN2)
library(mapview)
library(leafpop) # for popups in mapview
library(leafsync) # to sync obs and prediction maps

sessionInfo()

# attached base packages:
# [1] stats     graphics  grDevices utils     datasets  methods  
# [7] base     
# 
# other attached packages:
#  [1] leafsync_0.1.0     leafpop_0.1.0      mapview_2.11.2    
#  [4] SSN2_0.2.1         car_3.1-2          carData_3.0-5     
#  [7] performance_0.12.2 see_0.8.5          bestglm_0.37.3    
# [10] leaps_3.1          dummy_0.1.3        sf_1.0-16         
# [13] janitor_2.2.0      lubridate_1.9.3    forcats_1.0.0     
# [16] stringr_1.5.1      dplyr_1.1.4        purrr_1.0.2       
# [19] readr_2.1.5        tidyr_1.3.1        tibble_3.2.1      
# [22] ggplot2_3.5.1      tidyverse_2.0.0  

knitr::opts_chunk$set(message=FALSE, warning=FALSE,collapse = T)

```

# 1.0 Load SSN
```{r ssn_from_script_03}
load("outputs/j_ssn3.RData")

DFobs <- SSN2::ssn_get_data(j_ssn3) %>% clean_names(.)
names(DFobs)
```


# 2.0 Dummy Coding for Bestglm Models
The bestglm function only accepts dummy coding for factors.

# 2.1 Dummy Coding of Factors
Dummy coding for ecoregions, Virginia hydrologic unit subbasins, wetlands, impervious surfaces, and total phosphorus. See the description of dummy coding at <https://stats.oarc.ucla.edu/r/library/r-library-contrast-coding-systems-for-categorical-variables/>.
```{r dummy_code}

# remove geometry as bestglm does not want it
DFobsz <- st_set_geometry(DFobs, NULL)
#creating alternate dummy version for base case comparison if desired
DFobsz_alternate_dummy = DFobsz

# dummy code 5 ecoregions with base being Blue Ridge Mtns
eco_r <- dplyr::select(DFobsz, eco_region)
summary(eco_r)
glimpse(eco_r)

eco_d <- (data.frame(dummy(eco_r)))
class(eco_d)
dim(eco_d)
head(eco_d)
str(eco_d)
distinct(eco_d)
# 5 levels need only n-1 =4 dummy variables, removed base level of Blue Ridge Mountains by dropping first column
eco_d <- eco_d[c(-1)]
dim(eco_d)
head(DFobsz$eco_region)
distinct(eco_d)
head(eco_d)
str(eco_d)
class(eco_d)

DFobsz <- cbind(DFobsz,eco_d)
names(DFobsz)

# dummy code 5 vahusb with base being JU, James Upper
vahusb <- dplyr::select(DFobsz, vahusb)
summary(vahusb)
glimpse(vahusb)

vahusb_d <- (data.frame(dummy(vahusb)))
# 5 levels need only n-1 =4 dummy variables, removed base level of JU by dropping first column
vahusb_d <- vahusb_d[c(-1)]
dim(vahusb_d)
head(DFobsz$vahusb)
distinct(vahusb_d)
head(vahusb_d)
str(vahusb_d)
class(vahusb_d)

DFobsz <- cbind(DFobsz,vahusb_d)
names(DFobsz)


# dummy code 4 watershed wetland levels with base being none
wet_w_f <- dplyr::select(DFobsz, pct_wet_w_f)
summary(wet_w_f)
glimpse(wet_w_f)

wet_w_d <- (data.frame(dummy(wet_w_f)))
# 4 levels need only n-1 =3 dummy variables, removed base level of none by dropping first column
wet_w_d <- wet_w_d[c(-1)]
DFobsz <- cbind(DFobsz,wet_w_d)
names(DFobsz)

# dummy code 4 watershed-riparian wetland levels with base being none
wet_wr_f <- dplyr::select(DFobsz, pct_wet_wr_f)
summary(wet_wr_f)
glimpse(wet_wr_f)

wet_wr_d <- (data.frame(dummy(wet_wr_f)))
# 4 levels need only n-1 =3 dummy variables, removed base level of none by dropping first column
wet_wr_d <- wet_wr_d[c(-1)]
DFobsz <- cbind(DFobsz,wet_wr_d)
names(DFobsz)

# dummy code 4 catchment wetland levels with base being none
wet_c_f <- dplyr::select(DFobsz, pct_wet_c_f)
summary(wet_c_f)
glimpse(wet_c_f)

wet_c_d <- (data.frame(dummy(wet_c_f)))
# 4 levels need only n-1 =3 dummy variables, removed base level of none by dropping first column
wet_c_d <- wet_c_d[c(-1)]
DFobsz <- cbind(DFobsz,wet_c_d)
names(DFobsz)

# dummy code 3 catchment-riparian wetland levels (none, low, high) with base being none
wet_cr_f <- dplyr::select(DFobsz, pct_wet_cr_f)
summary(wet_cr_f)
glimpse(wet_cr_f)

wet_cr_d <- (data.frame(dummy(wet_cr_f)))
# 3 levels need only n-1 =2 dummy variables, removed base level of none by dropping first column
wet_cr_d <- wet_cr_d[c(-1)]
DFobsz <- cbind(DFobsz,wet_cr_d)
names(DFobsz)

# dummy code 3 catchment-riparian impervious surfaces levels with base being none
imp_cr_f <- dplyr::select(DFobsz, pct_imp_rp_c_f)
summary(imp_cr_f)
glimpse(imp_cr_f)

imp_cr_d <- (data.frame(dummy(imp_cr_f)))
# 3 levels need only n-1 =2 dummy variables, removed base level of none by dropping first column
imp_cr_d <- imp_cr_d[c(-1)]
DFobsz <- cbind(DFobsz,imp_cr_d)
names(DFobsz)

# dummy code 3 total phosphorus levels with base being low
tp_f <- dplyr::select(DFobsz, tp_f)
summary(tp_f)
glimpse(tp_f)

tp_d <- (data.frame(dummy(tp_f)))
# 3 levels need only n-1 =2 dummy variables, removed base level of none by dropping first column
tp_d <- tp_d[c(-1)]
DFobsz <- cbind(DFobsz,tp_d)
names(DFobsz)

# remove files not needed downstream
rm(eco_d, eco_r, imp_cr_d, imp_cr_f, tp_d, tp_f, vahusb, vahusb_d, wet_c_d, wet_c_f, wet_cr_d, wet_cr_f, wet_w_d, wet_w_f, wet_wr_d, wet_wr_f)
```

# 2.15 Dummy Coded to SF to SSN
The 23 dummy coded variables just created need to be put into an SF object and then, later, into an SSN object. 
```{r dummy_sf_ssn}
names(DFobsz)

DFobsz2 <- dplyr::select(DFobsz, c(station_id_2, eco_region_Central.Appalachian.Ridges.and.Valleys:tp_f_high))

# put dummy covariates in an SF object
DFobs2 <- full_join(DFobs, DFobsz2, by = join_by(station_id_2))
class(DFobs2)
names(DFobs2)


```





# 3.0 SSN Ws, WsRp, Cat, & CatRp Models
I need Neptune to vet code below across all 4 geographies. Before run check_model function have Plots tab fully expanded vertically. The check_model functions provides sisual check of various model assumptions. The added variable plots, also called partial regression plots, are made using the avPlots function and give a visualization of the slope estimates.

The candidate pool of 21 covariates came from these categories:  1) climate 2) landcover, 3) instream measurements, and 4) geophysical or hydrologic units.

The maximum number of covariates to be retained during bestglm selection was set to 4 before any modeling was done. This was done to produce a conservative model based on the number of observations relative to the number of parameters estimated in the model. Four slopes and 1 y-intercept are estimated with nvmax = 4. Then, a spatial stream network model could include up to 7 additional parameters, with 2 parameters estimated each for Euclidean, Tail Up, and Tail down autocovariances plus an independent error, or nugget, parameter. For 199 observations, there could be 12 model parameters or ~ 17 observations/parameter, and that seemed a good way to avoid overfitting the model as a rule of thumb is to have 10 observations/parameter.

# 4.0 Ws-Wq
Analysis done at watershed extent.
### 4.1 Bestglm and SSN0 Nugget Model
On 07/12/2024 I will simply keep precip_mm and tmean in the candidate pool for all 4 geographies. No change to original selected covariates after I replaced precip_mm with d_precip_m and tmean with d_t_mean, which are deviations from 30-year normals for precip and temperature.
On 04/22/2024 changed WsWq so that precip_mm and tmean are now included, but they were not selected by bestglm. The previous wetland factor with multiple levels was replaced with a binary absence/presence variable. 
The X1 matrix below includes that dummy factors with the station_id_2 enables the join the dummy factor to an sf object and then to the SSN. 


```{r ws_bestglm}

# candidate covariates
WsWq <- DFobsz|>
  dplyr::select(precip_mm, tmean, pct_imp_w_emplog, pct_for_w_emplog, pct_hay_w_emplog, pct_grs_w_emplog, elev_ws, do, p_h, tothab, l_spc, l_tds, l_turb, l_tn, vahusb_JM, vahusb_JR, vahusb_JA, vahusb_JL, bin_wet_w, tp_f_medium, tp_f_high)
str(WsWq)
# X can only contain numeric or factor
WsWq<- WsWq %>% mutate_at(c('vahusb_JM', 'vahusb_JR', 'vahusb_JA', 'vahusb_JL','tp_f_medium', 'tp_f_high'), as.numeric)

# bind matrix of covariates to response variable
WsWqy<-cbind.data.frame(WsWq,y=DFobsz$vsci)
ASRcum_wswq1.best<-bestglm(WsWqy,IC = "AIC", nvmax =4, TopModels=15)

# See what terms are in best models
ASRcum_wswq1.best$BestModels
ASRcum_wswq1.best$Subsets
summary(ASRcum_wswq1.best$BestModel)

# create mlr model so diagnostics can be run on that lm object
wswq_mlry2 <- lm(vsci ~ pct_imp_w_emplog + elev_ws + do + vahusb_JL, data = DFobsz)
summary(wswq_mlry2)
check_model(wswq_mlry2)
# if you wish to view avPlots uncomment below. But this hangs if you attempt to run the whole script as one block
# avPlots(wswq_mlry2)

# SSN0 nugget model
# station_id_2 lets join_by work
X1 <- DFobsz|>
   dplyr::select(station_id_2, vahusb_JM, vahusb_JR, vahusb_JA, vahusb_JL, pct_wet_w_f_low, pct_wet_w_f_medium, pct_wet_w_f_high, tp_f_medium, tp_f_high)

# X can only contain numeric or factor
X1 <- X1 %>% mutate_at(c('vahusb_JM', 'vahusb_JR', 'vahusb_JA', 'vahusb_JL', 'pct_wet_w_f_low', 'pct_wet_w_f_medium', 'pct_wet_w_f_high', 'tp_f_medium', 'tp_f_high'), as.numeric)
str(X1)
# put dummy covariates into SF object by full join to original DFobs
DFobs3a <- full_join(DFobs, X1, by = join_by(station_id_2))
names(DFobs3a)

# put SF object into SSN
j_ssn3 <-  SSN2::ssn_put_data(DFobs3a,j_ssn3)

# ssn0 Watershed
ssn0_wswq_mlr <- ssn_lm(
  formula = vsci ~ pct_imp_w_emplog + elev_ws + do + vahusb_JL,
  ssn.object = j_ssn3,
  tailup_type = "none",
  taildown_type = "none",
  euclid_type = "none",
  nugget_type = "nugget",
  estmethod = "reml" #default
)
summary(ssn0_wswq_mlr)
varcomp(ssn0_wswq_mlr)
loocv(ssn0_wswq_mlr)

#to run script continuously from top to bottom comment out plot statement below as requires repeated hitting of enter key, which skips script lines in doing so
# plot(ssn0_wswq_mlr, which = c(1:6))
```




### 4.2 SSN0 Nugget vs SSN1
This chunk compares nugget model, typical nonspatial MLR or SSN0, to the SSN1 model that specifies the spatial autocovariances that were identified in script 03_Y-Intercept_SSN_Model. The comparison is made using model fit diagnostics, AICc, and leave-one-out-cross-validation (LOOCV) metrics of the root mean squared prediction error (RMSPE) and the squared correlation between the observed value and its leave-one-out prediction (cor2).

Is Torgegram on flow-unconnected standardized residuals telling me that still some structured spatial autocorrelation left?

On 04/16/2024, re-ran R_vdeq_tutdeu_comparison.Rmd and models_yintercept AICc statistic and loocv(tdExp_eu_exp) RMSPE and cor2 showed that tail down exponential and Euclidean exponential was the best of the 15 spatial autocovariance models compared. For some reason I had been specifying tail down exponential and Euclidean Gaussian. I corrected so both TD and EU are exponential. The R_vdeq_tutdeu_comparison.Rmd file is at E:\R_vdeq_sci.
```{r nugget_vs_spatial_ws}
ssn_wswq_reml1 <- ssn_lm(
  formula = vsci ~ pct_imp_w_emplog + elev_ws + do + vahusb_JL,
  ssn.object = j_ssn3,
  tailup_type = "none",
  taildown_type = "exponential",
  euclid_type = "exponential",
  nugget_type = "nugget",
  estmethod = "reml",
  additive = "afv_area"
)
summary(ssn_wswq_reml1)
varcomp(ssn_wswq_reml1)
loocv(ssn_wswq_reml1)

#to run script continuously from top to bottom comment out plot statement below as requires repeated hitting of enter key, which skips script lines in doing so
# plot(ssn_wswq_reml1, which = c(1:6))

models_ssn0_ssn1ws <- glances(ssn0_wswq_mlr, ssn_wswq_reml1)
print(models_ssn0_ssn1ws)

# check with Travis if code below can be dropped. Not sure it is showing me much
aug_ssn_wswq_reml1 <- augment(ssn_wswq_reml1, drop = FALSE)
class(aug_ssn_wswq_reml1)
mapview(aug_ssn_wswq_reml1)

ggplot(aug_ssn_wswq_reml1, aes(x = pct_imp_w_emplog, y = .fitted, colour = vahusb)) + geom_point() + geom_smooth(method = "lm")

ggplot(aug_ssn_wswq_reml1, aes(x = vahusb, y = .fitted)) + geom_boxplot()

names(aug_ssn_wswq_reml1)

resid_ssn1 <- as.data.frame(aug_ssn_wswq_reml1)|>
   dplyr::select(station_id_2, .fitted, .resid, .std.resid)

class(resid_ssn1)
names(resid_ssn1)
# put dummy covariates in an SF object
DFobs4 <- full_join(DFobs3a, resid_ssn1, by = join_by(station_id_2))
names(DFobs4)

# put SF object into SSN
j_ssn4 <-  SSN2::ssn_put_data(DFobs4,j_ssn3)

ztg <- SSN2::Torgegram(.std.resid ~ 1, j_ssn4, type = c("flowuncon", "euclid"))

plot(ztg, main = "Torgegram Standardized Residuals")
class(ztg)
names(ztg)
ztg[["flowuncon"]]
ztg[["euclid"]]
# summary(ztg$flowcon)
summary(ztg$flowuncon)
summary(ztg$euclid)
plot(ztg)
# plot(ztg, type = "flowcon", main="Residuals")
plot(ztg, type = "flowuncon", main = "Torgegram Standardized Residuals")
plot(ztg, type = "euclid", main = "Torgegram Standardized Residuals")
class(ztg)
# plot(ztg, separate = TRUE, main = "TotHab")


```



## Synced map of Obs, Fitted, Standardized Residuals
Other than showing syncing can be done, why are these maps helpful?
```{r synced_maps}
# 
# mapview(aug_ssn_wswq_reml1, zcol = ".std.resid", cex = ".std.resid", alpha.regions = .8, legend = TRUE, layer.name = 'Standardized Residuals', popup = popupTable(aug_ssn_wswq_reml1, zcol = c("vsci", ".fitted", ".resid", ".std.resid")))
# 
# # observed vsci at 199 obs
# obs1 <- mapview(aug_ssn_wswq_reml1, zcol = "vsci", cex = "vsci", alpha.regions = .8, legend = TRUE, layer.name = 'VSCI', popup = popupTable(aug_ssn_wswq_reml1, zcol = c("vsci", ".fitted", ".resid", ".std.resid"))) 
# 
# # fitted vsci at 199 obs
# fit1 <- mapview(aug_ssn_wswq_reml1, zcol = ".fitted", cex = ".fitted", alpha.regions = .8, legend = TRUE, layer.name = 'Fitted', popup = popupTable(aug_ssn_wswq_reml1, zcol = c("vsci", ".fitted", ".resid", ".std.resid")))
# 
# # standardized residuals at 199 obs
# stdres1 <- mapview(aug_ssn_wswq_reml1, zcol = ".std.resid", cex = ".std.resid", alpha.regions = .8, legend = TRUE, layer.name = 'Std. Residual', popup = popupTable(aug_ssn_wswq_reml1, zcol = c("vsci", ".fitted", ".resid", ".std.resid")))
# 
# sync(obs1,fit1, ncol=1) # obs & predicted
# sync(fit1, stdres1, ncol=1)
# sync(obs1, fit1, stdres1, ncol=1)
```

### 4.3 Random Effects
The ssn_wswq_reml1 has the smallest AICc, but largest difference in AICc is only 2.4.

I want Neptune to check this code.
I believe the the random effect for station_id is coded correctly as it recognizes that there are 8 stations with repeated 4 to 5 observations whereas all the other stations have single observation. The st_type, for study type, I believe is not coded correctly as it only recognizes two groups, status or trend, stations so that code is commented out. Based on AICc and cor2, ssn_wswq_reml1 is still a better model. 
```{r random_effect_station_id}
ssn_wswq_rand1 <- ssn_lm(
  formula = vsci ~ pct_imp_w_emplog + elev_ws + do + vahusb_JL,
  ssn.object = j_ssn3,
  tailup_type = "none",
  taildown_type = "exponential",
  euclid_type = "exponential",
  nugget_type = "nugget",
  estmethod = "reml",
  additive = "afv_area",
  random = ~ (1 | station_id)
)

summary(ssn_wswq_rand1)
varcomp(ssn_wswq_rand1)
loocv(ssn_wswq_rand1)

# if want to run script continuously from top to bottom comment out plot statement below as requires repeated hitting of enter key
# plot(ssn_wswq_rand1, which = c(1:6))

models_rand1 <- glances(ssn0_wswq_mlr, ssn_wswq_reml1, ssn_wswq_rand1)
print(models_rand1)
# ssn_wswq_rand2 <- ssn_lm(
#   formula = y2 ~ pct_imp_w_emplog + elev_ws + do + vahusb_JL,
#   ssn.object = j_ssn3,
#   tailup_type = "none",
#   taildown_type = "exponential",
#   euclid_type = "gaussian",
#   nugget_type = "nugget",
#   estmethod = "reml",
#   additive = "afv_area",
#   random = ~ (1 | st_type)
# )
# 
# summary(ssn_wswq_rand2)
# varcomp(ssn_wswq_rand2)
# loocv(ssn_wswq_rand2)
# 
# models_rand2 <- glances(ssn0_wswq_mlr, ssn_wswq_reml1, ssn_wswq_rand1)


```

# 5.0 WsRp-Wq
Analysis done at Watershed-Riparian extent.
### 5.1 Bestglm and SSN0 Nugget Model
On 04/22/2024 changes WsWq so that precip_mm and tmean are now included. The previous wetland factor with multiple levels was replaced with a binary absence/presence variable. 
On 04/17/2024 caught typo in WsRpWQ as "pct_hay_w_emplog" listed and corrected to "pct_hay_wr_emplog".
On 04/16/2024 running models on vsci, and pct_imp_rp_w_emplog, elev_ws, do, and tothab selected.
```{r wsriparian_bestglm}
WsRpWq <- DFobsz|>
  dplyr::select(precip_mm, tmean, pct_imp_rp_w_emplog, pct_for_wr_emplog, pct_hay_wr_emplog, pct_grs_wr_emplog, elev_ws, do, p_h, tothab, l_spc, l_tds, l_turb, l_tn,  vahusb_JM, vahusb_JR, vahusb_JA, vahusb_JL, bin_wet_wr, tp_f_medium, tp_f_high)
str(WsRpWq)
# X can only contain numeric or factor
WsRpWq<- WsRpWq %>% mutate_at(c('vahusb_JM', 'vahusb_JR', 'vahusb_JA', 'vahusb_JL', 'tp_f_medium', 'tp_f_high'), as.numeric)

# bind matrix of covariates to response variable
WsRpWqy<-cbind.data.frame(WsRpWq,y=DFobsz$vsci)
ASRcum_wsrpwq1.best<-bestglm(WsRpWqy,IC = "AIC",nvmax =4, TopModels=15)

# See what terms are in best models
ASRcum_wsrpwq1.best$BestModels
ASRcum_wsrpwq1.best$Subsets
summary(ASRcum_wsrpwq1.best$BestModel)

wsrpwq_mlry2 <- lm(vsci ~ pct_imp_rp_w_emplog + elev_ws + do + tothab, data = DFobsz)
summary(wsrpwq_mlry2)
check_model(wsrpwq_mlry2)
# if you wish to view avPlots uncomment below. But this hangs if you attempt to run the whole script as one block
# avPlots(wsrpwq_mlry2)

# SSN0 nugget model
# station_id_2 lets join_by work
X2 <- DFobsz|>
  dplyr::select(station_id_2, vahusb_JM, vahusb_JR, vahusb_JA, vahusb_JL, pct_wet_wr_f_low, pct_wet_wr_f_medium, pct_wet_wr_f_high, tp_f_medium, tp_f_high)

# X can only contain numeric or factor
X2 <- X2 %>% mutate_at(c('vahusb_JM', 'vahusb_JR', 'vahusb_JA', 'vahusb_JL', 'pct_wet_wr_f_low', 'pct_wet_wr_f_medium', 'pct_wet_wr_f_high', 'tp_f_medium', 'tp_f_high'), as.numeric)

# put dummy covariates in an SF object
DFobs3a <- full_join(DFobs, X2, by = join_by(station_id_2))
names(DFobs3a)

# put SF object into SSN
j_ssn3 <-  SSN2::ssn_put_data(DFobs3a,j_ssn3)

# SSN0 Watershed-Riparian
ssn0_wsrpwq_mlr <- ssn_lm(
  formula = vsci ~ pct_imp_rp_w_emplog + elev_ws + do + tothab,
  ssn.object = j_ssn3,
  tailup_type = "none",
  taildown_type = "none",
  euclid_type = "none",
  nugget_type = "nugget",
  estmethod = "reml"
)
summary(ssn0_wsrpwq_mlr)
varcomp(ssn0_wsrpwq_mlr)
loocv(ssn0_wsrpwq_mlr)

#to run script continuously from top to bottom comment out plot statement below as requires repeated hitting of enter key, which skips script lines in doing so
# plot(ssn0_wsrpwq_mlr, which = c(1:6))
```

### 5.2 SSN0 Nugget vs SSN1
```{r nugget_vs_spatial_wsriparian}
ssn_wsrpwq_reml1 <- ssn_lm(
  formula = vsci ~ pct_imp_rp_w_emplog + elev_ws + do + tothab,
  ssn.object = j_ssn3,
  tailup_type = "none",
  taildown_type = "exponential",
  euclid_type = "exponential",
  nugget_type = "nugget",
  estmethod = "reml",
  additive = "afv_area"
)
summary(ssn_wsrpwq_reml1)
varcomp(ssn_wsrpwq_reml1)
loocv(ssn_wsrpwq_reml1)


#to run script continuously from top to bottom comment out plot statement below as requires repeated hitting of enter key, which skips script lines in doing so
# plot(ssn_wsrpwq_reml1, which = c(1:6))

models_ssn0_ssn1wsrp <- glances(ssn0_wsrpwq_mlr, ssn_wsrpwq_reml1)
print(models_ssn0_ssn1wsrp)

# does augmenting help?
aug_ssn_wsrpwq_reml1 <- augment(ssn_wsrpwq_reml1, drop = FALSE)
class(aug_ssn_wsrpwq_reml1)
mapview(aug_ssn_wsrpwq_reml1)


```


### 5.3 Random Effects: Ws_Riparian
```{r wsrpwq_random}
ssn_wsrpwq_rand1 <- ssn_lm(
  formula = vsci ~ pct_imp_rp_w_emplog + elev_ws + do + tothab,
  ssn.object = j_ssn3,
  tailup_type = "none",
  taildown_type = "exponential",
  euclid_type = "exponential",
  nugget_type = "nugget",
  estmethod = "reml",
  additive = "afv_area",
  random = ~ (1 | station_id)
)

summary(ssn_wsrpwq_rand1)
varcomp(ssn_wsrpwq_rand1)
loocv(ssn_wsrpwq_rand1)

#to run script continuously from top to bottom comment out plot statement below as requires repeated hitting of enter key, which skips script lines in doing so
# plot(ssn_wsrpwq_rand1, which = c(1:6))

models_rand2 <- glances(ssn0_wsrpwq_mlr, ssn_wsrpwq_reml1, ssn_wsrpwq_rand1)
print(models_rand2)

```

### ssn_wsrpwq_reml1 Mapview residuals 
Is map of observed, fitted, orstandardized residuals, helpful?

```{r maps_ssn_wsrpwq_reml1}
mapview(aug_ssn_wsrpwq_reml1, zcol = ".std.resid",  alpha.regions = .8, legend = TRUE, layer.name = 'Standardized Residuals', popup = popupTable(aug_ssn_wsrpwq_reml1, zcol = c("vsci", ".fitted", ".resid", ".std.resid")))
```

# 6.0 Cat-Wq
Analysis done at Catchment extent.
### 6.1 Bestglm and SSN0 Nugget Model
On 04/22/2024 changes CatWq so that precip_mm and tmean are now included. The previous wetland factor with multiple levels was replaced with a binary absence/presence variable as was grass cover. VSCI values at pct_grs_c_emplog at 0 encompasses range of VSCI values at pct_grs_c_emplog > 0. 
On 04/16/2024 running model on vsci and specifying exponential for tail down and Euclidean.
```{r cat_bestglm}
CatWq <- DFobsz|>
  dplyr::select(precip_mm, tmean, pct_imp_c_emplog, pct_for_c_emplog, pct_hay_c_emplog, bin_grs_c, elev_cat, do, p_h, tothab, l_spc, l_tds, l_turb, l_tn, vahusb_JM, vahusb_JR, vahusb_JA, vahusb_JL, bin_wet_c,  tp_f_medium, tp_f_high)
str(CatWq)
# X can only contain numeric or factor
CatWq<- CatWq %>% mutate_at(c('vahusb_JM', 'vahusb_JR', 'vahusb_JA', 'vahusb_JL', 'bin_wet_c','tp_f_medium', 'tp_f_high'), as.numeric)

# bind matrix of covariates to response variable
CatWqy<-cbind.data.frame(CatWq,y=DFobsz$vsci)
ASRcum_catwq1.best<-bestglm(CatWqy,IC = "AIC",nvmax =4, TopModels=15)

# See what terms are in best models
ASRcum_catwq1.best$BestModels
ASRcum_catwq1.best$Subsets
summary(ASRcum_catwq1.best$BestModel)

catwq_mlry2 <- lm(vsci ~ pct_imp_c_emplog + elev_cat + tothab + l_tn, data = DFobsz)
summary(catwq_mlry2)
check_model(catwq_mlry2)
# if you wish to view avPlots uncomment below. But this hangs if you attempt to run the whole script as one block
# avPlots(catwq_mlry2)

# No categorical covariates needed so commented out X3
# SSN0 nugget model
# st_id_tren lets join_by work
# X3 <- DFobsz|>
#  dplyr::select(st_id_tren, vahusb_JM, vahusb_JR, vahusb_JA, vahusb_JL, pct_wet_c_f_low, pct_wet_c_f_medium, pct_wet_c_f_high, tp_f_medium, tp_f_high)

# X can only contain numeric or factor
# X3 <- X3 %>% mutate_at(c('vahusb_JM', 'vahusb_JR', 'vahusb_JA', 'vahusb_JL', 'pct_wet_c_f_low', 'pct_wet_c_f_medium', 'pct_wet_c_f_high','tp_f_medium', 'tp_f_high'), as.numeric)

# put dummy covariates in an SF object
# DFobs3a <- full_join(DFobs2a, X3, by = join_by(st_id_tren))
# names(DFobs3a)
# put transformed covariates in an SF object
# DFobs2a <- full_join(DFobs2, DFobsz2, by = join_by(st_id_tren))

# put SF object into SSN
# j_ssn3 <-  SSN2::ssn_put_data(DFobs3a,j_ssn3)

# SSN0 Catchment
ssn0_catwq_mlr <- ssn_lm(
  formula = vsci ~ pct_imp_c_emplog + elev_cat + tothab + l_tn,
  ssn.object = j_ssn3,
  tailup_type = "none",
  taildown_type = "none",
  euclid_type = "none",
  nugget_type = "nugget",
  estmethod = "reml" #default
)
summary(ssn0_catwq_mlr)
varcomp(ssn0_catwq_mlr)
loocv(ssn0_catwq_mlr)

#to run script continuously from top to bottom comment out plot statement below as requires repeated hitting of enter key, which skips script lines in doing so
# plot(ssn0_catwq_mlr, which = c(1:6))
```

### 6.2 SSN0 Nugget vs SSN1
```{r nugget_vs_spatial_cat}
ssn_catwq_reml1 <- ssn_lm(
  formula = vsci ~ pct_imp_c_emplog + elev_cat + tothab + l_tn,
  ssn.object = j_ssn3,
  tailup_type = "none",
  taildown_type = "exponential",
  euclid_type = "exponential",
  nugget_type = "nugget",
  estmethod = "reml",
  additive = "afv_area"
)
summary(ssn_catwq_reml1)
varcomp(ssn_catwq_reml1)
loocv(ssn_catwq_reml1)


#to run script continuously from top to bottom comment out plot statement below as requires repeated hitting of enter key, which skips script lines in doing so
# plot(ssn_catwq_reml1, which = c(1:6))

models_ssn0_ssn1cat <- glances(ssn0_catwq_mlr, ssn_catwq_reml1)
print(models_ssn0_ssn1cat)

aug_ssn_catwq_reml1 <- augment(ssn_catwq_reml1, drop = FALSE)
class(aug_ssn_catwq_reml1)
mapview(aug_ssn_catwq_reml1)

```

### 6.3 Random Effects: Cat-Wq
Check random effects of station_id on ssn_catwq_rand1.
```{r catwq_random}
ssn_catwq_rand1 <- ssn_lm(
  formula = vsci ~ pct_imp_c_emplog + elev_cat + tothab + l_tn,
  ssn.object = j_ssn3,
  tailup_type = "none",
  taildown_type = "exponential",
  euclid_type = "exponential",
  nugget_type = "nugget",
  estmethod = "reml",
  additive = "afv_area",
  random = ~ (1 | station_id)
)
summary(ssn_catwq_rand1)
varcomp(ssn_catwq_rand1)
loocv(ssn_catwq_rand1)

#to run script continuously from top to bottom comment out plot statement below as requires repeated hitting of enter key, which skips script lines in doing so
# plot(ssn_catwq_rand1, which = c(1:6))

models_rand3 <- glances(ssn0_catwq_mlr, ssn_catwq_reml1,ssn_catwq_rand1)
print(models_rand3)
```

# 7.0 CatRp-Wq
Analysis done at Catchment-Riparian extent.
On 04/22/2024 changes CatRpWq so that precip_mm and tmean are now included. The previous wetland factor with multiple levels was replaced with a binary absence/presence variable as was grass cover.
With this smallest geography (by extent and configuration), I have to shift some land cover from continuous to categorical. Grass at this geography has a median value of 0 so I made an absence/presence variable for grass cover so a form of that covariate is in the model.
### 7.1 Bestglm and SSN0 Nugget Model
```{r catriparian_bestglm}
CatRpWq <- DFobsz|>
  dplyr::select(precip_mm, tmean, pct_imp_rp_c_emplog, pct_for_cr_emplog, pct_hay_cr_emplog, bin_grs_cr, elev_cat, do, p_h, tothab, l_spc, l_tds, l_turb, l_tn,  vahusb_JM, vahusb_JR, vahusb_JA, vahusb_JL, bin_wet_cr, tp_f_medium, tp_f_high)
str(CatRpWq)
# X can only contain numeric or factor
CatRpWq<- CatRpWq %>% mutate_at(c('vahusb_JM', 'vahusb_JR', 'vahusb_JA', 'vahusb_JL',
'tp_f_medium', 'tp_f_high'), as.numeric)

# bind matrix of covariates to response variable
CatRpWqy<-cbind.data.frame(CatRpWq,y=DFobsz$vsci)
ASRcum_catrpwq1.best<-bestglm(CatRpWqy,IC = "AIC",nvmax =4, TopModels=15)

# See what terms are in best models
ASRcum_catrpwq1.best$BestModels
ASRcum_catrpwq1.best$Subsets
summary(ASRcum_catrpwq1.best$BestModel)

catrpwq_mlry2 <- lm(vsci ~ pct_imp_rp_c_emplog + elev_cat + tothab + l_tn, data = DFobsz)
summary(catrpwq_mlry2)
check_model(catrpwq_mlry2)
# if you wish to view avPlots uncomment below. But this hangs if you attempt to run the whole script as one block
# avPlots(catrpwq_mlry2)

# No categorical predictors selected so comment out code
# SSN0 nugget model
# st_id_tren lets join_by work
# X4 <- DFobsz|>
#   dplyr::select(st_id_tren, vahusb_JM, vahusb_JR, vahusb_JA, vahusb_JL, pct_wet_cr_f_low, pct_wet_cr_f_high, tp_f_medium, tp_f_high)

# X can only contain numeric or factor
# X4 <- X4 %>% mutate_at(c('vahusb_JM', 'vahusb_JR', 'eco_region_Piedmont', 'vahusb_JA', 'pct_wet_cr_f_low', 'pct_wet_cr_f_high', 'tp_f_medium', 'tp_f_high'), as.numeric)

# put dummy covariates in an SF object
# DFobs3a <- full_join(DFobs2a, X4, by = join_by(st_id_tren))
# names(DFobs3a)
# put transformed covariates in an SF object
# DFobs2a <- full_join(DFobs2, DFobsz2, by = join_by(st_id_tren))

# put SF object into SSN
# j_ssn3 <-  SSN2::ssn_put_data(DFobs3a,j_ssn3)

# SSN0 Catchment
ssn0_catrpwq_mlr <- ssn_lm(
  formula = vsci ~ pct_imp_rp_c_emplog + elev_cat + tothab + l_tn,
  ssn.object = j_ssn3,
  tailup_type = "none",
  taildown_type = "none",
  euclid_type = "none",
  nugget_type = "nugget",
  estmethod = "reml" #default
)
summary(ssn0_catrpwq_mlr)
varcomp(ssn0_catrpwq_mlr)
loocv(ssn0_catrpwq_mlr)

#to run script continuously from top to bottom comment out plot statement below as requires repeated hitting of enter key, which skips script lines in doing so
# plot(ssn0_catrpwq_mlr, which = c(1:6))

```

### 7.2 SSN0 vs SSN1
```{r nugget_vs_spatial_catriparian}
ssn_catrpwq_reml1 <- ssn_lm(
  formula = vsci ~ pct_imp_rp_c_emplog + elev_cat + tothab + l_tn,
  ssn.object = j_ssn3,
  tailup_type = "none",
  taildown_type = "exponential",
  euclid_type = "exponential",
  nugget_type = "nugget",
  estmethod = "reml",
  additive = "afv_area"
)
summary(ssn_catrpwq_reml1)
varcomp(ssn_catrpwq_reml1)
loocv(ssn_catrpwq_reml1)

#to run script continuously from top to bottom comment out plot statement below as requires repeated hitting of enter key, which skips script lines in doing so
# plot(ssn_catrpwq_reml1, which = c(1:6))

models_ssn0_ssn1catrp <- glances(ssn0_catrpwq_mlr, ssn_catrpwq_reml1)
print(models_ssn0_ssn1catrp)

aug_ssn_catrpwq_reml1 <- augment(ssn_catrpwq_reml1, drop = FALSE)

# mapview(aug_ssn_catrpwq_reml1)

```

### 7.3 Random Effects: CatRp-Wq
```{r catrpwq_random}
ssn_catrpwq_rand1 <- ssn_lm(
  formula = vsci ~ pct_imp_rp_c_emplog + elev_cat + tothab + l_tn,
  ssn.object = j_ssn3,
  tailup_type = "none",
  taildown_type = "exponential",
  euclid_type = "exponential",
  nugget_type = "nugget",
  estmethod = "reml",
  additive = "afv_area",
  random = ~ (1 | station_id)
)
summary(ssn_catrpwq_rand1)
varcomp(ssn_catrpwq_rand1)
loocv(ssn_catrpwq_rand1)

#to run script continuously from top to bottom comment out plot statement below as requires repeated hitting of enter key, which skips script lines in doing so
# plot(ssn_catrpwq_reml1, which = c(1:6))

models_rand4 <- glances(ssn0_catrpwq_mlr, ssn_catrpwq_reml1, ssn_catrpwq_rand1)
print(models_rand4)
```

# 8.0 Compare Four Geography SSN Models with ML
This to compare different SSN models from the 4 geographies. Given that the models have different covariates then the estimation method needs to be maximum likelihood or ML for AICc comparisons. See An Introduction to Spatial Stream Network Modeling in R using SSN2 for a description of when REML vs ML is used. That discussion is the paragraph just above section 4.1.1 Model Summaries. Note to evaluate the four models in predictive performance the loocv function is used and comparisons of RMSPE and cor2 are made.

The watershed model, ssn_wswq_ml1, fits the data better and predicts better than the 3 other models. We will now use the watershed model for predictions in script 06.
```{r ssn_ml}
ssn_wswq_ml1 <- ssn_lm(
  formula = vsci ~ pct_imp_w_emplog + elev_ws + do + vahusb_JL,
  ssn.object = j_ssn3,
  tailup_type = "none",
  taildown_type = "exponential",
  euclid_type = "exponential",
  nugget_type = "nugget",
  estmethod = "ml",
  additive = "afv_area"
)

ssn_wsrpwq_ml1 <- ssn_lm(
  formula = vsci ~ pct_imp_rp_w_emplog + elev_ws + do + tothab,
  ssn.object = j_ssn3,
  tailup_type = "none",
  taildown_type = "exponential",
  euclid_type = "exponential",
  nugget_type = "nugget",
  estmethod = "ml",
  additive = "afv_area"
)

ssn_catwq_ml1 <- ssn_lm(
  formula = vsci ~ pct_imp_c_emplog + elev_cat + tothab + l_tn,
  ssn.object = j_ssn3,
  tailup_type = "none",
  taildown_type = "exponential",
  euclid_type = "exponential",
  nugget_type = "nugget",
  estmethod = "ml",
  additive = "afv_area"
)

ssn_catrpwq_ml1 <- ssn_lm(
  formula = vsci ~ pct_imp_rp_c_emplog + elev_cat + tothab + l_tn,
  ssn.object = j_ssn3,
  tailup_type = "none",
  taildown_type = "exponential",
  euclid_type = "exponential",
  nugget_type = "nugget",
  estmethod = "ml",
  additive = "afv_area"
)

models_ml_ssn <- glances(ssn_wswq_ml1, ssn_wsrpwq_ml1, ssn_catwq_ml1, ssn_catrpwq_ml1)

print(models_ml_ssn)

loocv(ssn_wswq_ml1)
loocv(ssn_wsrpwq_ml1)
loocv(ssn_catwq_ml1)
loocv(ssn_catrpwq_ml1)


```

# Neptune Code from TLH
This code is included for thoroughness, but does not need to be run unless needed.

In addition to reviewing R code, Neptune statistician Travis Linscome-Hatfield provided code chunks that evaluated 1) assigning a different reference group, rather than the default, and 2) did having a different reference group result in a different model being selected by the bestglm function. Running these two chunks resulted in the same model being selected based on the default group of JU.

## TLH Chunk 1 Dummy Coding of Factors-Travis
Travis Linscome-Hatfield at Neptune provide the code chunk below to see if changing the reference group altered the analysis results and conclusions. This evaluation was motivated by the observation that the default reference group for the Virginia hydrologic unit subbasins was JU, the upper James.
Dummy coding for ecoregions, Virginia hydrologic unit subbasins, wetlands, impervious surfaces, and total phosphorus.
```{r dummy_code-Travis}

names(DFobsz_alternate_dummy)
# the 3 commented lines below simply allow for variable definition for quick funciton adjustment.
dat = DFobsz_alternate_dummy
alt_dat = DFobsz_alternate_dummy
variable = "eco_region"
category_break = function(dat,alt_dat,variable,base1=NULL,base2=NULL){
  
  #pull the column we wish to dummy code
  var_r <- dplyr::select(dat, all_of(variable))
  # get the value with the highest number of instances
  var_b <- as.data.frame(table(var_r))
  var_b <- var_b[order(-var_b$Freq),]
  # if base cases aren't already defined
  if(is.null(base1)){
    # pick base case
    base1 = as.character(var_b[1,1])
  }
  
  # if base2 isn't defined 
  if(is.null(base2)){
    if(base1 == as.character(var_b[1,1])){
      base2 = as.character(var_b[2,1])
    }else{
      base2 = as.character(var_b[1,1])
    }
  }

  
  # dummying the variable
  var_d <- (data.frame(dummy(var_r)))
  # identify which column is the base case
  # turning column names into something that can be directly matched with the base case
  col_short = colnames(var_d)
  col_short = gsub(variable,"",col_short)
  col_short = gsub("_","",col_short)
  col_short = gsub("\\."," ",col_short)
  base1_col_id = which(base1 == col_short)
  base2_col_id = which(base2 == col_short)

  
  # remove the base case from var_d
  var_d1 = var_d[,-base1_col_id]
  var_d2 = var_d[,-base2_col_id]
  # binding into apropriate datasets
  dat = cbind(dat,var_d1)
  alt_dat = cbind(alt_dat,var_d2)
  
  out_list = list("Variable"=variable,
                  "base1"=base1,
                  "base2"=base2,
                  "dat"=dat,
                  "alt_dat"=alt_dat)
  return(out_list)
  
}

# dummy data for eco_region
df_dum_eco = category_break(dat= DFobsz_alternate_dummy,
                            alt_dat=DFobsz_alternate_dummy,
                            "eco_region")
# dummy data for vahusb
df_dum_vahusb = category_break(df_dum_eco$dat,
                               df_dum_eco$alt_dat,
                               "vahusb")
# dummy data for wetland (pct_wet_w_f)
levels(DFobsz_alternate_dummy$pct_wet_w_f)
df_dum_wet = category_break(dat = df_dum_vahusb$dat,
                            alt_dat = df_dum_vahusb$alt_dat,
                            variable = "pct_wet_w_f",
                            base1 = "none")
# dummy data for watershed-ripeaian wetland (pct_wet_wr_f)
levels(DFobsz_alternate_dummy$pct_wet_wr_f)
df_dum_wet_wr = category_break(df_dum_wet$dat,
                               df_dum_wet$alt_dat,
                               "pct_wet_wr_f",
                               base1 = "none")
# dummy data for catchment wetland levels (pct_wet_c_f)
levels(DFobsz_alternate_dummy$pct_wet_c_f)
df_dum_wet_c = category_break(df_dum_wet_wr$dat,
                              df_dum_wet_wr$alt_dat,
                              "pct_wet_c_f",
                              base1 = "none")
# dummy data for catchment-riparian wetland levels (pct_wet_cr_f)
levels(DFobsz_alternate_dummy$pct_wet_cr_f)
df_dum_wet_cr = category_break(df_dum_wet_c$dat,
                               df_dum_wet_c$alt_dat,
                               "pct_wet_cr_f",
                               base1 = "none")

# dummy code for catchment-riparian impervious surfaces levels (pct_imp_rp_c_f)
levels(DFobsz_alternate_dummy$pct_imp_rp_c_f)
df_dum_imp_rp = category_break(df_dum_wet_cr$dat,
                               df_dum_wet_cr$alt_dat,
                               "pct_imp_rp_c_f",
                               base1 = "none")

# dummy code for total phosphorus levels (tp_f)
levels(DFobsz_alternate_dummy$tp_f)
df_dum_tp = category_break(df_dum_imp_rp$dat,
                           df_dum_imp_rp$alt_dat,
                           "tp_f",
                           base1 = "low")

dat_dummy = df_dum_tp$dat
alt_dat_dummy = df_dum_tp$alt_dat

```

## TLH Chunk 2 Bestglm and SSN0 Nugget Model-Travis
On 01/23/2025 commented out TLH code as it showed obtained same model as chunk above.

The code below evaluates if a different reference group was chosen would a qualitatively different model result from the best glm selection. For example, if the reference watershed used was vahusb_JM instead of vahusb_JU.

```{r ws_bestglm-Travis}

#double checking our base cases for variables we're selecting for the bestglm
df_dum_vahusb$base1
df_dum_vahusb$base2
df_dum_tp$base1
df_dum_tp$base2
# initial selection
WsWq <- dat_dummy|>
  dplyr::select(precip_mm, tmean, pct_imp_w_emplog, pct_for_w_emplog, pct_hay_w_emplog, pct_grs_w_emplog, elev_ws, do, p_h, tothab, l_spc, l_tds, l_turb, l_tn, vahusb_JU, vahusb_JR, vahusb_JA, vahusb_JL, bin_wet_w, tp_f_medium, tp_f_high)
str(WsWq)
# X can only contain numeric or factor
WsWq<- apply(WsWq, 2, as.numeric)

WsWqy<-cbind.data.frame(WsWq,y=DFobsz$vsci)
ASRcum_wswq1.best<-bestglm(WsWqy,IC = "AIC", nvmax =4, TopModels=15)

# See what terms are in best models
ASRcum_wswq1.best$BestModels
ASRcum_wswq1.best$Subsets
summary(ASRcum_wswq1.best$BestModel)

wswq_mlry2 <- lm(vsci ~ pct_imp_w_emplog + elev_ws + do + vahusb_JL, data = DFobsz)
summary(wswq_mlry2)
check_model(wswq_mlry2)
# if you wish to view avPlots uncomment below. But this hangs if you attempt to run the whole script as one block
#avPlots(wswq_mlry2)

## testing alternate base case to see if any issues arise
WsWq2 <- alt_dat_dummy|>
  dplyr::select(precip_mm, tmean, pct_imp_w_emplog, pct_for_w_emplog, pct_hay_w_emplog, pct_grs_w_emplog, elev_ws, do, p_h, tothab, l_spc, l_tds, l_turb, l_tn, vahusb_JM, vahusb_JR, vahusb_JA, vahusb_JL, bin_wet_w, tp_f_low, tp_f_high)

WsWq2<- apply(WsWq2, 2, as.numeric)

WsWqy2<-cbind.data.frame(WsWq2,y=DFobsz$vsci)
ASRcum_wswq2.best<-bestglm(WsWqy2,IC = "AIC", nvmax =4, TopModels=15)
summary(ASRcum_wswq1.best$BestModel)

# Perfect we get the same model we had before, no issues there.
```