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.

01/24/2025 coded for status and trend station types so to evaluate ST_rand model

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(plotly)
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 01/24/2025 added code for st_type, status or trend station type.

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}
# use code to evaluate vsci by year for trend sites if a random effect for trend stations is needed
DFobs <- DFobs %>%
  mutate(
    st_type = case_when(
      station_id_2 == "2-JKS028.69_2004" ~ "status",
      station_id_2 != "2-JKS028.69_2004" & station_id_2 != station_id ~ "trend",
      station_id_2 != "2-JKS028.69_2004" & station_id_2 == station_id ~ "status",
    )
  )
summary(as.factor(DFobs$st_type))

ggplot(DFobs, aes(x = st_type, y = vscivcpmi)) + geom_boxplot() + labs(title = "Obs from 2001-2018 (n=199)")

select(DFobs, station_id_2, st_type, vscivcpmi) %>% arrange(st_type, vscivcpmi, station_id_2) %>% print(n =Inf)

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 01/24/2025 added code for st_type, status or trend station type.

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 <- preds3 %>%
  mutate(
    st_type = case_when(
    station_id ==  station_id_2 ~ "status",
    station_id != station_id_2 ~ "trend"
  )
)

summary(as.factor(preds3$st_type))

ggplot(preds3, aes(x = st_type, y = vscivcpmi)) + geom_boxplot() + labs(title = "Preds 2019-2022 (n=48)")

select(preds3, station_id_2, st_type, vscivcpmi) %>% arrange(st_type, station_id_2, vscivcpmi) %>% print(n =Inf)


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 12/26/2024 commented out models z1-z3

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

# TLH fitted to obs interval plot (LOOCV method)
On 12/26/2024 commenting this code because this is not the model I need for fitted interval obs plot.
```{r fitted_interval_plot_loocv}
# Alrighty, so this method mimics John's method
# which gives a loocv estimate of fitted value an se.
# This results in the squiggly lines because he's giving that 
# se at each individual value.
# cv_resids = loocv(z4ssn_wswq_reml1,
 #                 cv_predict = T,
#                  se.fit = T)
# pull sd of residuals and fitted values
# resid_Se = cv_resids$se.fit
# fitted_vals = cv_resids$cv_predict

# combine to a dataframe
# spat<-cbind(fitted_vals,resid_Se^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="LOOCV Spatial Model Observed vs. Predicted z4 reml")
# 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)

```

# Carson Code Points Outside LOOCV 90% Prediction Limit for n=199 Obs
Added this code on 12/24/2024 and commented out because this is not the final SSN model.
This code comes from local_SSN_FINAL_20160126_modelling.html 
```{r loocv_pts_outside}
# number of points outside 90% prediction limits
# os90pl.sp <-(ty < spat.Pred.pct[,"Q05"]) |  (ty > spat.Pred.pct[,"Q95"]) 
# summary(os90pl.sp)

# which observations are outside [5%, 95%] prediction limits?
# these are the Site IDs of the unusual values
# DFobsz[which(os90pl.sp),"station_id_2"]
```


# 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/26/2024 used geom_ribbon to make prediction interval.

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


# ST_rand Model vs Distance Model
On 01/24/2025 inserted this code chunk to run status/trend random effects model, ssn_wswq_rand1_st, and compare it to distance model, with both models estimated by ml given different covariates.
```{r st_vs_distance}
ssn_wswq_rand1_st <- ssn_lm(
  formula = vsci ~ pct_imp_w_emplog + elev_ws + do + vahusb_JL + dist_nearest_up_wb_km + st_type,
  ssn.object = j_ssn3,
  tailup_type = "none",
  taildown_type = "exponential",
  euclid_type = "exponential",
  nugget_type = "nugget",
  estmethod = "ml",
  additive = "afv_area",
  partition_factor = ~ju,
  random = ~ (1 | station_id)
)

summary(ssn_wswq_rand1_st)
varcomp(ssn_wswq_rand1_st)
loocv(ssn_wswq_rand1_st)

plot(ssn_wswq_rand1_st, which = c(1:6))

glances(a5ssn_wswq_ml1, ssn_wswq_rand1_st)


```

