all_eval.Rmd

---
output: 
  html_document:
    code_folding: hide
self_contained: yes
runtime: shiny
---

# Bayesian network tool {.tabset}

```{r globals, echo = T, message = F, warning = F}
knitr::opts_chunk$set(message = F, warning = F)

library(tidyverse)
library(ggmap)
library(lubridate)
library(geosphere)
library(stringi)
library(tibble)
library(bnlearn)
library(scales)
library(shiny)
library(sf)
library(sp)
library(gridExtra)
library(MuMIn)

data(restdat)
data(reststat)
data(wqdat)
data(wqstat)
data(map)
data(grdsres)

# source R files
source('R/get_chg.R')
source('R/get_clo.R')
source('R/get_dat.R')
source('R/get_lik.R')

# base map
# ext <- make_bbox(wqstat$lon, wqstat$lat, f = 0.2)
# get_stamenmap(ext, zoom = 11, maptype = "toner-lite", where = getwd())
pbase <- ggmap(map) +
  theme_bw() +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  )
```

## Interactive {.tabset}

```{r sliders}
## sliders

column(width = 12, 
  
  column(width = 6,
    
    # action button 
    h5('Hit once on startup, once for changing inputs:'),
    actionButton("apply", "Apply changes"),
    
    # model plot
    renderPlot({plot(net())}, width = 450, height = 400)

  ), 
  
  column(width = 6, 

    # which model to select
    selectInput("resgrp", label = 'Select model type:', 
                choices = list(
                  simple = 'top', 
                  complex = 'type'
                ), selected = 'simple'
              ),
    
    # wq year slice window
    numericInput("yrdf", "Choose year diff:", 5), 
    
    # wq site match upper limit
    numericInput("mtch", "Choose number of restoration matches:", 2),
    
    # year slider
    sliderInput("yrs", label = 'Select year ranges:',  
          min = 1974, max = 2016, 
          value = c(1974, 2016),
          sep = '', ticks = FALSE
        ), 
    
    # dynamic UI for chl thresh
    renderUI({
      
      req(uirngs())
      
      # ui rng inputs
      uirngs <- uirngs()
      
      # filter for chl
      chlrng <- uirngs %>% 
        filter(var == 'chval')
      
      # checker for default median, otherwise selected as input
      if (is.null(input$chlspl))
        valin <- chlrng$medv
      else
        valin <- input$chlspl

      # chl selector
      sliderInput("chlspl", label = 'Select chl lo/hi split:',
        min = chlrng$minv, max = chlrng$maxv,
        value = valin,
        sep = '', ticks = FALSE
        )

      }), 
    
    # dynamic UI for nit thresh
    renderUI({
      
      req('nival' %in% uirngs()$var)
      
      # ui rng inputs
      uirngs <- uirngs()
      
      # filter for nit
      nitrng <- uirngs %>% 
        filter(var == 'nival')
  
      # checker for default median, otherwise selected as input
      if (is.null(input$nitspl))
        valin <- nitrng$medv
      else
        valin <- input$nitspl
      
      # nit selector
      sliderInput("nitspl", label = 'Select nitrogen lo/hi split:',
        min = nitrng$minv, max = nitrng$maxv,
        value = valin,
        sep = '', ticks = FALSE
        )
    
      }), 
    
    # dynamic UI for sal thresh
    renderUI({
      
      req(uirngs())
      
      # ui rng inputs
      uirngs <- uirngs()
      
      # filter for sal
      salrng <- uirngs %>% 
        filter(var == 'saval')
      
      # checker for default median, otherwise selected as input
      if (is.null(input$salspl))
        valin <- salrng$medv
      else
        valin <- input$salspl
      
      # sal selector
      sliderInput("salspl", label = 'Select sal lo/hi split:',
        min = salrng$minv, max = salrng$maxv,
        value = valin,
        sep = '', ticks = FALSE
        )
    
      })
    
  )
  
)
```

