06_SSN_Model_Predict.Rmd

---
title: "James SSN Predictions"
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
---

Here's what I added in this script for the example!
# Outline
The purpose of this script is to make sure both obs and preds have the same covariates. This is required so that the SSN model built from the obs can then applied to make predictions at the preds. The 48 preds where withheld from from the SSN model building to provide an independent data set for measuring the predictive performance of the SSN model. The 48 preds were collected from spatially balanced probabilistic surveys done from 2019-2022.

Three modifications were made to the SSN WS-WQ model, ssn_wswq_reml1, from script 05 in attempts to get a better predictive model. Those modifications included:
1) Added a partition factor statement to the model based on parallel coordinate plot,
2) Based on the literature, tested 3 additional covariates to try and get a better fitting and predictive model,
3) Tested the hypothesis that monitoring sites having upstream waterbodies furthest away would have higher VSCI.

With those modification evaluated, then the final model, which included the partition factor and the nearest distance to upstream waterbody covariate, was evaluated by comparing observed vsci to predicted vsci, first, for the 199 observations from 2001-2018 using leave-one-out-cross-validation (loocv). The second evaluation used the 48 predictions sites from 2019-2022 with their observed vsci compared to the predicted vsci from the SSN model. Both of those evaluations were done graphs of the observed values on the y-axis and predicted values on the x-axis.

Log Keeping log below for now when reach final version will delete.

12/23/2024 Travis walked me through merging his interval_pred_upd branch with main repository on GitHub. Once that was done I pulled that from main repository to my local repository.


12/16/2024 trying to run distance hypothesis using binomial, numerically coded waterbody upstream, 0 is absent and 1 is present, and distance upstream to nearest waterbody as continuous covariate. 

11/25/2024 using James_071024_pluspreds as ssn object as has all 48 prediction points in one shapefile. And bringing additional covariates. Needing to assign 0's to preds_2019_2022 that initially had dist_nearest_up_wb_km set to NA for 31 of 48 prediction points.

11/19/2024 trying to run TLH code chunk for prediction interval graph

On 10/24/2024, went through the code with Travis and now runs without errors all chunks above Status vs Trend Predictions chunk. Pushed to repository.

Trying again on 08/14/2024 to push this code up to github repository.
On 08/13/2024, pushed to github. Trying push again as I don't see this code
on github. 

On 07/12/2024 need to make sure that james_071024.ssn, with PRISM data, gives same prediction output as previous ssn object. It does.

On 07/02/2024 R version 4.4.1 was installed.
On 07/01/2024 added code to include partition factor assignment to preds.

On 06/28/2024 added code to run a parition factor model.

On 06/17/2024 running code to get predictions for 2019-2020 that I put into a geopackage so I map and plot se vs prediction in QGIS.

On 05/14/2024 running code as updated R, RTools, and RStudio on 05/08/2024.

On 04/19/2024 also specify spatial autocovariance of tail down exponential and Euclidean exponential. Also, run MLR of mixed geography using impervious cover of all 4 geographies (to identify high collinearity), elevation at watershed, DO, tothab, and l_tn, and Lower James subbasin. Does that combination of variables better fit and outperform watershed-only model?

On 04/18/2024 now using vsci as response variable not y2.

On 04/17/2024 now using James041724.ssn as it has climate PRISM data out of StreamCat of precip_mm, Tmean, and Tmax.

On 04/09/2024 Mike Dumelle recommended:  1) use untransformed vsci, 2) include absence/presence binary land cover with percent land cover, and 3) include year.

On 04/04/2024 now using James040424.ssn as Ellen D'Amico from Pegasus updated Preds_2021_2022 points.

On 03/27/2024, Using code from step1_james_eda_v1.Rmd, but now specifying preds1 to bring in 22 points from 2019-2020 VDEQ monitoring and preds2 to bring in 26 points from 2021-2022 VDEQ monitoring using current SSN from James.ssn_032224.zip.
# Libraries
```{r setup,message=FALSE, warning=FALSE, collapse=TRUE}

library(SSN2)
library(tidyverse)
library(dplyr)
library(Hmisc)
library(sf)
library(moments)
library(fitdistrplus)
library(mapview)
library(leaflet)
library(leafpop) # for popups in mapview
library(leafsync) # sync mapview maps
library(units)
library(janitor)
library(dummy)
library(GGally)  # parallel coordinate plot
library(readxl)
library(bestglm)
library(car)
library(performance)
library(scales) # comma instead of scientific notation
library(corrr) # tidyverse correlation package
library(tmap)
library(gstat) # semivariogram cloud
library(dotwhisker)
# remotes::install_github("fsolt/dotwhisker") 
# library(dotwhisker) #coefficient plot did not install as archived
# consider https://github.com/ggobi/ggally/issues/372
# with function in Ggally of GGally::ggcoef() 
sessionInfo()

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

# SSN Observations and Preds
Notice this import explicitly brings in the prediction points. Also, this version of the SSN requires that both the obs and predpts to have the same covariates.
```{r ssn_obs_preds}
# now bringing predpts 2019-2022 combined
# current ssn
j_ssn1a <- SSN2::ssn_import("ssn_object/James_071024_pluspreds.ssn", predpts = c("Preds_2019_2022"),
overwrite = FALSE)

# previous ssn
# j_ssn1a <- SSN2::ssn_import("E:/R_vdeq_sci/Working/Data/neptune_analysis/ssn_objects/James041724.ssn", predpts = c("Preds_2019_2020", "Preds_2021_2022"),
#overwrite = FALSE)

# j_ssn1a <- SSN2::ssn_import("E:/R_vdeq_sci/Working/Data/neptune_analysis/ssn_objects/James.ssn", predpts = c("Preds_2019_2020", "Preds_2021_2022"),
# overwrite = FALSE)
names(j_ssn1a)
summary(j_ssn1a)

DFobs <- SSN2::ssn_get_data(j_ssn1a) %>% clean_names(.)

preds1 <- SSN2::ssn_get_data(j_ssn1a, name = "Preds_2019_2022") %>% clean_names(.)

# str(preds1)
# str(preds2)


```

## SFS Meeting Torgegram
```{r sfs_torgegram}

# torg2 <- ggplot(flowuncon, aes(x=dist, y=gamma)) + geom_point(size =3) + labs(x = "Flow-Unconnected Stream Distance (m)", y = "Semivariance", title = "VSCI Torgegram Stream Network Distance") + scale_x_continuous(labels=comma)
# png(file="figures_sfs/torg2.png",width=6,height=3,units="in",res=150)
#   torg2
# dev.off()


# esv2 <- ggplot(euclid, aes(x=dist, y=gamma)) + geom_point(size =3) + labs(x = "Euclidean Distance (m)", y = "Semivariance", title = "VSCI Semivariogram Euclidean Distance") + scale_x_continuous(labels=comma)
# png(file="figures_sfs/esv2.png",width=6,height=3,units="in",res=150)
#   esv2
# dev.off()
```


## Steps for working with preds based on obs

Note that both obs and preds need to be in the SSN object that is modeled. This results in the SSN model object having all the variables needed for the SSN predict function.
	1. Make factors of vahusb and year.
	2. Log total nitrogen. 
	3. Include additional covariates of wastewater treatment plant, depth to bedrock, upstream waterbody as categorical absent or present, and continuous covariate of distance to nearest upstream waterbody.
	4. Empirical logit transformations for impervious and forest
	5. Dummy code transformation of factors for vahusb.

On 03/04/2029, realized I needed to add in imputing of 2 observed sites for total habitat scores. I need tot_hab for catchment and catchment-riparian in the mixed geographies SSN analysis. Saw this when I made scatter plot of covariates.

## 1.20 Total Habitat (RBP) score
Bring Total Habitat Score (TotHab) in from Wadeable_ProbMon_2001-2018_Final_Final.xslx spreadsheet so it can be joined to DFobs and then j_ssn1a. These 2 stations:  2-JKS070.97 and 2-DDY000.75_2017 do not have tothab as Emma confirmed in her 11/24/2023 email. Both sites are in Central Appalachian Ridges and Valleys. Both have high VSCI scores of 73.8 and 84.5 (the latter is max VSCI), respectively. I will impute their tothab scores. For trend station 2-DDY000.75_2017, I will average the scores of 172.5, 173.5, and 178 from 2011, 2013, and 2015, respectively. For  2-JKS070.97 on 3rd order, I averaged nearby sites 2-JKS076.16, has tothab of 162.0 on 3rd order, is upstream of 2-JKS070.97, about 16 km apart. Also used Back Creek site, 2-BCC001.90 (has tothab of 189.0 on 2nd order),  that flows parallel to Jackson River, where the 2 sites are, and Back Creek site is near confluence to Jackson River.  2-BCC001.90 is about 9 stream km from 2-JKS070.97. Added new variable tothab.

```{r tothab}