# TLH fitted to obs interval plot for a5 model
01/14/2025 added Travis' code to make prediction of 199 sites vsci versus observed vsci as ggplot using same variance for each obs so get straight 90th percentile lines. 18 Stations fell outside of those lines.

Now that diagnostics checked out make an observed to fitted interval plot on model a5ssn_wswq_reml1 by augmenting to get diagnostics.
```{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,1))
head(round(aug_a5obs1ws_vsci$vsci,1))


# 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 out fitted values and sd of residuals
fitted_vals = a5ssn_wswq_reml1$fitted$response
resid_SD = sd(spatResids)


# combine to a dataframe note variance same for all obs
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)

# observed and median of our fitted values from 199 obs
ty<-aug_a5obs1ws_vsci$vsci
tx<-round(spat.Pred.median,1)

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

# put all components into data frame for ggplot
plot_dat = data.frame("vsci_pred" = tx, "vsci" = ty, "pred_5th" = op05, "pred_95th" = op95,"vsci_bound_order"=tx[ox],"station_id"=aug_a5obs1ws_vsci$station_id_2)


p = ggplot(data = plot_dat)+
  geom_point(aes(x = vsci_pred, y = vsci,label = station_id), color = "black")+
  geom_smooth(aes(x = vsci_pred, y = vsci), method = "loess", se = F, color = "blue")+
  geom_line(aes(x = vsci_bound_order, y = pred_5th), linetype = "dashed", color = "black")+
  geom_line(aes(x = vsci_bound_order, y = pred_95th), linetype = "dashed", color = "black")+
  geom_abline(intercept = 0, slope = 1, color = "red")+
  labs(x="vsci predicted",
       y ="vsci observed")
print(p)
ggplotly(p)

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

# number of points outside 90% prediction limits
os90pl.sp <-(ty < spat.Pred.pct[,"Q05"]) |  (ty > spat.Pred.pct[,"Q95"]) 
summary(os90pl.sp)

# which observations are outside [5%, 95%] prediction limits?
# these are the Site IDs of the unusual values
aug_a5obs1ws_vsci[which(os90pl.sp),"station_id_2"]

mapview(aug_a5obs1ws_vsci[which(os90pl.sp),"station_id_2"])

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

# plot using base plot with all components called separately into the 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="SSN Obs vs. Pred 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)

```

# TLH fitted to obs interval plot (LOOCV method)
On 01/14/2025 inserted TLH's code using observation-specific LOOCV variance with qnorm so get squiggly 90th percentile lines in ggplot.
Returns 19 obs sites falling outside LOOCV 90th percentile lines, and sorting by station_id_2 in mapview label highlight Tye River sites.
Uses se.fit from loocv with qnorm to make squiggly lines for 90th percentile.
On 12/16/2024 insert TLH's code so applied to a5ssn_wswq_reml1 model on 199 obs.
```{r fitted_interval_plot_loocv}
# Alrighty, so this method mimics John's method
# which gives a loocv estimate of fitted value an se.
# This results in the squiggly lines because he's giving that 
# se at each individual value.

cv_resids = loocv(a5ssn_wswq_reml1,
                  cv_predict = T,
                  se.fit = T)
# pull out fitted values and resids se
fitted_vals = cv_resids$cv_predict
resid_Se = cv_resids$se.fit


# combine to a dataframe and se.fit/variance differs for each observation
spat<-cbind(fitted_vals,resid_Se^2)
head(spat)

# spatial prediction percentiles
# qnorm arguments include
# p a vector of probabilities
# mean a vector of means, which are our fitted values
# sd a vector of standard deviations, the square root of the squared .se.fit
(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)

# observed and median of our prediction intervals
ty<-aug_a5obs1ws_vsci$vsci
tx<-round(spat.Pred.median,1)

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

# put all components into data frame for ggplot
plot_dat = data.frame("vsci_pred" = tx, "vsci" = ty, "pred_5th" = op05, "pred_95th" = op95,"vsci_bound_order"=tx[ox],"station_id"=aug_a5obs1ws_vsci$station_id_2)


p = ggplot(data = plot_dat)+
  geom_point(aes(x = vsci_pred, y = vsci,label = station_id), color = "black")+
  geom_smooth(aes(x = vsci_pred, y = vsci), method = "loess", se = F, color = "blue")+
  geom_line(aes(x = vsci_bound_order, y = pred_5th), linetype = "dashed", color = "black")+
  geom_line(aes(x = vsci_bound_order, y = pred_95th), linetype = "dashed", color = "black")+
  geom_abline(intercept = 0, slope = 1, color = "red")+
  labs(x="vsci predicted",
       y ="vsci observed")

print(p)
ggplotly(p)

# Carson Code Points Outside LOOCV 90% Prediction Limit for 199 obs. This code comes from local_SSN_FINAL_20160126_modelling.html 

# number of points outside 90% prediction limits
os90pl.sp <-(ty < spat.Pred.pct[,"Q05"]) |  (ty > spat.Pred.pct[,"Q95"]) 
summary(os90pl.sp)

# which observations are outside [5%, 95%] prediction limits?
# these are the Site IDs of the unusual values
aug_a5obs1ws_vsci[which(os90pl.sp),"station_id_2"]

mapview(aug_a5obs1ws_vsci[which(os90pl.sp),"station_id_2"])

rm(cv_resids, resid_Se, fitted_vals, spat, spat.Pred.pct, spat.Pred.median, ty, tx, ox, op05, op95, os90pl.sp, p)

# plot using base plot pulling in each component separately
# 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="LOOCV SSN Obs vs. Preds 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)

```