```{r reactives}
## reactives

# bn class reactive
net <- eventReactive(input$apply, { 
  
  # get net structure from model type
  resgrp <- input$resgrp
  if(resgrp == 'top')
    resgrp <- "[hab][wtr][salev|hab:wtr][chlev|hab:wtr:salev]"
  else 
    resgrp <- "[hab_enh][hab_est][hab_pro][non_src][pnt_src][salev|hab_enh:hab_est:hab_pro:non_src:pnt_src][nilev|salev:hab_enh:hab_est:hab_pro:non_src:pnt_src][chlev|nilev]"
                                     
  # create bn network object
  out <- model2network(resgrp)

  return(out)
  
})

# nominal category ranges from inputs, for dynamic ui selector of thresholds
uirngs <- eventReactive(input$apply, { 
  
  req(cdat())
  
  # inputs
  cdat <- cdat()
  
  # get ranges from inputs   
  rngs <- cdat %>% 
    select(matches('^saval$|^chval$|^nival$')) %>% 
    gather('var', 'val') %>% 
    group_by(var) %>% 
    summarize(
      minv = round(min(val, na.rm = T), 2),
      medv = round(quantile(val, 0.5, na.rm = T), 2),
      maxv = round(max(val, na.rm = T), 2)
    )
  
  return(rngs)
  
})

# reactives for restoration data subsets, data
restdat_sub <- eventReactive(input$apply, { 
  
  out <- restdat %>% 
    filter(date >= input$yrs[1] & date <= input$yrs[2])
  
  return(out)
  
})

# reactives for restoration data subsets, station lat/lon
reststat_sub <- eventReactive(input$apply, { 
  
  out <- reststat %>% 
    filter(id %in% restdat_sub()$id)
  
  return(out)
  
})

# wq stats matched to rest stats
allcdat <- eventReactive(input$apply, { 
  
  # reactive inputs
  mtch <- input$mtch
  yrdf <- input$yrdf
  resgrp <- input$resgrp
  chlspl <- input$chlspl
  nitspl <- input$nitspl
  salspl <- input$salspl
  restdat_sub <- restdat_sub()
  reststat_sub <- reststat_sub()

  # get cdat and wqmtch
  out <- get_dat(resgrp, restdat_sub, reststat_sub, wqstat, wqdat, mtch, yrdf, chlspl, nitspl, salspl)

  return(out)
  
})

# get conditional data for bn input
cdat <- eventReactive(input$apply, { 

  allcdat()$cdat
  
})

# get wqmtch for map
wqmtch <- eventReactive(input$apply, { 

  allcdat()$wqmtch

})

# fit simple bn model
cdat_mod <- eventReactive(input$apply, {

  # inputs
  cdat <- cdat()

  # format cdat for model input
  cdat_frm <- cdat %>%
    select_if(is.factor) %>%
    dplyr::select(-stat) %>% 
    na.omit %>%
    data.frame

  # fitted mod
  out <- bn.fit(net(), data = cdat_frm)

  return(out)

})

# get likelihood estimates from mod, bef/aft only (no salinity)
ests <- eventReactive(input$apply, {

  # inputs
  cdat <- cdat() # conditional data for mod input
  cdat_mod <- cdat_mod() # bn model

  # get likelihood
  out <- get_lik(cdat, cdat_mod)

  return(out)
  
})
```


### Spatial matches

```{r mapplo}
mapplo <- eventReactive(input$apply, {

  # reactives
  resgrp <- input$resgrp
  wqmtch <- wqmtch()
  restdat_sub <- restdat_sub()
  reststat_sub <- reststat_sub()
  
  # combine lat/lon for the plot
  toplo <- wqmtch %>%
    left_join(wqstat, by = 'stat') %>%
    left_join(reststat_sub, by = 'id') %>%
    rename(
      `Restoration type` = resgrp,
      `Distance (dd)` = dist
    )

  # restoration project grouping column
  restall <- left_join(restdat_sub, reststat_sub, by = 'id')
  names(restall)[names(restall) %in% resgrp] <- 'Restoration type'

  # outplot
  p <- pbase +
    geom_point(data = restall, aes(x = lon, y = lat, fill = `Restoration type`), size = 4, pch = 21) +
    geom_point(data = wqstat, aes(x = lon, y = lat), size = 2) +
    geom_segment(data = toplo, aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = -`Distance (dd)`, linetype = `Restoration type`), size = 1) +
    scale_alpha(guide = F) +
    theme_bw(base_family = 'serif', base_size = 16)

  return(p)

})

renderPlot({mapplo()}, width = 550, height = 650)

```