tothab_ds1 <- read_xlsx("data/Wadeable_ProbMon_2001-2018_Final_Final.xlsx", range = "Wadeable_ProbMon_2001-2018!D1:BK814")
                        
tothab_ds2 <- tothab_ds1 |>
  filter(SubBasin == "James") |>
  dplyr::select(StationID_Trend, TotHab) |>
  mutate_at(c('TotHab'), as.numeric)

summary(tothab_ds2$TotHab)

# 2-DDY000.75_2017 is on Daddy Run headwater of Calfpasture River
# https://stackoverflow.com/questions/32829358/dplyr-filter-with-sql-like-wildcard
dr_na <- filter(tothab_ds2, grepl("2-DDY000.75", StationID_Trend, fixed = TRUE))
summary(dr_na$TotHab)

tothab_ds3 <- tothab_ds2|>
  mutate(
    TotHab = case_when(StationID_Trend == "2-DDY000.75_2017" ~ 176.4,
       TRUE ~ TotHab))


# 2-JKS070.97 is on Jackson River
jr_na <- filter(tothab_ds3, StationID_Trend == "2-JKS076.16"| StationID_Trend == "2-BCC001.90")

summary(jr_na$TotHab)

tothab_ds4 <- tothab_ds3|>
  mutate(
    TotHab = case_when(StationID_Trend == "2-JKS070.97" ~ 175.5,
       TRUE ~ TotHab))
# no longer any NAs
summary(tothab_ds4$TotHab)

tothab_ds4 <- rename(tothab_ds4,c(tothab = TotHab, st_id_tren = StationID_Trend ))
head(tothab_ds4)

# add tothab as new variable
DFobs <- full_join(DFobs, tothab_ds4, by = join_by(station_id_2==st_id_tren))

names(DFobs)

# remove datasets not needed downstream
rm(tothab_ds1, tothab_ds2, dr_na, tothab_ds3, jr_na, tothab_ds4)
```

# Transform Obs Covariates
On 12/16/2024 added binomial wbc, wbc_bin, coded as numeric 0 or 1 similar to Dumelle et al. 2023 coding in National Lake assessment conductivity paper.
 
On 11/25/2024 added WWTP, STATSGo_Set2, and waterbody covariates.
0n 06/28/20204 added for ju - yes and ju - no to test if partitioning on ju helpful first by looking at separate torgegrams and then in modeling.
On 04/10/2024 added code for a binary presence covariate for impervious cover at watershed extent.
```{r transform_obs}

DFobs$year_f <- as.factor(DFobs$year)
DFobs$vahusb <- factor(DFobs$vahusb, levels = c("JU", "JM", "JR", "JA", "JL"))
summary(DFobs$vahusb)

DFobs <- DFobs %>%
  mutate(
    jl = case_when(
    vahusb != "JL" ~ 0,
    vahusb == "JL" ~1
  )
)

DFobs <- DFobs %>% 
  mutate(
    ju = case_when(
      vahusb == "JU" ~ "yes",
      .default = "no"
    )
  )
DFobs$ju <- factor(DFobs$ju, levels = c("yes","no"))
summary(DFobs$ju)

glimpse(DFobs)
# In ArcGIS these 2 fields are not numeric so have to mutate
DFobs2 <- DFobs %>% mutate_at(c('pct_for_c', 'pct_for_w'), as.numeric)

DFobs2$l_tn <- log(DFobs2$tn)

DFobs2$vsci <- round(DFobs2$vscivcpmi,1)

# WWTP data
wwtp_ds1 <- read.csv("data/WWTP_VA.csv") %>% clean_names(.)

DFobs2 <- left_join(DFobs2, wwtp_ds1, by = join_by(feature_id == comid))

DFobs2 <- DFobs2 %>% 
  mutate(
    wwtp = case_when(
      wwtp_all_dens_ws > 0 ~ "yes",
      wwtp_all_dens_ws == 0 ~ "no"
    )
  )

DFobs2$wwtp <- factor(DFobs2$wwtp, levels = c("yes","no"))
summary(DFobs2$wwtp)

# Statgo2 Set 2 variables, use RckDepWs
statsgo_set2 <- read.csv("data/STATSGO_Set2_VA.csv") %>% clean_names(.)

names(statsgo_set2)

DFobs2 <- left_join(DFobs2, statsgo_set2, by = join_by(feature_id == comid))

# Waterbody data
wb_ds1 <- read.csv("data/ObservationPoints_DistancesUpstream_100424.csv") %>% clean_names(.)
class(wb_ds1)
names(wb_ds1)
head(wb_ds1)
str(wb_ds1)
summary(wb_ds1$num_waterbody_up, na.rm = TRUE)
summary(wb_ds1$dist_nearest_up_wb_km, na.rm = TRUE)
wb_ds1$dist_nearest_up_wb_km <- round(wb_ds1$dist_nearest_up_wb_km, 1)

summary(wb_ds1$dist_nearest_up_wb_km, na.rm = TRUE)
# file below provides a link to waterbody, and subsequently to DFobs data.
StationIDs_UniqLocIDs <- read.csv("data/StationIDs_UniqLocIDs.csv") %>% clean_names(.)

names(StationIDs_UniqLocIDs)

wb_ds2 <- left_join(StationIDs_UniqLocIDs, wb_ds1, by=join_by(uniq_loc_id))
names(wb_ds2)

# from help by = join_by(name == artist)
# station_id_2 and station_id are both long station trend names
DFobs2 <- left_join(DFobs2, wb_ds2, by = join_by(station_id_2==station_id))
names(DFobs2)
DFobs2 <- DFobs2 %>% 
  mutate(
    wbc = case_when(
      incl_nhdwb == "yes" ~ "present",
      incl_nhdwb != "yes" ~ "absent"
    )
  )

class(DFobs2$wbc)

DFobs2$wbc <- factor(DFobs2$wbc, levels = c("present", "absent"))
summary(DFobs2$wbc)

# adding binomial wbc, wbc_bin, coded as numeric 0 or 1 similar to Dumelle et al. 2023 coding in National Lake assessment conductivity paper.

DFobs2 <- DFobs2 %>% 
  mutate(
    wbc_bin = case_when(
      incl_nhdwb == "yes" ~ 1,
      incl_nhdwb != "yes" ~ 0
    )
  )

class(DFobs2$wbc_bin)
DFbos2 <- as.numeric(DFobs2$wbc_bin)
class(DFobs2$bin_wbc)

summary(DFobs2$dist_nearest_up_wb_km)

DFobs2$dist_nearest_up_wb_km[is.na(DFobs2$dist_nearest_up_wb_km)] <- 0

summary(DFobs2$dist_nearest_up_wb_km)

names(DFobs2)

class(DFobs2)
# note ssn_put_data requires sf object and SSN2 object
j_ssn2 <-  SSN2::ssn_put_data(DFobs2,j_ssn1a)
# just doing this assignment so not have to rename objects
j_ssn3 <- j_ssn2