# Ribbon LOOCV Fitted Interval
From cv_resids, use fitted_vals and resid_Se to make ribbon for LOOCV fitted interval so to identify obs having interesting vsci values. Uses +/- 1.96*se.fit from loocv to make squiggly ribbon.
```{r ribbon_loocv_fitted_interval}

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

cv_resids = loocv(a5ssn_wswq_reml1,
                  cv_predict = T,
                  se.fit = T)

spat_loocv1 <- cbind(cv_resids$cv_predict, cv_resids$se.fit, aug_a5obs1ws_vsci$vsci)
colnames(spat_loocv1) <- c("cv_predict", "se.fit", "vsci")

spat_loocv1 <- as.data.frame(spat_loocv1)
head(spat_loocv1)

spat_loocv1$cv_predict <- round(spat_loocv1$cv_predict,1)

station_id_2 <-dplyr::select(aug_a5obs1ws_vsci,station_id_2)

# put sf object first so sticky geometry retained
spat_loocv2 <- bind_cols(station_id_2, spat_loocv1)


p <- ggplot(spat_loocv2, aes(x = cv_predict, y = vsci, label = station_id_2)) + geom_point() + geom_abline() + geom_ribbon(aes(ymin = cv_predict - (1.96*se.fit), ymax = cv_predict + (1.96*se.fit)), fill = "grey", alpha = .8) + geom_line(aes( y = cv_predict)) + labs(x = "VSCI LOOCV Fitted 2001-2018", y = "VSCI Observed 2001-2018", title = "SSN LVOOCV Obs vs Fitted a5 reml n=199")

print(p)
ggplotly(p)

# number of points outside 90% loocv fitted limits (ofl)
ofl_stations1 <- filter(spat_loocv2, vsci < (cv_predict - (1.96*se.fit)) | vsci > (cv_predict + (1.96*se.fit)))
ofl_stations1

mapview(ofl_stations1)

edges <- st_read("ssn_object/James_071024_pluspreds.ssn/edges.shp")

# now see 2 trend stations on Tye River
# when zoom in
# https://stackoverflow.com/questions/36469379/multiple-markers-on-same-coordinate
st_jitter(ofl_stations1, factor = 0.001) %>% 
  mapview(.) + mapview(edges)

rm(cv_resids, spat_loocv1, spat_loocv2, p, ofl_stations1)
```