### Chl response

```{r chlplo}
chlplo <- eventReactive(input$apply, {

  #inputs
  ests <- ests()

  # format for plot
  toplo <- ests %>%
    mutate(event = gsub('^wtr\\_|^hab\\_|^hab\\_enh\\_|^hab\\_est\\_|^hab\\_pro\\_|^non\\_src\\_|^pnt\\_src\\_', '', event)) %>%
    mutate(
      chlev = factor(chlev, levels = c('lo', 'hi')),
      salev = factor(salev, levels = c('lo', 'hi'), labels = c('Salinity lo', 'Salinity hi'))
      ) %>%
    spread(event, est) %>%
    mutate(
      chg = 100 * (aft - bef),
      chg = round(chg, 1)
      )
  
  # barplots of % change bef/aft
  p <- ggplot(toplo, aes(x = chlev, y = chg, group = project)) +
    geom_bar(colour = 'black', stat = 'identity', aes(fill = chg)) +
    theme_bw(base_family = 'serif', base_size = 16) +
    theme(
      strip.background = element_blank(),
      panel.grid = element_blank(),
      legend.position = 'none'
    ) +
    geom_hline(yintercept = 0) +
    xlab('Chlorophyll') +
    ylab('% Change after - before') +
    scale_fill_gradient2('% change', low = 'lightgreen', mid = "white", high = 'lightblue', midpoint = 0) +
    facet_grid(salev ~ project)

  return(p)
  
})

renderPlot({chlplo()}, width = 600, height = 600)

```

## Static

```{r sliders2}

column(width = 12, 
       
  # select model type
  column(width = 6, 
    
    # select model type
    selectInput("resgrp2", label = 'Select model type:', choices = c('simple', 'complex'), selected = 'simple'),
    
    # select first mtch    
    renderUI({
      
      # inputs
      resgrp2 <- input$resgrp2
      yrgrp <- input$yrgrp
      
      # get yrs from model type selection
      mtch <- grdsres %>% 
        filter(resgrp == resgrp2) %>% 
        filter(yrs %in% yrgrp) %>% 
        pull(mtch) %>% 
        unique
    
      selectInput("mtch1", label = 'First number of restoration matches:', choices = mtch, selected = max(mtch))
      
    })
    
  ),
  
  column(width = 6, 
         
    # select year group     
    renderUI({
      
      # inputs
      resgrp2 <- input$resgrp2
      
      # get yrs from model type selection
      yrs <- grdsres %>% 
        filter(resgrp == resgrp2) %>% 
        pull(yrs) %>% 
        unique
      
      selectInput("yrgrp", label = 'Select year group:', choices = yrs)
      
    }),
    
    # select second mtch
    renderUI({

      # inputs
      resgrp2 <- input$resgrp2
      yrgrp <- input$yrgrp

      # get yrs from model type selection
      mtch <- grdsres %>%
        filter(resgrp == resgrp2) %>%
        filter(yrs %in% yrgrp) %>%
        pull(mtch) %>%
        unique

      selectInput("mtch2", label = 'Second number of restoration matches:', choices = mtch, selected = max(mtch))

    })
    
  )
  
)
  
```