# VARIABLE ADJUSTMENT ZONE 4
### Variables to apply empirical logit transformation
emplog_vars <- c("pct_for_w","pct_imp_w","pct_crop_w","pct_hay_w","pct_grs_w","pct_shrb_w","pct_for_wr","pct_imp_rp_w","pct_crop_wr","pct_hay_wr","pct_grs_wr","pct_shrb_wr","pct_for_c","pct_imp_c","pct_crop_c","pct_hay_c","pct_grs_c", "pct_shrb_c", "pct_for_cr","pct_imp_rp_c","pct_crop_cr","pct_hay_cr","pct_grs_cr","pct_shrb_cr")

# remove geometry so empirical logit can be applied
DFobsz <- st_set_geometry(DFobs2, NULL)
################################################################
################################################################


## transform these variables and put the new values into new columns in DFobs
for(var in emplog_vars){
  ### create new tranformed data column to preserve the original
  new_nm <- paste0(var,"_emplog")
  dat_vec_obs <- DFobsz[,var]
  # dat_vec_preds <- DFpreds[,var]
  
  #converting to 0-1 range
  dat_vec_obs <- dat_vec_obs/100
  # dat_vec_preds <- dat_vec_preds/100
  
  # dat_vec[dat_vec == 1] <- .9999
  # dat_vec[dat_vec == 0] <- .0001
  
  if(any(dat_vec_obs > 1 | dat_vec_obs < 0)){
    cat("ERROR: percentage variables outside logical bounds")
  }
  
  small_dat_vec_obs <- dat_vec_obs[dat_vec_obs <1 & dat_vec_obs >0]
  # small_dat_vec_preds <- dat_vec_preds[dat_vec_preds <1 & dat_vec_preds >0]
  op1_obs <- small_dat_vec_obs
  op2_obs<- 1-small_dat_vec_obs
  # op1_preds <- small_dat_vec_preds
  # op2_preds <- 1-small_dat_vec_preds
  
  ## minimum of op1 op2
  delt_obs <- min(c(op1_obs,op2_obs))
  # delt_preds <- min(c(op1_preds,op2_preds))
  
  ## getting set of frequencies
  freqs_obs <- NULL
  for(i in 1:length(dat_vec_obs)){
    if(dat_vec_obs[i] <= delt_obs){
      freqs_obs[i] <- delt_obs/2
    }else if(dat_vec_obs[i] >= 1- delt_obs){
      freqs_obs[i] <- 1-(delt_obs/2)
    }else{
      freqs_obs[i] <-dat_vec_obs[i]
    }
  }
  
 # freqs_preds <- NULL
 # for(i in 1:length(dat_vec_preds)){
 #   if(dat_vec_preds[i] <= delt_preds){
 #     freqs_preds[i] <- delt_preds/2
 #   }else if(dat_vec_preds[i] >= 1- delt_preds){
 #     freqs_preds[i] <- 1-(delt_preds/2)
 #   }else{
 #     freqs_preds[i] <-dat_vec_preds[i]
  #  }
 # }
  
  ##getting logits
  logits_obs <- log(freqs_obs/(1-freqs_obs))
  DFobsz[,new_nm] <- logits_obs
  
  # logits_preds <- log(freqs_preds/(1-freqs_preds))
  # DFpreds[,new_nm] <- logits_preds
}

DFobsz <- DFobsz %>%
  mutate(
    imp_w_pres = case_when(
    pct_imp_w_emplog <= -8.111428 ~ 0,
    pct_imp_w_emplog > -8.111428 ~1
  )
)

head(DFobsz$imp_w_pres)
head(DFobsz$pct_imp_w)

names(DFobsz)
DFobsz2 <- dplyr::select(DFobsz, c(station_id_2, pct_for_w_emplog:imp_w_pres))
# put transformed covariates in an SF object
DFobs2a <- full_join(DFobs2, DFobsz2, by = join_by(station_id_2))

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

# 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)

X1 <- DFobsz|>
   dplyr::select(station_id_2, vahusb_JM, vahusb_JR, vahusb_JA, vahusb_JL) %>%
  mutate_at(c('vahusb_JM', 'vahusb_JR', 'vahusb_JA', 'vahusb_JL'), as.numeric)
# X1 can only contain numeric or factor

str(X1)

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

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

# Transform Pred Covariates
On 12/16/2024 added binomial wbc, wbc_bin, coded as numeric 0 or 1 similar to Dumelle et al. 2023 coding in National Lake assessment conductivity paper.

On 11/25/2024 have to assign 0's for 31 of 48 preds that had distances of NA and additional covariates.


