03_StabilityAnalysis.Rmd

---
title: "3_StabilityAnalysis"
author: "Simon Topp"
date: "9/12/2020"
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
library(tidyverse)
library(lubridate)
library(purrr)
library(furrr)
library(data.table)
library(feather)
library(sf)
library(dtwclust)
library(networkD3)
library(tictoc)
knitr::opts_chunk$set(echo = TRUE)
```

## Make the sankey Network

```{r}
clusters.sf <- st_read(paste0('data/out/',iterations,'_clustersSF.shp'))

links.base <- clusters.sf %>% st_set_geometry(NULL) %>%
  select(Hylak_id, period, cluster) %>%
  arrange(period) %>%
  pivot_wider(names_from = period, values_from = cluster) 

p1Tab <- as.data.frame(table('source' = links.base[[2]],'target' = links.base[[3]])) %>%
  mutate(source = paste0('p1_',source), target = paste0('p2_',target))

p2Tab <- as.data.frame(table('source' = links.base[[3]],'target' = links.base[[4]])) %>%
  mutate(source = paste0('p2_',source), target = paste0('p3_',target))

p3Tab <- as.data.frame(table('source' = links.base[[4]],'target' = links.base[[5]])) %>%
  mutate(source = paste0('p3_',source), target = paste0('p4_',target))

p4Tab <- as.data.frame(table('source' = links.base[[5]],'target' = links.base[[6]])) %>%
  mutate(source = paste0('p4_',source), target = paste0('p5_',target))

p5Tab <- as.data.frame(table('source' = links.base[[6]],'target' = links.base[[7]])) %>%
  mutate(source = paste0('p5_',source), target = paste0('p6_',target))

links <- bind_rows(p1Tab,p2Tab,p3Tab,p4Tab,p5Tab)

nodes <- data.frame(
  name=c(as.character(links$source), as.character(links$target)) %>% 
    unique()
  ) %>%
  mutate(nodeGroup = ifelse(grepl('Spring',name), 1, 
                            ifelse(grepl('Summer',name), 2, 
                                   ifelse(grepl('Bimodal',name),3,
                                          ifelse(grepl('(Blue)',name),4,
                                                 ifelse(grepl('(Green)',name),5,NA))))),
         nodeGroup = factor(nodeGroup))

groupLabels <- tibble(nodeGroup = factor(c(1:5)), label = c('Spring Greening', 'Summer Greening', 'Bimodal', 'Aseasonal (Blue)', 'Aseasonal (Green)'))

nodes <- nodes %>% left_join(groupLabels)
 
# With networkD3, connection must be provided using id, not using real name like in the links dataframe. So we need to reformat it.
links$IDsource <- match(links$source, nodes$name)-1 
links$IDtarget <- match(links$target, nodes$name)-1
 
# Make the Network
# Define the color scale
cs <- 'd3.scaleOrdinal() .range(["#440154FF", "#2A788EFF", "#7AD151FF"])'
p <- sankeyNetwork(Links = links, Nodes = nodes,
                     Source = "IDsource", Target = "IDtarget",
                     Value = "Freq", NodeID = "name", NodeGroup = 'nodeGroup',
                     sinksRight=F, nodeWidth = 10, nodePadding = 30,
                   colourScale = cs, fontSize = 12, fontFamily = 'Arial')

p


## Try it with ggalluvial
library(ggalluvial)

links.alluvial <- links.base %>%
  select(Hylak_id:`(2014,2020]`) %>%
  group_by(`(1984,1990]`, `(1990,1996]`, `(1996,2002]`, `(2002,2008]`, `(2008,2014]`, `(2014,2020]`) %>%
  summarise(count = n()) %>%
  ungroup() %>%
  mutate(alluvium = row_number()) %>%
  pivot_longer(-c(alluvium, count), names_to = 'Period', values_to = 'Cluster') %>%
  mutate(Cluster = factor(Cluster, levels = c('Spring Greening', 'Summer Greening', 'Bimodal', 'Aseasonal (Blue)', 'Aseasonal (Green)')))