```{r reactives2}
# plot data matches
plostamtc <- reactive({
  
  req(input$mtch1)
  
  # inputs
  resgrp2 <- input$resgrp2
  yrgrp <- input$yrgrp
  mtch1 <- input$mtch1
  mtch2 <- input$mtch2

  # filter data to plot
  plosta <- grdsres %>% 
    filter(resgrp == resgrp2) %>% 
    filter(yrs == yrgrp) %>% 
    filter(mtch %in% as.numeric(c(mtch1, mtch2))) %>% 
    mutate(
      salev = factor(salev, levels = c('lo', 'hi'), labels = c('Sal lo', 'Sal hi')),
      chlev = factor(chlev, levels = c('lo', 'hi'), labels = c('Chl lo', 'Chl hi'))
    )
  
  return(plosta)
  
})

# plot data yrdf
plostatim <- reactive({
  
  req(input$yrdf1)
  
  # inputs
  resgrp2 <- input$resgrp2
  yrgrp <- input$yrgrp
  yrdf1 <- input$yrdf1
  yrdf2 <- input$yrdf2

  # filter data to plot
  plosta <- grdsres %>% 
    filter(resgrp == resgrp2) %>% 
    filter(yrs == yrgrp) %>% 
    filter(yrdf %in% as.numeric(c(yrdf1, yrdf2))) %>% 
    mutate(
      salev = factor(salev, levels = c('lo', 'hi'), labels = c('Sal lo', 'Sal hi')),
      chlev = factor(chlev, levels = c('lo', 'hi'), labels = c('Chl lo', 'Chl hi'))
    )
  
  return(plosta)
  
})
```

Category splits for chlorophyll fixed at 11 ug/L, nitrogen at 0.5 mg/L, salinity at 27 psu.

```{r plostamtc}
renderPlot({
  
  # input
  plostamtc <- plostamtc()

  p <- ggplot(plostamtc, aes(x = yrdf, y = chg, colour = chlev)) +
    # geom_line(aes(linetype = factor(mtch))) +
    geom_point(size = 4, aes(shape = factor(mtch)), alpha = 0.7) +
    stat_smooth(method = 'lm', se = F, aes(linetype = factor(mtch))) +
    facet_grid(salev ~ project) +
    theme_bw(base_size = 16) +
    theme(
      strip.background = element_blank(), 
      legend.position = 'top'
    ) + 
    scale_colour_discrete("Chlorophyll zone") +
    scale_linetype_discrete('Number of restoration site matches') + 
    scale_shape_manual('Number of restoration site matches', values = c(16, 17)) + 
    scale_x_continuous("Year difference") +
    scale_y_continuous("% Change after - before") +
    geom_hline(yintercept = 0) +
    ggtitle('Results by time difference')
  
  return(p)
    
}, height = 500, width = 800)

```

```{r sliders3}

column(width = 12, 
       
  # select model type
  column(width = 6, 
    
    # select first yrdf
    renderUI({

      # inputs
      resgrp2 <- input$resgrp2
      yrgrp <- input$yrgrp

      # get yrdf from model type selection
      yrdf <- grdsres %>%
        filter(resgrp == resgrp2) %>%
        filter(yrs %in% yrgrp) %>%
        pull(yrdf) %>%
        unique

      selectInput("yrdf1", label = 'First number of year differences:', choices = yrdf, selected = max(yrdf))

    })
    
  ),
  
  column(width = 6, 
    
    # select second yrdf
    renderUI({

      # inputs
      resgrp2 <- input$resgrp2
      yrgrp <- input$yrgrp

      # get yrdf from model type selection
      yrdf <- grdsres %>%
        filter(resgrp == resgrp2) %>%
        filter(yrs %in% yrgrp) %>%
        pull(yrdf) %>%
        unique

      selectInput("yrdf2", label = 'Second number of year differences:', choices = yrdf, selected = max(yrdf))

    })
    
    
  )
  
)
```

```{r plostatim}
renderPlot({
  
  # input
  plostatim <- plostatim()

  p <- ggplot(plostatim, aes(x = mtch, y = chg, colour = chlev)) +
    # geom_line(aes(linetype = factor(yrdf))) +
    geom_point(size = 4, aes(shape = factor(yrdf)), alpha = 0.7) +
    stat_smooth(method = 'lm', se = F, aes(linetype = factor(yrdf))) +
    facet_grid(salev ~ project) +
    theme_bw(base_size = 16) +
    theme(
      strip.background = element_blank(), 
      legend.position = 'top'
    ) + 
    scale_colour_discrete("Chlorophyll zone") +
    scale_linetype_discrete('Time differences') + 
    scale_shape_manual('Time differences', values = c(16, 17)) + 
    scale_x_continuous("Number of site matches") +
    scale_y_continuous("% Change after - before") +
    geom_hline(yintercept = 0) +
    ggtitle('Results by site match differences')
  
  return(p)
    
}, height = 500, width = 800)

```