```{r transform_preds}
# did this assignment so would not have to change all the object names below
preds3 <- preds1

# create same name as in obs models
preds3$tothab <- preds3$tot_hab

preds3$year_f <- as.factor(preds3$year)
preds3$vahusb <- factor(preds3$vahusb, levels = c("JU", "JM", "JR", "JA", "JL"))
summary(preds3$vahusb)

preds3 <- preds3 %>%
  mutate(
    jl = case_when(
    vahusb != "JL" ~ 0,
    vahusb == "JL" ~1
  )
)

preds3 <- preds3 %>% 
  mutate(
    ju = case_when(
      vahusb == "JU" ~ "yes",
      .default = "no"
    )
  )
preds3$ju <- factor(preds3$ju, levels = c("yes","no"))
summary(preds3$ju)


preds3$l_tn <- log(preds3$tn)

preds3$vsci <- round(preds3$vscivcpmi,1)

# WWTP data loaded from obs so comment  out
# wwtp_ds1 <- read.csv("E:/R_vdeq_nhdplus/WWTP_VA.csv") %>% clean_names(.)

preds3 <- left_join(preds3, wwtp_ds1, by = join_by(feature_id == comid))

preds3 <- preds3 %>% 
  mutate(
    wwtp = case_when(
      wwtp_all_dens_ws > 0 ~ "yes",
      wwtp_all_dens_ws == 0 ~ "no"
    )
  )

preds3$wwtp <- factor(preds3$wwtp, levels = c("yes","no"))
summary(preds3$wwtp)

# loaded from obs so commented out
# Statgo2 Set 2 variables, use RckDepWs
#statsgo_set2 <- read.csv("E:/R_vdeq_nhdplus/STATSGO_Set2_VA.csv") %>% clean_names(.)

names(statsgo_set2)

preds3 <- left_join(preds3, statsgo_set2, by = join_by(feature_id == comid))

# Preds Waterbody data
preds1wb_ds1 <- read.csv("data/PredictionPoints_2019_2020_DistancesUpstream_091924.csv") %>% clean_names(.)

preds2wb_ds1 <- read.csv("data/PredictionPoints_2021_2022_DistancesUpstream_091924.csv") %>% clean_names(.)

predsallwb_ds1 <- bind_rows(preds1wb_ds1, preds2wb_ds1)

class(predsallwb_ds1)
names(predsallwb_ds1)
head(predsallwb_ds1)
str(predsallwb_ds1)
summary(predsallwb_ds1$num_waterbody_up, na.rm = TRUE)
summary(predsallwb_ds1$dist_nearest_up_wb_km, na.rm = TRUE)
predsallwb_ds1$dist_nearest_up_wb_km <- round(predsallwb_ds1$dist_nearest_up_wb_km, 1)

summary(predsallwb_ds1$dist_nearest_up_wb_km, na.rm = TRUE)

summary(predsallwb_ds1$dist_nearest_up_wb_km)

predsallwb_ds1$dist_nearest_up_wb_km[is.na(predsallwb_ds1$dist_nearest_up_wb_km)] <- 0

summary(predsallwb_ds1$dist_nearest_up_wb_km)

# file below provides a link to waterbody, and subsequently to DFobs data so commented out for pred
# StationIDs_UniqLocIDs <- read.csv("E:/R_vdeq_sci/Working/StationIDs_UniqLocIDs.csv") %>% clean_names(.)
# 
# names(StationIDs_UniqLocIDs)
# 
# wb_ds2 <- left_join(StationIDs_UniqLocIDs, wb_ds1, by=join_by(uniq_loc_id))
# names(wb_ds2)

# from help by = join_by(name == artist)
# station_id_2 and station_id are both long station trend names
preds3 <- left_join(preds3, predsallwb_ds1, by = join_by(station_id_2 == station_id))
names(preds3)

preds3 <- preds3 %>% 
  mutate(
    wbc = case_when(
      incl_nhdwb == "yes" ~ "present",
      incl_nhdwb != "yes" ~ "absent"
    )
  )

preds3$wbc <- factor(preds3$wbc, levels = c("present", "absent"))
summary(preds3$wbc)

preds3 <- preds3 %>% 
  mutate(
    wbc_bin = case_when(
      incl_nhdwb == "yes" ~ 1,
      incl_nhdwb != "yes" ~ 0
    )
  )

preds3$wbc_bin <- as.numeric(preds3$wbc_bin)

names(preds3)

class(preds3)
# note ssn_put_data requires sf object and SSN2 object
j_ssn3 <-  SSN2::ssn_put_data(preds3,j_ssn3, name = "Preds_2019_2022", resize_data = FALSE)
summary(j_ssn3)

# saveRDS(j_ssn3, file = "j_ssn3.rds")

# VARIABLE ADJUSTMENT ZONE 4
### Variables to apply empirical logit transformation
# emplog_vars <- c("pct_for_w","pct_imp_w","pct_crop_w","pct_hay_w","pct_grs_w","pct_shrb_w","pct_for_wr","pct_imp_rp_w","pct_crop_wr","pct_hay_wr","pct_grs_wr","pct_shrb_wr","pct_for_c","pct_imp_c","pct_crop_c","pct_hay_c","pct_grs_c", "pct_shrb_c", "pct_for_cr","pct_imp_rp_c","pct_crop_cr","pct_hay_cr","pct_grs_cr","pct_shrb_cr")

# remove geometry so empirical logit can be applied
predsz <- st_set_geometry(preds3, NULL)
################################################################
################################################################


## transform these variables and put the new values into new columns in preds
for(var in emplog_vars){
  ### create new tranformed data column to preserve the original
  new_nm <- paste0(var,"_emplog")
  dat_vec_obs <- predsz[,var]
  # dat_vec_preds <- DFpreds[,var]
  
  #converting to 0-1 range
  dat_vec_obs <- dat_vec_obs/100
  # dat_vec_preds <- dat_vec_preds/100
  
  # dat_vec[dat_vec == 1] <- .9999
  # dat_vec[dat_vec == 0] <- .0001
  
  if(any(dat_vec_obs > 1 | dat_vec_obs < 0)){
    cat("ERROR: percentage variables outside logical bounds")
  }
  
  small_dat_vec_obs <- dat_vec_obs[dat_vec_obs <1 & dat_vec_obs >0]
  # small_dat_vec_preds <- dat_vec_preds[dat_vec_preds <1 & dat_vec_preds >0]
  op1_obs <- small_dat_vec_obs
  op2_obs<- 1-small_dat_vec_obs
  # op1_preds <- small_dat_vec_preds
  # op2_preds <- 1-small_dat_vec_preds
  
  ## minimum of op1 op2
  delt_obs <- min(c(op1_obs,op2_obs))
  # delt_preds <- min(c(op1_preds,op2_preds))
  
  ## getting set of frequencies
  freqs_obs <- NULL
  for(i in 1:length(dat_vec_obs)){
    if(dat_vec_obs[i] <= delt_obs){
      freqs_obs[i] <- delt_obs/2
    }else if(dat_vec_obs[i] >= 1- delt_obs){
      freqs_obs[i] <- 1-(delt_obs/2)
    }else{
      freqs_obs[i] <-dat_vec_obs[i]
    }
  }
  
 # freqs_preds <- NULL
 # for(i in 1:length(dat_vec_preds)){
 #   if(dat_vec_preds[i] <= delt_preds){
 #     freqs_preds[i] <- delt_preds/2
 #   }else if(dat_vec_preds[i] >= 1- delt_preds){
 #     freqs_preds[i] <- 1-(delt_preds/2)
 #   }else{
 #     freqs_preds[i] <-dat_vec_preds[i]
  #  }
 # }
  
  ##getting logits
  logits_obs <- log(freqs_obs/(1-freqs_obs))
  predsz[,new_nm] <- logits_obs
  
  # logits_preds <- log(freqs_preds/(1-freqs_preds))
  # DFpreds[,new_nm] <- logits_preds
}

predsz <- predsz %>%
  mutate(
    imp_w_pres = case_when(
    pct_imp_w_emplog <= -8.111428 ~ 0,
    pct_imp_w_emplog > -8.111428 ~1
  )
)

head(predsz$imp_w_pres)
head(predsz$pct_imp_w)
names(predsz)
predsz2 <- dplyr::select(predsz, c(station_id_2, pct_for_w_emplog:imp_w_pres))
# put transformed covariates in an SF object
preds3a <- full_join(preds3, predsz2, by = join_by(station_id_2))
summary(preds3a$pct_imp_w_emplog)
class(preds3a)

# saveRDS(preds3a, file = "preds3a.rds")

# put SF object into SSN
j_ssn3 <-  SSN2::ssn_put_data(preds3a,j_ssn3, name = "Preds_2019_2022", resize_data = FALSE)
summary(j_ssn3$preds$Preds_2019_2022$pct_imp_w_emplog)


# dummy code 5 vahusb with base being JU, James Upper
vahusb <- dplyr::select(predsz, 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(predsz$vahusb)
distinct(vahusb_d)
head(vahusb_d)
str(vahusb_d)
class(vahusb_d)

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

P1 <- predsz|>
   dplyr::select(station_id_2, vahusb_JM, vahusb_JR, vahusb_JA, vahusb_JL) %>%
  mutate_at(c('vahusb_JM', 'vahusb_JR', 'vahusb_JA', 'vahusb_JL'), as.numeric)
str(P1)

# P1 can only contain numeric or factor for 22 preds

# put dummy covariates in an SF object
preds3b <- full_join(preds3a, P1, by = join_by(station_id_2))
names(preds3b)

# put SF object into SSN
j_ssn3 <-  SSN2::ssn_put_data(preds3b,j_ssn3, "Preds_2019_2022", resize_data = FALSE)

summary(j_ssn3)

# Check that afv_area is in SSN object
str(j_ssn3$preds$Preds_2019_2022$afv_area)
```

# Create Distance Matrix
When using a new SSN object, such as James_071024_pluspreds.ssn have to create distance matrices. Ran first set of code and now see a distance matrix folder in where the SSN object is stored at:
ssn_object/James_071024_pluspreds.ssn and see sub folders for obs and Preds_2019_2022 subfolder.
Already run so commented out.
```{r distance_matrix}

# Distance matrix for obs and first set of prediction points from 2019-2022
# SSN2::ssn_create_distmat(j_ssn3, predpts = "Preds_2019_2022", overwrite = TRUE)


```


## Torgegrams DO, TN and Partitioning
In the coefficent plot, the flip-flop in instream stressor, being dissolved oxygen at Watershed and Watershed-Riparian, and then total nitrogen for Catchment and Catchment-Riparian suggests that different stressors have their effects on stream condition index at different spatial extents. I thought of following that up in two ways. First, see if semivariograms and Torgegrams of dissolved oxygen and total nitrogen show different ranges, the distance at which observations are independent. My hypothesis is that dissolved oxygen will have a larger range than total nitrogen. Second, I want to see if I update the watershed model to include all 3 instream stressors does that beat the current watershed model having only dissolved oxygen.