links.alluvial %>%
  ggplot(aes(x = factor(Period), y = count, stratum = Cluster, label = Cluster, 
             alluvium = alluvium, fill = Cluster)) +
  geom_flow() +
  geom_stratum(alpha = .8) +
  #scale_fill_manual(values = c("#440154FF", "#25848EFF","#38B977FF", "#3E4B89FF", "#BBDF27FF")) +
  scale_fill_manual(values = c("#5F4690", "#EDAD08","#CC503E", "#38A6A5", "#73AF48")) +
  labs(x = 'Period', y = 'Lake Proportion') + 
  theme_classic() +
  theme(legend.position = 'top',
        axis.ticks = element_blank(),
        axis.line = element_blank(),
        axis.text.y = element_blank(),
        legend.title = element_blank())

ggsave('figures/transitionPlotGG_V2.png', width = 6.5, height = 4.5, units = 'in')

```

## Look at stability metrics from cluster transitions

```{r}
## Pull out summary stats
links.summary <- links %>%
  group_by(source) %>%
  mutate(total = sum(Freq),
         PChange = Freq/total) %>%
  ungroup() %>%
  left_join(nodes, by = c('source' = 'name')) %>%
  rename(nodeSource = nodeGroup) %>% left_join(nodes, by = c('target' = 'name')) %>%
  rename(nodeTarget = nodeGroup) %>%
  group_by(label.x, label.y) %>%
  summarise(mean = mean(PChange),
            sd = sd(PChange))
  
## Get inter vs intra numbers for statistical analysis
interVintra <- links %>%
  group_by(source) %>%
  mutate(total = sum(Freq),
         PChange = Freq/total) %>%
  ungroup() %>%
  left_join(nodes, by = c('source' = 'name')) %>%
  rename(nodeSource = nodeGroup) %>%left_join(nodes, by = c('target' = 'name')) %>%
  rename(nodeTarget = nodeGroup) %>%
  mutate(group = ifelse(nodeSource == nodeTarget, 'intra', 'inter'),
         group = paste0(group, nodeSource))


#Take a look at them
interVintra %>% group_by(group, label.x) %>%
  summarise(mean = mean(PChange),
            sd = sd(PChange))

transition.summaries <- interVintra %>% group_by(label.x, label.y) %>%
  summarise(mean = mean(PChange),
            sd = sd(PChange))


dt <- FSA::dunnTest(PChange~group, data = interVintra %>% filter(grepl('intra', group)), method = 'bonferroni')
print(dt,  dunn.test.results = T)


kruskal.test(PChange ~ group, interVintra)
dt <- FSA::dunnTest(PChange~group, data = interVintra, method = 'bonferroni')
print(dt,  dunn.test.results = T)

# Look at total number of changes
links.base <- links.base %>%
  mutate(change1 = ifelse(`(1984,1990]` == `(1990,1996]`, 0,1),
         change2 = ifelse(`(1990,1996]` == `(1996,2002]`, 0,1),
         change3 = ifelse(`(1996,2002]` == `(2002,2008]`, 0,1),
         change4 = ifelse(`(2002,2008]` == `(2008,2014]`, 0,1),
         change5 = ifelse(`(2008,2014]` == `(2014,2020]`, 0,1),
         totalChange = factor(change1 + change2 + change3 + change4 + change5))

summary(links.base$totalChange)

# Look at number of unique states
stability <- links.base %>% 
  pivot_longer(`(1984,1990]`:`(2014,2020]`, names_to = 'Period', 
               values_to = 'Cluster') %>% 
  group_by(Hylak_id) %>% 
  summarize(States = factor(length(unique(Cluster)), 
                            levels = c(1:5), labels = c(1:5))) %>%
  left_join(links.base %>% select(Hylak_id, totalChange))
 
summary(stability$States)

ggplot(stability, aes(x = factor(totalChange), fill = States)) + 
  geom_bar() +
  scale_fill_viridis_d(end = .8) +
  labs(x = 'Number of Cluster Transitions', y = 'Lake Count', 
       fill = 'Unique\nClusters', title = 'Lake Transitions and States') +
  theme_bw()