# Preds for Distance Hypothesis
On 01/23/2025 the abrupt change from saw tooth prediction interval to
steady prediction interval is driven by the shorter range of vsci obs and vsci predictions for trend stations compared to the status status stations.
On 01/21/2025 added station type, st_type, of status or trend station as noticed first 11 stations with smallest .se.fit were trend stations.

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
# note level set to 0.90
aug_a5predict1ws_vsci <- augment(a5ssn_wswq_reml1, newdata = "Preds_2019_2022", se_fit = TRUE, interval = c("prediction"), level = 0.90)
# class(aug_a5predict1ws_vsci)
# names(aug_a5predict1ws_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)
head(aug_a5predict1ws_vsci$.se.fit)
# Pineiro 2008 put observed on y and predict on x
# label specific points
# https://stackoverflow.com/questions/15624656/label-points-in-geom-point

select(aug_a5predict1ws_vsci, station_id_2, st_type, vsci, vsci_pred, .se.fit) %>% arrange(.se.fit) %>% print(n =48)

# aug_a5predict1ws_vsci <- aug_a5predict1ws_vsci %>%
#   mutate(
#     st_type = case_when(
#     station_id ==  station_id_2 ~ "status",
#     station_id != station_id_2 ~ "trend"
#   )
# )

# select(aug_a5predict1ws_vsci, station_id_2, st_type, vsci, vsci_pred, .se.fit) %>% arrange(st_type, .se.fit) %>% print(n =48)

ggplot(aug_a5predict1ws_vsci, aes(x = st_type, y = vsci)) + geom_boxplot() +labs(title = "Observed VSCI at n = 48 Preds Sites")

ggplot(aug_a5predict1ws_vsci, aes(x = st_type, y = vsci_pred)) + geom_boxplot() + labs(title = "Predicted VSCI at n = 48 Preds Sites")

ggplot(aug_a5predict1ws_vsci, aes(x = st_type, y = .se.fit)) + geom_boxplot() + labs(title = "Predicted VSCI at n = 48 Preds Sites")

ggplot(aug_a5predict1ws_vsci, aes(x = vsci_pred, y = vsci, size = .se.fit, shape = st_type)) + geom_point()

ggplot(aug_a5predict1ws_vsci, aes(x = vsci_pred, y = vsci, size = .se.fit, shape = st_type)) + geom_point()+ geom_smooth(aes(x = vsci_pred, y = vsci), method = "lm", se = F, color = "blue")

# using linewidth instead of size b/c ggplot geom_smooth waring
ggplot(aug_a5predict1ws_vsci, aes(x = vsci_pred, y = vsci, linewidth = .se.fit, shape = st_type)) + geom_point()+ geom_smooth(aes(x = vsci_pred, y = vsci), method = "lm", se = F, color = "blue")

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


```


# Preds for ST_rand model reml
Using reml estimate as testing predictions.
```{r st_rand1_hypothesis}

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

summary(ssn_wswq_rand1_st_reml)
varcomp(ssn_wswq_rand1_st_reml) # now space varcomp small
loocv(ssn_wswq_rand1_st_reml)

# augment returns an sf object
# note level set to 0.90
aug_st5predict1ws_vsci <- augment(ssn_wswq_rand1_st_reml, newdata = "Preds_2019_2022", se_fit = TRUE, interval = c("prediction"), level = 0.90)
# class(aug_a5predict1ws_vsci)
# names(aug_a5predict1ws_vsci)
aug_st5predict1ws_vsci$vsci_pred <- round(aug_st5predict1ws_vsci$.fitted,1)
aug_st5predict1ws_vsci$vsci_lwr <- round(aug_st5predict1ws_vsci$.lower,1)
aug_st5predict1ws_vsci$vsci_upr <- round(aug_st5predict1ws_vsci$.upper,1)
head(aug_st5predict1ws_vsci$.se.fit)
# Pineiro 2008 put observed on y and predict on x
# label specific points
# https://stackoverflow.com/questions/15624656/label-points-in-geom-point

select(aug_st5predict1ws_vsci, station_id_2, st_type, vsci, vsci_pred, .se.fit) %>% arrange(.se.fit) %>% print(n =48)

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

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

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