On 05/03/2024, not seeing that much spatial structure in DO or Log(TN). On 07/12/2024 still need to test if 3 instream stressors in a model outperforms just DO in the model. Try subsetting ssn objectby ju factor and then making separate Torgegrams.
```{r do_tn_toregrams}
# 
# class(j_ssn3)
# ztg <- SSN2::Torgegram(vscivcpmi ~ 1, j_ssn3, type = c("flowcon", "flowuncon", "euclid"))
# # SSN2::plot.Torgegram(ztg)
# torg <- ztg[["euclid"]]
# names(torg)
# class(torg)
# ggplot(torg, aes(x=dist, y=gamma,size=np)) + geom_point()
# 
# summary(j_ssn3$obs$ju)
# glimpse(j_ssn3$obs$ju)
# 
# #ju yes has 71 obs and no has 128 obs
# # esv_separate <- map(c(yes, no), ~ SSN2::Torgegram(vscivcpmi ~ 1, j_ssn3 |> filter(ju == .x)))
# # 
# # esv_separate
# # 
# # # plot(esv_separate[[1]]) not helpful
# # 
# # esv_ju_yes <- esv_separate[[1]] 
# # esv_ju_no <- esv_separate[[2]]
# # 
# # names(esv_ju_yes)
# # 
# # esv1 <- ggplot(esv_ju_yes, aes(x=dist, y=gamma, size=np)) + geom_point() + labs(x = "Euclidean Distance (m)", y = "Semivariance", title = "JU yes") + scale_x_continuous(labels=comma)
# # 
# # esv2 <- ggplot(esv_ju_no, aes(x=dist, y=gamma, size=np)) + geom_point() + labs(x = "Euclidean Distance (m)", y = "Semivariance", title = "JU no") + scale_x_continuous(labels=comma)
# # 
# 
# 
# 
# ztg <- SSN2::Torgegram(l_tn ~1, j_ssn3, type = c("flowcon", "flowuncon", "euclid"))
# plot(ztg, main = "Torgegram Log(TN)")
# plot(ztg, type = "flowcon", main = "Torgegram Log(TN)")
# plot(ztg, type = "flowuncon", main = "Torgegram Log(TN)")
# plot(ztg, type = "euclid", main = "Torgegram Log(TN)")
# plot(ztg, separate = TRUE, main = "Torgegram Log(TN)")
# 
# flowcon <- ztg[["flowcon"]]
# ggplot(flowcon, aes(x=dist, y=gamma)) + geom_point()
# flowcon$dist_type <- "flow_conn"
# 
# flowuncon <- ztg[["flowuncon"]]
# names(flowuncon)
# class(flowuncon)
# ggplot(flowuncon, aes(x=dist, y=gamma)) + geom_point()
# flowuncon$dist_type <- "flow_unconn"
# 
# euclid <- ztg[["euclid"]]
# names(euclid)
# class(euclid)
# ggplot(euclid, aes(x=dist, y=gamma)) + geom_point()
# euclid$dist_type <- "euclid"
# 
# torg <- bind_rows(flowuncon,euclid)
# ggplot(torg, aes(x=dist, y=gamma, colour = dist_type)) + geom_point()
# 
# # test if Torgegram would work using similar code below
# # esv_separate <- map(c(2001, 2006), ~ esv(log_Zn ~ log_dist2road, moss |> filter(year == .x)))
# 

```


# Ws-Wq Model
This is the SSN Ws-Wq model that was evaluated in script 05. Here is it run again and augmented with diagnostic and fitted values. The parallel coordinate plot of the augmented data frame was grouped by the 5 vahusb. The clustering of the JU variables in that plot is what suggested using the partition factor. See the description of partition factor under section 4.1.3 Advanced Features in An Introduction to Spatial Stream Network Modeling in R using SSN2.

On 04/19/2024 kept tail down exponential and changed Euclidean to exponential based on models_yintercept.csv table.
On 04/10/2023 evaluated model using vsci, untransformed, and added imp_w_pres, a binary presence indicator of impervious cover. This was approach was used by Dumelle et al. 2023 in their spatial modeling of NLA conductivity.
```{r wswq_ssn_obs_fitted}
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)

#plot(ssn_wswq_reml1, which = c(1:6))

tidy(ssn_wswq_reml1)

aug_ssn_wswq_reml1 <- augment(ssn_wswq_reml1, drop = FALSE)
class(aug_ssn_wswq_reml1)
names(aug_ssn_wswq_reml1)

mapview(aug_ssn_wswq_reml1)

mapview(aug_ssn_wswq_reml1, zcol = ".std.resid", cex = ".std.resid", alpha.regions = .8, legend = TRUE, popup = popupTable(aug_ssn_wswq_reml1, zcol = c(".std.resid")))# + mapview(stream)

# parallel coordinate plot
# names(aug_ssn_wswq_reml1)
pcpobs <- ggparcoord(data = aug_ssn_wswq_reml1, columns = c(136,155,38,161,165,204,166,232), groupColumn = "vahusb", scale = "std", showPoint = TRUE, title = "Observed Sites", alphaLines = 0.8, boxplot = FALSE)
pcpobs

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()

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

class(resid_ssn1)
names(resid_ssn1)

ssn_wswq_fit_reml1 <- dplyr::select(aug_ssn_wswq_reml1, c(station_id, station_id_2, year, vahusb, do, tn, tothab, l_tn, vscivcpmi, pct_imp_c, pct_imp_w, elev_ws, pct_imp_w_emplog, vsci, .fitted, .std.resid))

saveRDS(ssn_wswq_fit_reml1, file = "outputs/ssn_wswq_fit_reml1.rds")

st_write(ssn_wswq_fit_reml1, dsn = file.path(getwd(), "ssn_wswq_fit_reml1.gpkg"), layer = "ssn_wswq_fit_reml1", driver = "GPKG", quiet = FALSE, append = FALSE)

# put selected 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)

res_tg1 <- SSN2::Torgegram(.std.resid ~ 1, j_ssn4, type = c("flowuncon", "euclid"))
plot(res_tg1, main = "Torgegram of Standardized Residuals from Watershed SSN")
```

# Ws-Wq Partitioned
The z4 model below has the partitioned factor of ju, Upper James. The models z1-z3 tested different coding and how different fixed and random effects compared to ssn_wswq_reml1. Finally, z4 model had a better fit, based on AIcc, than ssn_wswq_reml1. Also, the RMSPE is slightly smaller for the z4 model, and its cor2 is slightly larger.