ggsave('figures/transitionCounts.png', width = 3, height = 3, units = 'in')

## Make a final stability dataset linked to hydrolakes and GLCP

stability <- stability %>% 
  inner_join(hl %>% mutate(lwRatio = Lake_area/Wshd_area)) %>% 
  st_as_sf() 

## Bring in the GLCP
## Available at https://environmentaldatainitiative.org/edis-featured-data-contributions/glcp-dataset/
# glcp <- read_csv('data/in/glcp.csv')
# glcp.filt <- as.data.table(glcp)[Hylak_id %in% stability$Hylak_id]
# write_feather(glcp.filt, 'data/in/glcpFiltered.feather')

glcp <- read_feather('data/in/glcpFiltered.feather') %>%
  mutate(spRatio = seasonal_km2/permanent_km2)

# Since timespans don't match up, just calculate mean as cv values for each variable
glcp.avg <- glcp %>% group_by(Hylak_id) %>% summarise_at(vars(mean_monthly_precip_mm:spRatio), c(mean = mean, cv = raster::cv), na.rm = T)

stability <- stability %>% left_join(glcp.avg) 

write_feather(stability %>% st_set_geometry(NULL), paste0('data/out/',iterations,'_stability.feather'))
```

## Look at the relationships between stability and lake and landscape variables

```{r}
stability <- read_feather(paste0('data/out/',iterations,'_stability.feather'))

names(stability)

stability.long <- stability %>% group_by(totalChange) %>%
  summarise_at(vars(Lake_area:Wshd_area, lwRatio, mean_monthly_precip_mm_mean:spRatio_cv), median, na.rm = T)
  
stability.long <- stability.long %>% pivot_longer(Lake_area:spRatio_cv, names_to = 'Variable')

lms <- stability.long %>% 
  mutate(totalChange = as.numeric(as.character(totalChange))) %>%
  filter(!Variable %in% c('Lake_area', 'Vol_src', 'Vol_res', 'total_precip_mm_mean', 'total_precip_mm_cv')) %>% #remove redundant variables and the total precip from GLCP cause it's redundant with monthly
  group_by(Variable) %>% 
  nest() %>%
  mutate(lm = purrr::map(data, ~lm(totalChange ~ value, .)),
         summary = purrr::map(lm, broom::tidy)) %>%
  select(-c(data,lm)) %>%
  unnest(summary)

view(lms %>% filter(term == 'value'))

unique(lms$Variable)

renamer <- tibble(Variable = c(
"Shore_len"                  , "Shore_dev"                  , "Vol_total"              ,    
"Depth_avg"                  , "Dis_avg"                    , "Res_time"               ,    
"Elevation"                  , "Slope_100"                  , "Wshd_area"              ,    
"lwRatio"                    , "mean_monthly_precip_mm_mean", "mean_annual_temp_k_mean", 
"pop_sum_mean"               , "seasonal_km2_mean"          , "permanent_km2_mean"     , 
"total_km2_mean"             , "spRatio_mean"               , "mean_monthly_precip_mm_cv",
"mean_annual_temp_k_cv"      , "pop_sum_cv"                 , "seasonal_km2_cv"          , 
"permanent_km2_cv"           , "total_km2_cv"               , "spRatio_cv"),

namesNew = c(
"Shore Length (km)"        , "Shore Development"    ,"Volume (cu. km)"      , 
"Depth (m)"            , "Discharge (cu. m/s)"      , "Residence Time (days)"    ,
"Elevation (m)"        , "Surrounding Slope"        , "Watershed Area (sq. km)"  ,
"Lake/Watershed Ratio" , 'Mean Monthly Precip (mm)' , "Mean Temperature (k)" , 
"Population"           , "Seasonal Water (sq. km)"  , "Permanent Waters (sq. km)", 
"Total Area (sq. km)"  , "Seasonal/Permanent Ratio" , "CV of Monthly Precip" ,
"CV of Mean Annual Temp", "CV of Population"        , "CV of Seasonal Waters"    , 
"CV of Permanent Water", "CV of Total Area"         , "CV of Seasonal/Permanent Ratio"))