```


# TLH Fitted interval preds
On 01/10/2025 merged Travis' branch to make this as a ggplotly plot
Based on 12/24/2024 email this graph of preds is using standard error of model so that gives straight prediction interval lines
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

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

# combine to a dataframe note variance/se is same for all obs
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<-round(aug_a5predict1ws_vsci$vsci,1)
tx<-round(spat.Pred.median,1)

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

# put all components into data frame for ggplot
plot_dat = data.frame("vsci_pred" = tx, "vsci" = ty, "pred_5th" = op05, "pred_95th" = op95,"vsci_bound_order"=tx[ox],"station_id"=aug_a5predict1ws_vsci$station_id_2)


p = ggplot(data = plot_dat)+
  geom_point(aes(x = vsci_pred, y = vsci,label = station_id), color = "black")+
  geom_smooth(aes(x = vsci_pred, y = vsci), method = "loess", se = F, color = "blue")+
  geom_line(aes(x = vsci_bound_order, y = pred_5th), linetype = "dashed", color = "black")+
  geom_line(aes(x = vsci_bound_order, y = pred_95th), linetype = "dashed", color = "black")+
  geom_abline(intercept = 0, slope = 1, color = "red")+
  labs(x="vsci predicted",
       y ="vsci observed",
       title = "SSN Withheld Obs vs Pred a5 reml n=48 interval 0.90, tlh method")
print(p)
ggplotly(p)

# identify does not work in RStudio
# identify(tx, ty, labels = station_id_2, plot=TRUE)

# number of points outside 90% prediction limits
os90pl.sp <-(ty < spat.Pred.pct[,"Q05"]) |  (ty > spat.Pred.pct[,"Q95"]) 
summary(os90pl.sp)

# which observations are outside [5%, 95%] prediction limits?
# these are the Site IDs of the unusual values
aug_a5predict1ws_vsci[which(os90pl.sp),"station_id_2"]

mapview(aug_a5predict1ws_vsci[which(os90pl.sp),"station_id_2"])

rm(fitted_vals, resid_SD, spat, spat.Pred.pct, spat.Pred.median, ty, tx, ox, op05, op95, plot_dat, p, os90pl.sp)

# plot using base plot with all components called separately into the 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 Withheld Obs vs. Pred a5 Model SE 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)
```



# TLH Fitted interval preds (quasi-LOOCV method)
On 01/10/2025 merged Travis'branch to make this as a ggplotly plot

On 12/27/2024, is ribbon prediction interval in chunk 20, Preds for Distance Hypothesis, equivalent to chunk 22, fitted_interval_plot_loocv-preds, below? Doing head(aug_a5predict1ws_vsci$vsci_lwr) and same to vsci_upr show ribbon prediction intervals wider than Q05 and Q95. Consequently, the narrower Q05 and Q95 prediction intervals identify more points outside the intervals.