On 11/25/2024 why does z4ssn_wswq_ml1 model give much more reasonable range estimates than reml model?
On 07/03/2024 need to look at diagnostic plots of partition model. Diagnostic plots are saved as Word files at E:\R_vdeq_sci\Working\Data\neptune_analysis\scripts_by_basin.
On 06/28/2024, I coded for a partition factor base on JU subbasin versus non-JU subasins , partly because parallel coordinate plot showed JU to differ in covariates and vsci from rest of subbasins. Results from dummy coded model above match factor coded model below. Also, on 06/20/2024 compared the Watershed ssn model to a model with a random effect of year. The Watershed ssn model had a much lower AICc than the model with a random effect of year.
```{r wswq_ssn_partition}
# z1 model shows that jl as factor gives same results at dummy coded vahusb_JL
z1ssn_wswq_reml1 <- ssn_lm(
  formula = vsci ~ pct_imp_w_emplog + elev_ws + do + jl,
  ssn.object = j_ssn3,
  tailup_type = "none",
  taildown_type = "exponential",
  euclid_type = "exponential",
  nugget_type = "nugget",
  estmethod = "reml",
  additive = "afv_area"
)
summary(z1ssn_wswq_reml1)

# random effect of year
z2ssn_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",
   random = ~ year_f
)

summary(z2ssn_wswq_reml1)
varcomp(z2ssn_wswq_reml1)
loocv(z2ssn_wswq_reml1)

#plot(z2ssn_wswq_reml1, which = c(1:6))
# 
glances(ssn_wswq_reml1,z2ssn_wswq_reml1)
# AICc ssn_ws_reml1 1415 vs. AICc z2ssn_wswq_reml1 1424

# fixed effect of year
z3ssn_wswq_reml1 <- ssn_lm(
   formula = vsci ~ pct_imp_w_emplog + elev_ws + do + vahusb_JL + year_f,
   ssn.object = j_ssn3,
   tailup_type = "none",
   taildown_type = "exponential",
   euclid_type = "exponential",
   nugget_type = "nugget",
   estmethod = "reml",
   additive = "afv_area"
 )

summary(z3ssn_wswq_reml1)
varcomp(z3ssn_wswq_reml1)
loocv(z3ssn_wswq_reml1)

#plot(z3ssn_wswq_reml1, which = c(1:6))

glances(ssn_wswq_reml1,z2ssn_wswq_reml1,z3ssn_wswq_reml1)
# AICc ssn_ws_reml1 1415 vs. AICc z2ssn_wswq_reml1 1423 vs AICc z3ssn 1327
# RMSPE/cor2 ssn_ws_reml1 8/38/0.41 vs. RMSPE/cor2 z2ssn_wswq_reml1 8.46/0.40 vs RMSPE/cor2 z3ssn 8.70/0.378

# On 06/28/2024 partition factor on ju motivated by parallel coordinate plots of obs and preds.  Partition factor, from help, affects covariance portion of the model while fixed effects affect the mean portion. Therefore I think I can compare 2 models based on AICc derived from reml. 

z4ssn_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",
  partition_factor = ~ ju
)
summary(z4ssn_wswq_reml1)
varcomp(z4ssn_wswq_reml1)
loocv(z4ssn_wswq_reml1) %>% print(n = Inf)

glances(ssn_wswq_reml1,z4ssn_wswq_reml1)

#plot(z4ssn_wswq_reml1, which = c(1:6))

tidy(z4ssn_wswq_reml1)

aug_z4_wswq_reml1 <- augment(z4ssn_wswq_reml1, drop = FALSE)
summary(aug_z4_wswq_reml1$.std.resid)
mapview(aug_z4_wswq_reml1, zcol = ".std.resid", cex = ".std.resid", alpha.regions = .8, legend = TRUE, popup = popupTable(aug_ssn_wswq_reml1, zcol = c(".std.resid"))) 

aug_z4_wswq_reml1 <- dplyr::select(aug_z4_wswq_reml1, c(station_id, station_id_2, year, vahusb, do, tn, tothab, l_tn, vscivcpmi, pct_imp_c, pct_imp_w, elev_ws, pct_imp_w_emplog, vsci, .fitted, .std.resid))

saveRDS(aug_z4_wswq_reml1, file = "outputs/aug_z4_wswq_reml1.rds")

st_write(aug_z4_wswq_reml1, dsn = file.path(getwd(), "aug_z4_wswq_reml1.gpkg"), layer = "aug_z4_wswq_reml1", driver = "GPKG", quiet = FALSE, append=FALSE)

# B/c want to compare predictions from models having different covariates specifying estimation method as ml

z4ssn_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",
  partition_factor = ~ ju
)
summary(z4ssn_wswq_ml1)
varcomp(z4ssn_wswq_ml1)
loocv(z4ssn_wswq_ml1) %>% print(n = Inf)



```

# TLH fitted to obs interval plot
On 12/23/2024 commenting this code out b/c this is not the model I need for fitted interval obs plot.

On 12/16/2024 first using a4 reml model and will compare to z5 reml model with distance covariate as that provides a better fit based on smaller AICc value.
```{r fitted_interval_obs_plot}
# z4 model was run above

## fit model
# z4ssn_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",
#   partition_factor = ~ju
# )
# summary(z4ssn_wswq_reml1)

# how do residuals and fitted values below compare to what comes out with augment function?



# class(aug_z4_wswq_reml1)
# names(aug_z4_wswq_reml1)
# 
# head(round(aug_z4_wswq_reml1$.fitted))
# head(round(aug_z4_wswq_reml1$vsci))


# get spatial residuals from model
# spatResids<-residuals(z4ssn_wswq_reml1)
# zspatResids<-residuals(z4ssn_wswq_reml1)
# head(aug_z4ssn_wswq_reml1$.resid)
# the three statements above all return same results

# pull sd of residuals and fitted values
# resid_SD = sd(spatResids)
# fitted_vals = z4ssn_wswq_reml1$fitted$response

# combine to a dataframe
# spat<-cbind(fitted_vals,resid_SD^2)
# head(spat)

# spatial prediction percentiles
# (spat.Pred.pct<-t(apply(spat,1,function(x) qnorm(c(Q05=0.05,Q10=0.1,Q25=0.25,Median=0.5,Q75=0.75,Q90=0.9,Q95=0.95),x[1],sqrt(x[2])))))

# spat.Pred.median = spat.Pred.pct[,4]
# head(spat.Pred.median)

# make sure I've got the right obj
# DFobs<-ssn_get_data(j_ssn3)
# observed and median of our prediction intervals
# ty<-DFobs$vsci
# tx<-spat.Pred.median

#percentiles
# ox<-order(tx)
# op05<-spat.Pred.pct[ox,"Q05"]
# op95 <- spat.Pred.pct[ox,"Q95"]

# plot
# plot(tx,ty,xlim=range(na.omit(tx)),ylim=range(na.omit(ty),na.omit(op95)),xlab="VSCI Predicted Value",ylab="VSCI Observed Value",main="Spatial Model Observed vs. Predicted z4 reml n = 199")
# abline(0,1,col="red")
# lines(lowess(na.omit(tx),na.omit(ty)),col='blue')

# 5th percentile of predictive distribution
# lines(tx[ox],op05,lty=2,col='black')
# 95th percentile of predictive distribution
# lines(tx[ox],op95,lty=2,col='black')
#legend
# legend("topleft",lty=c(1,1,2,2),col=c("red","blue","black","black"),legend=c("y = x","lowess","5% Pred","95% Pred"), x.intersp = 0.7, y.intersp = 0.7)

# cor(tx,ty, method = "pearson")

# these object names are used again downstream for preds so remove here

# rm(spatResids,resid_SD,fitted_vals,spat,spat.Pred.pct,spat.Pred.median,ty,tx,ox,op05,op95)
```


# Additional Hypotheses for a4 model
On 12/16/2024 in near future will likely comment out this code chunk as no support found for including these additional covariates, which are included in model a4ssn_wswq_ml1.

On 11/25/2024 From Miller et al. 2017's results, suggested hypotheses about WWTPs and upstream reservoirs. Also, I will test a depth to bedrock soils hypothesis. The goal of adding these covariates is to get improvements in predictions for 2019-2022 data. This is referred to as the a4 model, for additional hypotheses. Note because of new covariates using ml estimation (Ver Hoef 2014 page 23). Commented out SSN diagnostic plot which statement as seemed cause error with following residual assignment.
```{r additional_hypotheses}

a4ssn_wswq_ml1 <- ssn_lm(
  formula = vsci ~ pct_imp_w_emplog + elev_ws + do + vahusb_JL + rck_dep_ws + wwtp + wbc,
  ssn.object = j_ssn3,
  tailup_type = "none",
  taildown_type = "exponential",
  euclid_type = "exponential",
  nugget_type = "nugget",
  estmethod = "ml",
  additive = "afv_area",
  partition_factor = ~ju
)
summary(a4ssn_wswq_ml1)
varcomp(a4ssn_wswq_ml1)
loocv(a4ssn_wswq_ml1) %>% print(n = Inf)

glances(z4ssn_wswq_ml1,a4ssn_wswq_ml1)
loocv(z4ssn_wswq_ml1)

# plot(a4ssn_wswq_ml1, which = c(1:6))

# tidy(a4ssn_wswq_reml1)

# inserted r in file name to indicate residuals now included
a4rssn_wswq_ml1 <- augment(a4ssn_wswq_ml1, drop = FALSE)
# note on 10/04/2024 fixed mapview so popupTable using correct object, had been pointing to ssn_wswq_reml1
summary(a4rssn_wswq_ml1$.std.resid)
mapview(a4rssn_wswq_ml1, zcol = ".std.resid", cex = ".std.resid", alpha.regions = .8, legend = TRUE, popup = popupTable(a4rssn_wswq_ml1, zcol = c(".std.resid"))) 

a4rssn_wswq_ml1 <- dplyr::select(a4rssn_wswq_ml1, c(station_id, station_id_2, year, vahusb, do, tn, tothab, l_tn, vscivcpmi, pct_imp_c, pct_imp_w, elev_ws, pct_imp_w_emplog, vsci, .fitted, .std.resid))

# saveRDS(a4rssn_wswq_ml1, file = "a4rssn_wswq_ml1.rds")
# 
# st_write(a4rssn_wswq_ml1, dsn = file.path(getwd(), "a4rssn_wswq_ml1.gpkg"), layer = "a4rssn_wswq_ml1", driver = "GPKG", quiet = FALSE)

```