# After talking to Michael from GLCP, Total_precip isn't what we thought it was, so remove it.
#"total_precip_mm_cv" ,"total_precip_mm_mean"   ,
#"CV of Annual Precip" ,"Mean Annual Precip (mm)"  , 


# Look at cross-correlations
png(filename = 'figures/CorPlot.png', type = 'windows', width = 11, height = 10, units = 'in', res = 200)#, width = 6.5, units = 'in', type = 'windows', res = 72) 

corPlot <- stability %>%
  select(Hylak_id, renamer$Variable) %>% 
  pivot_longer(-Hylak_id, names_to = 'Variable') %>%
  left_join(renamer) %>%
  select(namesNew, value, Hylak_id) %>%
  pivot_wider(everything(), names_from = namesNew) %>%
  select(-Hylak_id) %>%
  cor(., use = "complete.obs") %>%
  corrplot::corrplot(., method = 'square', type="upper", order="hclust", 
         tl.col="black", tl.srt=45, #Text label color and rotation
         # Combine with significance
         #p.mat = p.mat, sig.level = 0.01, insig = "blank", 
         # hide correlation coefficient on the principal diagonal
         diag=FALSE )
dev.off()

lms <- lms %>% filter(term == 'value') %>% left_join(renamer) %>%
  arrange(p.value)

stab.sum <- stability %>%
  select(totalChange,renamer$Variable) %>%
  pivot_longer(-totalChange, names_to = 'Variable') %>%
  group_by(totalChange, Variable) %>%
  summarise(quant25 = quantile(value, .25, na.rm = T),
            quant75 = quantile(value, .75, na.rm = T),
            median = median(value, na.rm = T)) %>% ungroup()

stab.sum %>%
  filter(Variable %in% lms$Variable[lms$p.value < 0.05]) %>% #[1:6]) %>% #
  left_join(renamer) %>%
  filter(namesNew != 'Permanent Waters (sq. km)') %>% # This is super reduntant so lets remove it.
  ggplot(aes(x = totalChange, y = median)) +
  #geom_point() +
  #geom_pointrange(aes(ymin = quant25, ymax = quant75)) +
  geom_crossbar(aes(ymin = quant25, ymax = quant75), fill = 'grey95') +
  #scale_y_continuous(trans = 'log10') +
  theme_bw()+
  theme(panel.grid.minor = element_blank()) +
  labs(x = 'Stability Class (0 = More Stable, 5 = Less Stable)', y = 'Metric Medians' ) +
  facet_wrap(~namesNew, scales = 'free', labeller = labeller(namesNew = label_wrap_gen(15)))

ggsave('figures/SigCors_v2.png', width = 4, height = 4, units = 'in')

lms %>% ungroup() %>%
  select(Variable = namesNew, Coefficient = estimate, Std.Error = std.error, p.value) %>%
  arrange(p.value) %>%
  mutate_at(vars(Coefficient:p.value), ~round(.,3)) %>%
  knitr::kable(., 'html') %>%
  kableExtra::column_spec(1:4,width = '1in') %>%
  kableExtra::kable_styling("striped") %>%
  kableExtra::as_image(width = 4, file = 'figures/lmTable.png')
```


## Make Figure showing the spatial distribution and stability of each cluster
## Overall
```{r}
usa <- maps::map('usa', plot = F) %>% st_as_sf() %>% st_transform(102003) 

grid <- st_make_grid(usa, cellsize = c(75000,75000), square = F) %>% st_as_sf() %>% mutate(ID = row_number())
 
grid <- grid %>% st_join(clusters.sf %>% st_transform(st_crs(grid)), left = F)

Modes <- function(x) {
  ux <- unique(x)
  tab <- tabulate(match(x, ux))
  mode <- ux[tab == max(tab)]
  mode <- ifelse(length(mode) > 1, 'Mixed', as.character(mode))
  return(as.factor(mode))
}