Based on 12/24/2024 email from Travis this uses standard error of each prediction point so get squiggly lines.
```{r fitted_interval_plot_loocv-preds}
# Alrighty, so this method mimics John's method
# which gives a loocv estimate of fitted value an se.
# However, as we're looking at new predicted data it doesn't actually refit the model n times so we don't have proper loocv fitted values.
# This results in the squiggly lines because he's giving that 
# se at each individual value.
fitted_vals = aug_a5predict1ws_vsci$.fitted
resid_Se = aug_a5predict1ws_vsci$.se.fit

# combine to a dataframe and se.fit/variance differs for each observation
spat<-cbind(fitted_vals,resid_Se^2)
head(spat)

# spatial prediction percentiles
# qnorm arguments include
# p a vector of probabilities
# mean a vector of means, which are our fitted values
# sd a vector of standard deviations, the square root of the squared .se.fit
(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)


# observed and median of our prediction intervals
ty<-round(aug_a5predict1ws_vsci$vsci,1)
tx<-round(spat.Pred.median,1) ## these just are the fitted values

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

# put all components into data frame for ggplot
plot_dat = data.frame("vsci_pred" = tx, "vsci" = ty, "pred_5th" = op05, "pred_95th" = op95,"vsci_bound_order"=tx[ox],"station_id"=aug_a5predict1ws_vsci$station_id_2)


p = ggplot(data = plot_dat)+
  geom_point(aes(x = vsci_pred, y = vsci,label = station_id), color = "black")+
  geom_smooth(aes(x = vsci_pred, y = vsci), method = "loess", se = F, color = "blue")+
  geom_line(aes(x = vsci_bound_order, y = pred_5th), linetype = "dashed", color = "black")+
  geom_line(aes(x = vsci_bound_order, y = pred_95th), linetype = "dashed", color = "black")+
  geom_abline(intercept = 0, slope = 1, color = "red")+
  labs(x="vsci predicted",
       y ="vsci observed",
       title = "SSN Withheld Obs vs Pred a5 reml n=48 interval 0.90, jhc method")

print(p)
ggplotly(p)

# Carson Code Points Outside LOOCV 90% Prediction Limit for n=48 Withheld Obs
# This code comes from local_SSN_FINAL_20160126_modelling.html that was used in the Scown et al. 2017 paper.

# number of points outside 90% prediction limits
os90pl.sp <-(ty < spat.Pred.pct[,"Q05"]) |  (ty > spat.Pred.pct[,"Q95"]) 
summary(os90pl.sp)

# which observations are outside [5%, 95%] prediction limits?
# these are the Site IDs of the unusual values
aug_a5predict1ws_vsci[which(os90pl.sp),"station_id_2"]

mapview(aug_a5predict1ws_vsci[which(os90pl.sp),"station_id_2"])

rm(fitted_vals, resid_Se, spat, spat.Pred.pct, spat.Pred.median, ty, tx, ox, op05, op95, plot_dat, p, os90pl.sp)

# plot using base plot with all components called in separatel to 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 Withheld Obs vs. Pred a5 reml se.fit 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")

```

# Ribbon Prediction Interval
This is producing a 90% prediction interval because that was the level set when augmented predictions were made. Is the spatial prediction percentiles using qnorm from TLH and JC code more accurate?
```{r ribbon_pred_interval}
# ribbon prediction interval based on + or - 1.96*.se.fit
# ggplot(aug_a5predict1ws_vsci, aes(x = vsci_pred, y = vsci, label = station_id_2)) + geom_point() + geom_abline() + geom_ribbon(aes(ymin = vsci_pred - (1.96*.se.fit), ymax = vsci_pred + (1.96*.se.fit)), fill = "grey", alpha = .8) + geom_line(aes( y = vsci_pred)) + labs(x = "VSCI Predicted 2019-2022", y = "VSCI Observed 2019-2022", title = "SSN Withheld Obs vs Pred a5 reml n=48")

# ribbon using lower and upper SSN2 prediction intervals
# seems to give same results as 1.96*.se.fit above
ribbon_pi1 <- ggplot(aug_a5predict1ws_vsci, aes(x = vsci_pred, y = vsci, label = station_id_2)) + geom_point() + geom_abline() + geom_ribbon(aes(ymin = vsci_lwr, ymax = vsci_upr), fill = "grey", alpha = .8) + geom_line(aes( y = vsci_pred)) + labs(x = "VSCI Predicted 2019-2022", y = "VSCI Observed 2019-2022", title = "SSN Withheld Obs vs Pred a5 reml n=48 interval 0.90, mgm method")

print(ribbon_pi1)
ggplotly(ribbon_pi1)

# number of points outside 90% prediction limits (opl)
# opl_stations1 <- filter(aug_a5predict1ws_vsci, vsci < (vsci_pred - (1.96*.se.fit)) | vsci > (vsci_pred + (1.96*.se.fit)))

opl_stations2 <- filter(aug_a5predict1ws_vsci, vsci < vsci_lwr | vsci > vsci_upr)

class(opl_stations2)
mapview(opl_stations2)
opl_stations2
# given comment "Target and Sampled (no fish or PHAB in 2019?) Dry,no staff", should this site 2-DCK003.94_2019 been sampled?
```