# Distance Hypothesis
On 12/16/2024 testing distance hypothesis of upstream waterbodies by now using wbc_bin and 
dist_nearest_up_wb_km with those covariates entered into the model as 
wbc_bin
wbc_bin:dist_nearest_up_wb_km.
This follows coding as used in Dumelle et al. 2023 paper. Using that coding produced a poorer fitting more complicated model than z4. Using just the continuous covariate of distance, as done on 11/19/2024, produces a better fitting model a5ssn_wswq_ml1. Model a5ssn_wswq_ml1 has slightly better values for RMSPE and cor2 than z4ssn_wswq_ml1. What does added variable plot of that model look like?


Run ML estimation when comparing to models with different covariates. Run with REML when comparing to models with different spatial autocovariances and for prediction.
On 11/19/2024 Tried running interaction of wbc*dist_nearest_up_wb_km, but got error of singular matrix. Consequently, only tested model using continuous covariate dist_nearest_up_wb_km.
Again, commented out plot(a5ssn_wswq_ml1, which = c(1:6)) as seems to prevent code running straight through. The a5ssn_wswq_ml1 model with only continuous covariate dist_nearest_up_wb_km has smaller AICc than z4 model.
```{r distance_hypothesis}
# test distance covariate in model a5

ggplot(DFobs3a, aes(x = dist_nearest_up_wb_km, y = vsci)) + geom_point() + facet_wrap(vars(wbc))  + geom_smooth(method = "lm", se =TRUE)

ggplot(DFobs3a, aes(x = dist_nearest_up_wb_km, y = vsci)) + geom_point() + geom_smooth(method = "lm", se =TRUE)

a5ssn_wswq_ml1 <- ssn_lm(
  formula = vsci ~ pct_imp_w_emplog + elev_ws + do + vahusb_JL + dist_nearest_up_wb_km,
  ssn.object = j_ssn3,
  tailup_type = "none",
  taildown_type = "exponential",
  euclid_type = "exponential",
  nugget_type = "nugget",
  estmethod = "ml",
  additive = "afv_area",
  partition_factor = ~ju
)

summary(a5ssn_wswq_ml1)
varcomp(a5ssn_wswq_ml1)
loocv(a5ssn_wswq_ml1)

glances(z4ssn_wswq_ml1,a5ssn_wswq_ml1)

# plot(a5ssn_wswq_ml1, which = c(1:6))


#### TLH lm discussion


# test = ssn_get_data(j_ssn3)
# fit = lm(vsci ~pct_imp_w_emplog + elev_ws + do + vahusb_JL + wbc*dist_nearest_up_wb_km - dist_nearest_up_wb_km,data = test)
# summary(fit)
# check_model(fit)
# avPlots(fit)
### as you can see we get a nice output with wbcabsent:dist_nearest_up_wb_km NA'd out. So this is what I'm trying to replicate in the ssn2 formulation. 

# I added test1 code with only distance covariate

# test = ssn_get_data(j_ssn3)
# fit1 = lm(vsci ~pct_imp_w_emplog + elev_ws + do + vahusb_JL + dist_nearest_up_wb_km,data = test)
# summary(fit1)
# check_model(fit1)
# avPlots(fit1)
```

# Check Nonspatial Diagnostics on a5 model
Get checks on collinearity, added variable plots and other diagnostics.
Diagnostics look good.
```{r a5_mlr_diagnostics}
a5_mlr1 <- lm(vsci ~ pct_imp_w_emplog + elev_ws + do + vahusb_JL + dist_nearest_up_wb_km,
  data = DFobs3a)

summary(a5_mlr1)
check_model(a5_mlr1)
avPlots(a5_mlr1)
```

# TLH fitted to obs interval plot for a5 model

Now that diagnostics checked out make an observed to fitted interval plot on model a5ssn_wswq_reml1 by augmenting to get diagonstics.
```{r fitted_interval_obs_plot}
# z4 model was run above

a5ssn_wswq_reml1 <- ssn_lm(
  formula = vsci ~ pct_imp_w_emplog + elev_ws + do + vahusb_JL + dist_nearest_up_wb_km,
  ssn.object = j_ssn3,
  tailup_type = "none",
  taildown_type = "exponential",
  euclid_type = "exponential",
  nugget_type = "nugget",
  estmethod = "reml",
  additive = "afv_area",
  partition_factor = ~ju
)

summary(a5ssn_wswq_reml1)
varcomp(a5ssn_wswq_reml1)
loocv(a5ssn_wswq_reml1)

# how do residuals and fitted values below compare to what comes out with augment function?

aug_a5obs1ws_vsci <- augment(a5ssn_wswq_reml1, drop = FALSE)

class(aug_a5obs1ws_vsci)
names(aug_a5obs1ws_vsci)

head(round(aug_a5obs1ws_vsci$.fitted))
head(round(aug_a5obs1ws_vsci$vsci))


# get spatial residuals from model
spatResids<-residuals(a5ssn_wswq_reml1)
# zspatResids<-residuals(a5ssn_wswq_reml1)
# head(aug_a5obs1ws_vsci$.resid)
# the three statements above all return same results

# pull sd of residuals and fitted values
resid_SD = sd(spatResids)
fitted_vals = a5ssn_wswq_reml1$fitted$response

# combine to a dataframe
spat<-cbind(fitted_vals,resid_SD^2)
head(spat)

# spatial prediction percentiles
(spat.Pred.pct<-t(apply(spat,1,function(x) qnorm(c(Q05=0.05,Q10=0.1,Q25=0.25,Median=0.5,Q75=0.75,Q90=0.9,Q95=0.95),x[1],sqrt(x[2])))))

spat.Pred.median = spat.Pred.pct[,4]
head(spat.Pred.median)

# make sure I've got the right obj
DFobs<-ssn_get_data(j_ssn3)
# observed and median of our prediction intervals
ty<-DFobs$vsci
tx<-spat.Pred.median

#percentiles
ox<-order(tx)
op05<-spat.Pred.pct[ox,"Q05"]
op95 <- spat.Pred.pct[ox,"Q95"]

# plot
plot(tx,ty,xlim=range(na.omit(tx)),ylim=range(na.omit(ty),na.omit(op95)),xlab="VSCI Predicted Value",ylab="VSCI Observed Value",main="Spatial Model Observed vs. Predicted a5 reml n = 199")
abline(0,1,col="red")
lines(lowess(na.omit(tx),na.omit(ty)),col='blue')

# 5th percentile of predictive distribution
lines(tx[ox],op05,lty=2,col='black')
# 95th percentile of predictive distribution
lines(tx[ox],op95,lty=2,col='black')
#legend
legend("topleft",lty=c(1,1,2,2),col=c("red","blue","black","black"),legend=c("y = x","lowess","5% Pred","95% Pred"), x.intersp = 0.7, y.intersp = 0.7)

cor(tx,ty, method = "pearson")

# these object names are used again downstream for preds so remove here

rm(spatResids,resid_SD,fitted_vals,spat,spat.Pred.pct,spat.Pred.median,ty,tx,ox,op05,op95)

```