modalClust <- grid %>% st_set_geometry(NULL) %>%
  group_by(ID) %>%
  summarise(modalClust = Modes(cluster))

gridClusters <- grid %>% inner_join(modalClust %>% 
                                      mutate(modalClust = factor(modalClust, levels = c('Spring Greening','Summer Greening','Bimodal','Aseasonal', 'Mixed'))))
           
ggplot() + 
  geom_sf(data = usa) +
  geom_sf(data = gridClusters, aes(fill = modalClust)) +
  geom_sf(data = Ecoregs, fill = 'transparent', color = 'red', size = 1) +
  scale_fill_viridis_d('Modal\nCluster') +
  ggthemes::theme_map(base_size = 11) +
  theme(legend.position = 'bottom')


ggsave('figures/DominantClusterv2.png', width = 4.5, height = 4, units = 'in')


## Stability
grid <- st_make_grid(usa, cellsize = c(50000,50000), square = F) %>% st_as_sf() %>% mutate(ID = row_number())
 
grid <- grid %>% 
  st_join(stability %>% inner_join(hl) %>% st_as_sf() %>% st_transform(st_crs(usa)), left = F)

Modes <- function(x) {
  ux <- unique(x)
  tab <- tabulate(match(x, ux))
  ux[tab == max(tab)][1]
}

modalStability <- grid %>% st_set_geometry(NULL) %>%
  group_by(ID) %>%
  summarise(modalStab = Modes(totalChange))

gridStab <- grid %>% inner_join(modalStability) %>% mutate(modalStab = as.numeric(as.character(modalStab)))

ggplot() + 
  geom_sf(data = usa) +
  geom_sf(data = gridStab, aes(fill = modalStab)) +
  scale_fill_gradient(low = '#108dc7', high = '#ef8e38', 'Mode State \nChanges') +
  ggthemes::theme_map(base_size = 11) +
  theme(legend.position = 'bottom')

ggsave('figures/ModalStateChange.png', width = 3.5, height = 3.5, units = 'in')

```

## Distance frequency distribution to examine spatial autocorrelation

```{r}
## Figure out modal cluster for each lake

Modes <- function(x) {
  ux <- unique(x)
  tab <- tabulate(match(x, ux))
  ux[tab == max(tab)][1]
}


modeCluster <- clusters.sf %>% st_set_geometry(NULL) %>%
  group_by(Hylak_id) %>%
  summarise(modeClust = Modes(cluster))

## Turn it into meters so we can make sense of distances, we'll use albers equal area USA
modeCluster <- modeCluster %>% 
  inner_join(
    clusters.sf %>% distinct(Hylak_id, .keep_all = T) %>% select(Hylak_id, Pour_lat, Pour_long)) %>% st_as_sf() %>%
  st_transform(102003)


## Calculate the distance matrix between a sample of clusters (this scales exponentially, so
## memory limits get hit pretty fast)
disSamp <- modeCluster %>% sample_frac(.3)
dist <- st_distance(disSamp)/1000
dist[lower.tri(dist, diag = T)] <- NA
units(dist) <- NULL
dist <- round(dist)

## Create a matrix where each cell shows if to clusters are the same or different from each other
sameV <- matrix(data = rep(disSamp$modeClust, nrow(disSamp)), ncol = nrow(disSamp))
sameH <- matrix(data = rep(disSamp$modeClust, nrow(disSamp)), nrow = nrow(disSamp), byrow = T)
same <- sameV == sameH
same[lower.tri(same, diag = T)] <- NA
rm(sameH,sameV)

## Map a function calculating the frequency of same cluster pairs to different cluster pairs
## at 50km intervals
window <- seq(50,4000,50)

spatialSim <- function(distance){
  tibble(distance = distance, 
         sameC = sum(same & dist <= distance & dist > distance -50, na.rm = T),
         difC = sum(!same & dist <= distance & dist > distance -50, na.rm = T))
}

distFreq <- window %>% map_dfr(spatialSim)

write_feather(distFreq, paste0('data/out/',iterations,'_DistanceFrequencies.feather'))

## Get rid of stuff cause it takes up a lot of memory
rm(same, dist, disSamp)

distFreq %>% pivot_longer(-distance, names_to = 'ClustType', values_to = 'Frequency') %>%
  ggplot(aes(x = distance, y = Frequency, color = ClustType)) + geom_point()

distFreq %>% mutate(pSame = sameC/(sameC + difC), pDif = difC/(sameC + difC)) %>% select(-sameC, -difC) %>%
  pivot_longer(-distance, names_to = 'ClustType', values_to = 'Frequency') %>%
  ggplot(aes(x = distance, y = Frequency, color = ClustType)) + geom_col(position = 'stack') +
  geom_hline(aes(yintercept = .2))
```


## Look at spatial distribution of clusters

```{r}
## By period
usa <- maps::map('usa', plot = F) %>% st_as_sf() %>% st_transform(102003) 

grid <- st_make_grid(usa, cellsize = c(100000,100000), square = F) %>% st_as_sf() %>% mutate(ID = row_number())
 
grid <- grid %>% st_join(clusters.sf %>% st_transform(102003), left = F)

Modes <- function(x) {
  ux <- unique(x)
  tab <- tabulate(match(x, ux))
  mode <- ux[tab == max(tab)]
  mode <- ifelse(length(mode) > 1, 'Mixed', as.character(mode))
  return(as.factor(mode))
}

modalClust <- grid %>% st_set_geometry(NULL) %>%
  group_by(ID, period) %>%
  summarise(modalClust = Modes(cluster)) %>%
  mutate(modalClust = factor(modalClust, levels = c('Spring Greening', 'Summer Greening', 'Bimodal', 'Aseasonal (Blue)', 'Aseasonal (Green)', 'Mixed')))

grid <- grid %>% inner_join(modalClust)

p1 <- ggplot() + 
  geom_sf(data = usa, fill = 'grey90') +
  geom_sf(data = grid, aes(fill = modalClust)) +
  scale_fill_manual(values = c("#5F4690", "#EDAD08","#CC503E", "#38A6A5", "#73AF48", 'grey70'))+
  ggthemes::theme_map(base_size = 12) +
  labs(fill = 'Modal Cluster', tag = 'a') +
  theme(legend.position = 'top') +
  facet_wrap(~period)

p1
ggsave('figures/DominantCluster.png', width = 6.5, height = 4, units = 'in')


## They look real bad together
p2 <- distFreq %>% mutate(pDif = difC/(sameC + difC), pSame = 1-pDif) %>% 
  ggplot(aes(x = distance)) + 
  geom_ribbon(aes(ymax = pSame, ymin = 0, fill = 'Same\nCluster')) +
  geom_ribbon(aes(ymax = 1, ymin = pSame, fill = 'Different\nCluster')) +
  scale_fill_viridis_d(option = 'plasma', end =.6) +
  coord_cartesian(xlim = c(50,3000)) +
  geom_segment(aes(x = 50, xend = 3500, y  = .2, yend = .2, color = 'Expected\nRandom\nFrequency'), size = 1, linetype = 4) +
  scale_color_manual(values = 'black') +
  theme_classic() +
  theme(legend.title = element_blank(),
        legend.direction = 'horizontal',
        legend.box = 'horizontal') +
  scale_y_continuous(labels = scales::percent) +
  scale_x_continuous(trans = 'log10') +
  labs(x= 'Distance Between Lakes (km)', y = 'Frequency', fill = 'Cluster Type', tag = 'b')


g <- gridExtra::grid.arrange(p1,p2, nrow = 2, heights = c(6,1.5))

ggsave('figures/DominantClusterPlusFreq_v2.png', plot = g, width = 6.5, height = 6, units = 'in')


summary(lm(pSame~distance, data = distFreq %>% mutate(pSame = sameC/(sameC+difC))))
```