# Preds for Distance Hypothesis
Note for predictions want estimation method to be reml.
Looking to see if distance to nearest upstream waterbody covariate pulls outliers closer to 1:1 line, which is the a5 model.
```{r distance_hypothesis_preds}
# check stats on pred covariates to make sure corrects preds ssn data frame being used
# summary(j_ssn3$preds$Preds_2019_2022$pct_imp_w_emplog)
# summary(j_ssn3$preds$Preds_2019_2022$elev_ws)
# summary(j_ssn3$preds$Preds_2019_2022$do)
# summary(j_ssn3$preds$Preds_2019_2022$vahusb_JL)
# summary(j_ssn3$preds$Preds_2019_2022$dist_nearest_up_wb_km)
# 
# 
# summary(j_ssn3)
# str(j_ssn3$obs$pct_imp_w_emplog)

a5ssn_wswq_reml1 <- ssn_lm(
  formula = vsci ~ pct_imp_w_emplog + elev_ws + do + vahusb_JL + dist_nearest_up_wb_km,
  ssn.object = j_ssn3,
  tailup_type = "none",
  taildown_type = "exponential",
  euclid_type = "exponential",
  nugget_type = "nugget",
  estmethod = "reml",
  additive = "afv_area",
  partition_factor = ~ju
)
summary(a5ssn_wswq_reml1)
varcomp(a5ssn_wswq_reml1)
loocv(a5ssn_wswq_reml1)

# augment returns an sf object
aug_a5predict1ws_vsci <- augment(a5ssn_wswq_reml1, newdata = "Preds_2019_2022", se_fit = TRUE, interval = c("prediction"))
# class(predict1ws_vsci)
# names(predict1ws_vsci)
aug_a5predict1ws_vsci$vsci_pred <- round(aug_a5predict1ws_vsci$.fitted,1)
aug_a5predict1ws_vsci$vsci_lwr <- round(aug_a5predict1ws_vsci$.lower,1)
aug_a5predict1ws_vsci$vsci_upr <- round(aug_a5predict1ws_vsci$.upper,1)

# Pineiro 2008 put observed on y and predict on x
# label specific points
# https://stackoverflow.com/questions/15624656/label-points-in-geom-point

# plotting without 1:1 abline shows how weak correlation is.
ggplot(aug_a5predict1ws_vsci, aes(x = vsci_pred, y = vsci)) + geom_point() + geom_abline()

a5preds2019_2022_v1 <- ggplot(aug_a5predict1ws_vsci, aes(x = vsci_pred, y = vsci, label = station_id_2)) + geom_point() + geom_abline() + geom_text(aes(label =ifelse(vsci<45, as.character(station_id_2), '')), hjust = -0.1, vjust = 0)+ labs(x = "Predicted VSCI", y = "Observed VSCI, 2019-2022", title = "a5reml Watershed SSN JU Partitioned & Upstream Waterbody Distance") + xlim(35,80) + ylim(35,80)

png(file="figures_sfs/a5preds2019_2022_v1.png",width=6,height=3,units="in",res=150)
a5preds2019_2022_v1
dev.off()

cor(aug_a5predict1ws_vsci$vsci_pred,aug_a5predict1ws_vsci$vsci, method = "pearson")

#2019-2022 a5reml correlation with partition factor of ju and distance covariate is 0.57

as.data.frame(aug_a5predict1ws_vsci)|>
  dplyr::select(vsci_pred,vsci)|>
  correlate()

#check as outlier differs in 2019-20, 2-BSC001.58, vs 2021-22, 2-WPK002.93
# as.data.frame(aug_a5predict1ws_vsci)|>
#   dplyr::filter(station_id_2 != "2-WPK002.93" & station_id_2 != "2-BSC001.58") |>
#   dplyr::select(vsci_pred,vsci)|>
#   correlate()
# dropping outliers a5reml gives corr 0.49 for 2019-2022

# ggplot(aug_a5predict1ws_vsci, aes(x = vsci_pred, y = vsci)) + geom_point() + geom_abline() + geom_errorbar(aes(xmax = vsci_upr, xmin = vsci_lwr)) + labs(x = "Predicted VSCI", y = "Observed VSCI, 2019-2022")
# 
# mapview(aug_a5predict1ws_vsci)

# gather covariates, preds, write out geopackage for 2 prediction datasets. Comment out code based if preds1 or preds2 is used

# a4preds_cov_19_20 <- dplyr::select(a4predict1ws_vsci, c(station_id, st_id_tren, year, vahusb, do, tn, tothab, l_tn, vscivcpmi, pct_imp_c, elev_ws, pct_imp_w_emplog, vsci, vahusb, vsci_pred, .se.fit, vsci_lwr,vsci_upr))
# 
# saveRDS(a4preds_cov_19_20, file = "a4preds_cov_19_20.rds")
# 
# st_write(a4preds_cov_19_20, dsn = file.path(getwd(), "a4preds_cov_1920.gpkg"), layer = "a4preds_cov_1920", driver = "GPKG", quiet = FALSE)

# a4preds_cov_19_22 <- dplyr::select(a4predict1ws_vsci, c(station_id, st_id_tren, year, vahusb, do, tn, tothab, l_tn, vscivcpmi, pct_imp_c, elev_ws, pct_imp_w_emplog, vahusb, vsci, vsci_pred, .se.fit, vsci_lwr,vsci_upr))
# 
# saveRDS(a4preds_cov_19_22, file = "a4preds_cov_19_22.rds")
# 
# st_write(a4preds_cov_19_22, dsn = file.path(getwd(), "preds_cov_19_22.gpkg"), layer = "preds_cov_1922", driver = "GPKG", quiet = FALSE)


```

### TLH Fitted interval preds
On 12/23/2024, merged Travis' code so that predictions from 48 pred sites also used residuals from the 48 predictions to build the fitted interval.
On 12/19/2025 removed a5ssn_wswq_reml1 and replaced with aug_a5predict1ws_vsci, with that SSN model based on the 48 pred sites.
I believe this coding is fitting to 199 obs because using a5ssn_wswq_reml1. I think it should be using aug_a5predict1ws_vsci as it has 48 predicted sites from 2019-2022.
```{r fitted_interval_preds_plot}


# create residuals in the new prediction frame as they don't provide from the prediction augment call
aug_a5predict1ws_vsci$.residuals <- aug_a5predict1ws_vsci$vsci-aug_a5predict1ws_vsci$.fitted

# pull sd of residuals and fitted values
resid_SD = sd(aug_a5predict1ws_vsci$.residuals)
fitted_vals = aug_a5predict1ws_vsci$.fitted


# combine to a dataframe
spat<-cbind(fitted_vals,resid_SD^2)
head(spat)
# spatial prediction percentiles
(spat.Pred.pct<-t(apply(spat,1,function(x) qnorm(c(Q05=0.05,Q10=0.1,Q25=0.25,Median=0.5,Q75=0.75,Q90=0.9,Q95=0.95),x[1],sqrt(x[2])))))

spat.Pred.median = spat.Pred.pct[,4]

# observed and median of our prediction intervals
ty<-aug_a5predict1ws_vsci$vsci
tx<-spat.Pred.median

#percentiles
ox<-order(tx)
op05<-spat.Pred.pct[ox,"Q05"]
op95 <- spat.Pred.pct[ox,"Q95"]

# plot
plot(tx,ty,xlim=range(na.omit(tx)),ylim=range(na.omit(ty),na.omit(op95)),xlab="VSCI Predicted Value 2019-2022",ylab="VSCI Observed Value 2019-2022",main="SSN Obs vs. Pred a5 reml JU Partitioned & Upstream Waterbody Distance n =48")
abline(0,1,col="red")
lines(lowess(na.omit(tx),na.omit(ty)),col='blue')

# 5th percentile of predictive distribution
lines(tx[ox],op05,lty=2,col='black')
# 95th percentile of predictive distribution
lines(tx[ox],op95,lty=2,col='black')
#legend
legend("topleft",lty=c(1,1,2,2),col=c("red","blue","black","black"),legend=c("y = x","lowess","5% Pred","95% Pred"), x.intersp = 0.7, y.intersp = 0.7)

cor(tx,ty, method = "pearson")

```