diff --git a/.gitignore b/.gitignore
index d714431..1d73627 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,3 +3,6 @@
 .Rhistory
 *.Rproj
 .DS_Store
+Data/
+Figures/
+Supp Figures/
\ No newline at end of file
diff --git a/16dayRecursion_HomeRange_Script.R b/16dayRecursion_HomeRange_Script.R
index 248c631..da29a35 100644
--- a/16dayRecursion_HomeRange_Script.R
+++ b/16dayRecursion_HomeRange_Script.R
@@ -120,7 +120,7 @@ polymeta <- names(tumaps_sf) %>%
   rename(id = V1, year = V2, date_id = V3) %>%
   mutate(
     file = paste0( # match processed MODIS file names
-      "Processed/", year,
+      "~/Dropbox/Processed/", year,
       formatC(as.numeric(date_id), width = 3, flag = "0"), ".tif"
     ),
     list_name = names(tumaps_sf)
@@ -133,22 +133,11 @@ ordered_files <- polymeta[match(names(tumaps_sf), polymeta$list_name), ]
 # names(tumaps_sf) == ordered_files$list_name # yay! the same!
 
 plan(multiprocess)
-extracted_df <- future_pmap_dfr(
+biweekly <- future_pmap_dfr(
   list(polys = tumaps_sf, file = ordered_files$file),
   extract_NDVI
 )
 
-# Etosha Salt Pan results in missing values in some satellite imagery: set these to -10000
-# and handle appropriately for summary statistics
-biweekly <- extracted_df %>%
-  mutate(
-    meanNDVI = ifelse(is.na(meanNDVI), -10000, meanNDVI),
-    sdNDVI = ifelse(is.na(sdNDVI), 0, meanNDVI),
-    maxNDVI = ifelse(is.na(maxNDVI), -10000, maxNDVI),
-    minNDVI = ifelse(is.na(minNDVI), -10000, minNDVI),
-    medianNDVI = ifelse(is.na(medianNDVI), -10000, medianNDVI)
-  )
-
 ### Summary Statistics
 biweekly %>% select(id, date_id) %>% distinct() %>% nrow()
 biweekly %>% select(id) %>% distinct() %>% nrow()
@@ -158,7 +147,8 @@ biweekly %>%
   summarise(
     mean_nsv = mean(nsv.43200), sd_nsv = sd(nsv.43200),
     mean_mlnv = mean(mnlv.43200), sd_mnlv = sd(mnlv.43200),
-    mean_ndvi = mean(meanNDVI * .0001), sd_NDVI = sd(sdNDVI * .0001)
+    mean_ndvi = mean(meanNDVI * .0001, na.rm = T), 
+    sd_NDVI = sd(sdNDVI * .0001, na.rm = T)
   )
 
 ########  16 day Home Range constuction & NDVI Analysis ########
@@ -226,7 +216,7 @@ polymeta <- names(polys) %>%
     interval = paste0( year,
     formatC(as.numeric(date_id), width = 3, flag = "0")),
     file = paste0(
-      "Processed/", year,
+      "~/Dropbox/Processed/", year,
       formatC(as.numeric(date_id), width = 3, flag = "0"), ".tif"
     ),
     list_name = names(polys)
@@ -239,14 +229,14 @@ extracted_UDs <- pmap_df(list(polys = polys_by_interval,
                              file = unique(polymeta$file)), 
                         extract_NDVI) 
 
-scaled <- extracted_UDs %>% 
+scaledk <- extracted_UDs %>% 
   mutate(iso.level = as.factor(iso.level), 
          date_id = as.character(date_id),
          meanNDVI = meanNDVI *.0001,
          logarea = log(area))%>% 
   rename(met = method)
 
-scaled_k90 <- scaled %>% filter(met == "klocoh", iso.level == 90)
+scaled_k90 <- scaledk %>% filter(met == "klocoh", iso.level == 90)
 
 GLM <- gls(logarea ~ meanNDVI, data = scaled_k90, method = "ML", na.action = na.omit)
 lmm <- lme(logarea ~ meanNDVI, data = scaled_k90,
diff --git a/24hrDisplacement_Script.R b/24hrDisplacement_Script.R
index 561a5e6..1f1a4fd 100644
--- a/24hrDisplacement_Script.R
+++ b/24hrDisplacement_Script.R
@@ -7,7 +7,12 @@ source("DataParsing_Script.R")
 library(caret)
 
 ## Time of Day displacement analysis
-ToDfilter <- all_rhinos  %>% st_set_geometry(NULL) %>% 
+sampling_days <- c(seq.Date(from = ymd("2011-10-19"), to = ymd("2014-01-24"), 1), 
+                   seq.Date(from = ymd("2017-06-20"), to = ymd("2018-11-29"), 1))
+intervals <- interval(start = ymd_hm(paste(sampling_days - 1, "22:30")), 
+                      end = ymd_hm(paste(sampling_days, "22:30")))
+
+ToDfilter <- all_rhinos %>% st_set_geometry(NULL) %>% 
   mutate(rnd_hour = hour(round_date(date, unit = "1 hour")), 
          ToD = case_when(
            rnd_hour %in% c(23, 0:3) ~ "mid-night",
@@ -23,29 +28,43 @@ ToDfilter <- all_rhinos  %>% st_set_geometry(NULL) %>%
            rnd_hour %in% 17:21 ~ 19,
            TRUE ~ 0
          ),
-         offset_hour = ymd_hms(paste(local_date, local_time)) - 
-           ymd_h(paste(local_date, rnd_hour)),
-         offset_ToD = 
-           ymd_hms(paste(local_date, local_time)) - 
-           ymd_h(paste(local_date, ToD_est))) %>% 
+         offset_hour = case_when(
+           rnd_hour == 0 & hour(date) == 23 ~ difftime(ymd_hms(paste(local_date, local_time)), 
+                                                       ymd_h(paste(local_date + 1, rnd_hour)), 
+                                                       units = "secs"), 
+           TRUE ~ difftime(ymd_hms(paste(local_date, local_time)),
+                           ymd_h(paste(local_date, rnd_hour)),
+                           units = "secs") 
+         ),
+         offset_ToD = case_when(
+           rnd_hour %in% c(0, 23) & hour(date) != 0 ~ difftime(ymd_hms(paste(local_date, local_time)), 
+                                                               ymd_h(paste(local_date + 1, ToD_est)), 
+                                                               units = "secs"), 
+           TRUE ~ difftime(ymd_hms(paste(local_date, local_time)),
+                           ymd_h(paste(local_date, ToD_est)),
+                           units = "secs") 
+         )) %>% 
   filter(ToD != "other") %>% 
-  group_by(id, local_date, ToD) %>% 
+  mutate(interval = map_dbl(date, ~which(.x %within% intervals))) %>% 
+  group_by(id, interval, ToD) %>%  
   mutate(closest =  abs(as.numeric(offset_ToD)) == min(abs(as.numeric(offset_ToD)))) %>% 
   filter(closest) %>% ungroup() %>% 
   mutate(ToD = readr::parse_factor(ToD, ordered = T, levels = c("mid-night", "dawn", "mid-day", "dusk")))
 
 ToD <- ToDfilter %>% split(.$id) %>% 
-  map(., ~.x %>% 
-        complete(local_date = seq.Date(from = min(local_date, na.rm = T), to = max(local_date, na.rm = T), 1)) %>% 
-        complete(ToD, nesting(local_date)) %>% filter(!is.na(ToD)) %>% 
-        arrange(local_date, ToD)) 
+  map(., ~.x %>%
+        complete(interval = seq(from = min(interval, na.rm = T), to = max(interval, na.rm = T), 1)) %>% 
+        complete(ToD, nesting(interval)) %>% filter(!is.na(ToD)) %>% 
+        fill(id) %>% 
+        arrange(interval, ToD)) 
+
+# test that the ordering happened correctly, crucial to next step:
+sum(map_lgl(ToD,~ all(.x$ToD == rep_len(levels(ToDfilter$ToD), length(.x$ToD))))) == 59
 
 lag4 <- ToD %>%
-  keep(., function(x) {
-    na.omit(unique(x$id)) != "SAT645"
-  }) %>% # drop 645, too short, throws error.
   map_df(., function(x) {
-    mutate(x,
+    
+    df <- mutate(x,
            dx = c(diff(x), NA),
            dy = c(diff(y), NA),
            dt = sqrt(dx^2 + dy^2),
@@ -54,11 +73,38 @@ lag4 <- ToD %>%
            dt2 = sqrt(dx2^2 + dy2^2),
            dx4 = c(diff(x, lag = 4), NA, NA, NA, NA),
            dy4 = c(diff(y, lag = 4), NA, NA, NA, NA),
-           dt4 = sqrt(dx4^2 + dy4^2),
-           n = n()
+           dt4 = sqrt(dx4^2 + dy4^2)
     )
+    
+    # extracting true time between points
+    # being explicit about units is important. 
+    s <- diff(df$date)
+    s2 <- diff(df$date, lag = 2)
+    s4 <- diff(df$date, lag = 4)  
+    units(s) <- "hours"
+    units(s2) <- "hours"
+    units(s4) <- "hours"
+    
+    df$time <- c(as.vector(s), NA)
+    df$time2 <- c(as.vector(s2), NA, NA)
+    df$time4 <- c(as.vector(s4), NA, NA, NA, NA)
+    
+    return(df)
   })
 
+### make sure time diffferences are as expected:
+mean(lag4$time, na.rm = T)
+range(lag4$time, na.rm = T)
+sd(lag4$time, na.rm = T)
+
+mean(lag4$time2, na.rm = T)
+sd(lag4$time2, na.rm = T)
+range(lag4$time2, na.rm = T)
+
+mean(lag4$time4, na.rm = T)
+range(lag4$time4, na.rm = T)
+sd(lag4$time4, na.rm = T)
+
 # identify individuals with consistent 12 hr fixes 
 con_12hr <- map(ToD, function(df) {
   df <- filter(df, ToD %in% c("dawn", "dusk"))
@@ -69,31 +115,85 @@ con_12hr <- map(ToD, function(df) {
 ids <- names(map_dbl(con_12hr, nrow)[map_dbl(con_12hr, nrow) > 60])
 
 chosen9 <- con_12hr[ids] %>%
-  map_df(., ~ mutate(.x,
+  map_df(., function(x){
+    
+    df <- mutate(x,
                      dx = c(diff(x), NA),
                      dy = c(diff(y), NA),
                      dt = sqrt(dx^2 + dy^2),
                      dx2 = c(diff(x, lag = 2), NA, NA),
                      dy2 = c(diff(y, lag = 2), NA, NA),
                      dt2 = sqrt(dx2^2 + dy2^2),
-                     n = n()
-  )) %>% 
+                     n = n())
+    
+    s <- diff(df$date)
+    s2 <- diff(df$date, lag = 2)
+    units(s) <- "hours"
+    units(s2) <- "hours"
+    
+    df$time <- c(as.vector(s), NA)
+    df$time2 <- c(as.vector(s2), NA, NA)
+    
+    return(df)
+  }) %>% 
   group_by(id) %>% 
-  mutate(index = seq(1, n(), 1))
+  mutate(index = seq(1, n(), 1), 
+         disphr = dt2/time2)
+
+### Periodicity test!
+chosen9 %>% split(.$id) %>% 
+  map_dbl(nrow)
+# a couple have 105-155 relocations but most < 100
+
+chosen9 %>% split(.$id) %>% 
+  map(function(x){ptest::ptestg(na.omit(x$disphr))})
+# rejected the null hypothesis of periodicity (@ alpha = .05) for 1 individual SAT2372.  
 
-## KS test
-collapsed <- filter(lag4, !is.na(dt4)) %>%
-  mutate(ToD2 = forcats::fct_collapse(ToD, dawn = "dawn", group_other = T))
-classes <- collapsed %>% dplyr::select(ToD2, dt4) %$% downSample(dt4, ToD2)
-dawn_class <- filter(classes, Class == "dawn") %>% pull(x)
-other_class <- filter(classes, Class == "Other") %>% pull(x)
-ks.test(dawn_class, other_class)
+# Pairwise:
+#dumbest way i can figure this out tonight... 
+x <- expand.grid(1:4, 1:4) %>% filter(Var1 != Var2)
+x[x[,2]>=x[,1],] 
+
+comb <- expand.grid(unique(levels(lag4$ToD)), unique(levels(lag4$ToD))) %>% 
+  filter(Var1 != Var2) %>% .[x[,2]>=x[,1],] %>% mutate_all(as.character())
+
+
+ToD_disp <- filter(lag4, !is.na(dt4)) %>%
+  mutate(disphr = dt4/time4) %>% 
+  select(ToD, disphr) %>% 
+  split(.$ToD) %>%
+  map(~pull(.x, disphr))
+
+results <- pmap_dbl(list(a = comb$Var1, b = comb$Var2), 
+                    function(a,b){
+                      suppressWarnings(ks.test(ToD_disp[[a]], ToD_disp[[b]])$p.value)
+                    }) %>% p.adjust(method = "bonferroni") 
+
+names(results) <- paste(comb$Var1, comb$Var2, sep =".")
+
+results %>% round(3) %>% knitr::kable("latex")
 
 ## Time of day points near water?
-buffers <- st_buffer(water_utm, 500)
-water_pts <- st_intersection(buffers, 
-                             st_as_sf(lag4, 
-                                      coords = c("x", "y"), 
-                                      crs = 32733, na.fail = F))
-water_pts %>% select(-n) %>% group_by(ToD) %>% tally()
+water_sens <- function(dist) {
+  buffers <- st_buffer(water_utm, dist)
+  water_pts <- st_intersection(buffers, 
+                               st_as_sf(lag4, 
+                                        coords = c("x", "y"), 
+                                        crs = 32733, na.fail = F))
+  water_pts %>% st_set_geometry(NULL) %>% group_by(ToD) %>% tally()
+}
+
 
+Cairo::CairoPNG(filename = "waterbuffer.png", width = 1200, height = 1000,
+                res = 180)
+map_df(c("100" = 100, "250" = 250, "500" = 500, "750" = 750,  "1000" = 1000), 
+       water_sens, .id = "dist") %>% 
+  ggplot(., aes(x = as.numeric(dist), y = n, color = ToD)) + 
+  geom_point() +
+  geom_path() + 
+  theme_minimal() +
+  labs(x = "Buffer distance (m) from watering hole", 
+       y = "Number of fixes",
+       title = "Influence of time of day on distance to watering hole") +
+  theme(plot.title = element_text(hjust = .5))
+dev.off()
diff --git a/AnnualRecursion_Script.R b/AnnualRecursion_Script.R
index 0235438..4bd8ff2 100644
--- a/AnnualRecursion_Script.R
+++ b/AnnualRecursion_Script.R
@@ -82,10 +82,11 @@ tumap_stats <- tibble(id = names(revis_prop), revis_prop, visited, avg_revis, sd
 
 cor(tumap_stats[,-1])
 
-cor(revis_prop, visited)
-cor(revis_prop[-6], visited[-6])  # 280 is an outlier... not sure why. 
+cor.test(revis_prop, visited)
+cor.test(revis_prop[-6], visited[-6])  # 280 is an outlier... not sure why. 
 
-cor(med_revis, visited)
+cor.test(med_revis, visited)
+cor.test(sd_revis, visited)
 
 # range size has a negative relationship with proportion revisted cells
 # i.e. smaller ranges are likely to revisited more of their range. 
@@ -98,4 +99,11 @@ top25rev <- map_dbl(tumaps_sf, ~filter(.x, nsv.604800 > 0)  %>%
 
 # the proportion of cells visited that gain the top 25% of returns, 
 # is fairly equal across individuals, and is not strongly correlated with range size. 
-cor(top25rev, visited)
\ No newline at end of file
+top25rev/visited 
+cor(top25rev, visited)
+
+
+stats <- tibble(visited, revisited,  top25rev, revis_prop, 
+                avg_revis, sd_revis, med_revis, id = names(visited))
+
+summary(lm(visited ~ revis_prop, data = stats))
diff --git a/DataParsing_Script.R b/DataParsing_Script.R
index fcf7941..63eb305 100644
--- a/DataParsing_Script.R
+++ b/DataParsing_Script.R
@@ -12,14 +12,14 @@ conflicted::conflict_prefer("filter", "dplyr")
 
 # Metadata regarding file index and location of MODIS imagery
 # MODIS imagery is freely available from NASA.gov
-filemetadata <- read_csv("MODIS_metadata.csv")
+filemetadata <- read_csv("Data/MODIS_metadata.csv")
 
 # downloaded from http://ecoregions2017.appspot.com/, Licensed under CC-BY 4.0
-ecoregions <- read_sf("Ecoregions2017.shp")
+ecoregions <- read_sf("Data/Ecoregions2017.shp")
 
 # functional water data - available upon reasonable request
 # not public in order to protect endangered species
-water_utm <- st_read("functional water.shp", crs = 4326) %>%
+water_utm <- st_read("Data/functional water.shp", crs = 4326) %>%
   st_transform(32733)
 
 # rhino data - available upon reasonable request
diff --git a/Figures_Script.R b/Figures_Script.R
index 032e5fc..a86cfb5 100644
--- a/Figures_Script.R
+++ b/Figures_Script.R
@@ -3,12 +3,11 @@
 ### (submitted to Movement Ecology, BioMove Special Edition)
 ### Author: Dana Paige Seidel
 
-####### Necessary Libraries#######
+####### Necessary Libraries #######
 library(stmove) # package not available on CRAN, see github.com/dpseidel/stmove
 library(ggsn) # for spatial scale and coordinate arrow
 
 ########### Source Analysis Scripts ##########
-source("DataParsing_Script.R") # defines all_rhinos, namib_eco
 source("24hrDisplacement_Script.R") # defines: chosen9, lag4, ToD
 source("16dayRecursion_HomeRange_Script.R") # defines: biweekly, tumaps
 source("AnnualRecursion_Script.R") #  defines: df_3visits
@@ -42,6 +41,8 @@ ggplot() +
 stmove::plot_timeline(all_rhinos)
 
 #### Figure 3 - Dawn displacement time series
+Cairo::CairoPNG(filename = "Figures/distplacement.png", width = 1500, height = 1500,
+                res = 180)
 chosen9 %>%
   ungroup() %>%
   filter(ToD == "dawn") %>%
@@ -58,21 +59,25 @@ chosen9 %>%
     nudge = sd(yday(date))
   ) %>%
   ggplot() +
-  geom_path(aes(dt2, x = yday(date), color = sex), na.rm = T) +
+  geom_path(aes(disphr, x = yday(date), color = sex), na.rm = T) +
   scale_color_manual(values = c("black", "#000080")) +
   # geom_text(aes(x = x_lab + nudge/3, y = 15500,
   #              label = paste(round(n/2), "days")), size = 3) +
-  coord_cartesian(ylim = c(0, 16000)) +
+  coord_cartesian(ylim = c(0, 700)) +
   facet_wrap(~id_f, scales = "free_x", ncol = 3) +
   theme_minimal() +
   labs(
-    title = "Dawn-Dawn 24-displacement through time",
-    x = "Julian day", y = " 24-hr displacement (m)",
-    caption = "Note: A single observation approaching 30 km displacement was cut off from SAT277's graph to maintain comparable scales"
+    title = "Dawn-Dawn standardized 24-hr displacement through time",
+    x = "Julian day", y = "24-hr displacement per hour (m/hr)",
+    caption = "Note: A single observation approaching 1250 m/hr displacement was cut off from SAT277's graph to maintain comparable scales"
   ) +
   theme(plot.title = element_text(hjust = .5), legend.position = "none", panel.spacing = unit(1, "lines"))
+dev.off()
 
 #### Figure 4. Facet Distributions -- distributions
+Cairo::CairoPNG(filename = "Figures/distributions.png", width = 1500, height = 1500,
+                res = 180)
+
 chosen9 %>%
   ungroup() %>%
   filter(ToD == "dawn") %>%
@@ -89,36 +94,40 @@ chosen9 %>%
     nudge = sd(yday(date))
   ) %>%
   ggplot() +
-  geom_density(aes(dt2, color = sex), na.rm = T) +
+  geom_density(aes(disphr, color = sex), na.rm = T) +
   scale_color_manual(values = c("black", "#000080")) +
   geom_text(aes(
-    x = 25000, y = .0006,
+    x = 900, y = .0125,
     label = paste(round(n / 2), "days")
   ), size = 3) +
   facet_wrap(~id_f, ncol = 3) +
   theme_minimal() +
   labs(
-    title = "Dawn-Dawn 24-displacement through time",
-    x = " 24-hr displacement (m)"
+    title = "Dawn-Dawn standardized 24-hr displacement through time",
+    x = "24-hr displacement per hour (m/hr)"
   ) +
   theme(plot.title = element_text(hjust = .5), legend.position = "none", panel.spacing = unit(1, "lines"),
         plot.margin = margin(10, 15, 10, 10))
 dev.off()
 
+
 ####  Fig. 5
+Cairo::CairoPNG(filename = "Figures/density.png", width = 1500, height = 1500,
+                res = 180)
 ggplot(lag4) +
-  geom_density(aes(x = dt4, color = ToD, fill = ToD), alpha = .1, na.rm = T) +
+  geom_density(aes(x = dt4/time4, color = ToD, fill = ToD), alpha = .1, na.rm = T) +
   theme_minimal() +
   theme(legend.position = "top", plot.title = element_text(hjust = .5)) +
-  labs(title = "24-hr displacement distributions across time of day", x = "24-hr displacement (m)") +
+  labs(title = "Standardized 24-hr displacement distributions across time of day", 
+       x = "24-hr displacement per hour (m/hr)") +
   scale_color_manual("Time of Day", values = c("#26384f", "#57beff", "#4f8a99", "#00578b")) +
   scale_fill_manual("Time of Day", values = c("#26384f", "#57beff", "#4f8a99", "#00578b"))
+dev.off()
 
 #### Figure 6.
 # T-Locoh  TimeUse Plot.
 par(mfrow = c(1, 2))
 plot(tumaps$`SAT189-2011-305`$result)
-dev.off()
 
 #### Figure 7
 filter(biweekly, nsv.43200 > 1) %>%
@@ -140,6 +149,7 @@ filter(biweekly, nsv.43200 > 1) %>%
   ) +
   theme_minimal() +
   theme(plot.title = element_text(hjust = .5))
+dev.off()
 
 #### Figure 8:
 ggplot(filter(df_3visits, sv == T) %>% group_by(id.x) %>% 
@@ -151,5 +161,5 @@ ggplot(filter(df_3visits, sv == T) %>% group_by(id.x) %>%
                     values = colorRampPalette(c("#08306b", "#9ecae1"))(9)) +
   theme_minimal() +
   labs(x = "Time between visits (days)", 
-       title = "Annual recursion: time to to return across individuals") +
+       title = "Annual recursion: time to return across individuals") +
   theme(plot.title = element_text(hjust = .5))
diff --git a/README.md b/README.md
index 3d15392..4f0e470 100755
--- a/README.md
+++ b/README.md
@@ -6,14 +6,15 @@ Mesoscale movement and recursion behaviors of Namibian black rhinos
 This repository contains the source code of analyses contained within my 
 manuscript with Wayne Linklater, Werner Kilian, Pierre du Preez, and Wayne M. Getz, 
 submitted for consideration to Movement Ecology's BioMove Special Edition in
-May 2019. 
+May 2019. Revisions were completed in July 2019, introducing small changes to
+the 24-hr displacement estimation and adding several new supplementary analyses. 
 
 - Author: Dana Seidel 
 - License: [CC-BY 4.0](https://creativecommons.org/licenses/by/4.0/)
 
 ## A note on data availability
 This repository includes cleaned analysis scripts but no associated data and is
-therefore *not* completely reproducible.The datasets generated and/or analyzed 
+therefore *not* completely reproducible. The datasets generated and/or analyzed 
 during this study are not publicly available in order to protect the locations 
 and resources of a critically endangered species but may be made privately 
 available by me, the corresponding author, on reasonable request. The MODIS 
diff --git a/SupplementaryAnalyses/BiweekelyNullModel.R b/SupplementaryAnalyses/BiweekelyNullModel.R
new file mode 100644
index 0000000..c1898fc
--- /dev/null
+++ b/SupplementaryAnalyses/BiweekelyNullModel.R
@@ -0,0 +1,211 @@
+##### Null Distribution of NDVI -- for 16 day revisitation grids?
+
+### randomly sample -- the same number of points from each NDVI file
+### and compare to the biweekly extracted points. 
+
+### do we want to randomly sample within a home range?
+set.seed(4236)
+
+biweekly_null <- biweekly %>% 
+  group_by(id, year, date_id) %>% 
+  filter(nsv.43200 > 0, !is.na(meanNDVI)) %>% 
+  mutate(nullNDVI_idint = sample(meanNDVI, replace = F)) %>% 
+  group_by(id) %>% 
+  mutate(nullNDVI_id = sample(meanNDVI, replace = F)) %>% 
+  ungroup() %>% 
+  mutate(nullNDVI = sample(meanNDVI, replace = F))
+
+
+Cairo::CairoPNG(filename = "ndvidist.png", width = 1200, height = 1000,
+                res = 180)
+ggplot() + 
+  geom_histogram(aes(x = meanNDVI * .0001, 
+                     fill = ifelse(nsv.43200 > 0, "visited", "not-visited" )), 
+                 data = biweekly, bins = 100) + 
+  scale_fill_manual("Visitation", values = c("#696969", "#4169e1")) +
+  theme_minimal() + labs(
+    title = "Distribution of extracted NDVI values", 
+    x = "Average NDVI within extracted 1km^2 grid cell"
+  )  +
+  theme(plot.title = element_text(hjust = .5))
+dev.off()
+
+library(patchwork)
+p_org = filter(biweekly_null, nsv.43200 > 1) %>%
+  ggplot(aes(y = nsv.43200, x = meanNDVI * .0001)) +
+   stat_bin2d(aes(fill = log(stat(count))), binwidth = c(.05, 1)) +
+  # # stat_bin2d(aes(color = ..count..), size = 1, bins = 30) +
+   scale_fill_gradient("Count",
+                       low = "#c6dbef", high = "#08306b",
+                       labels = function(x) {
+                         round(exp(1)^x)
+                       }
+   ) +
+  # geom_jitter(size = .6, alpha = .3) +
+  #geom_smooth(method = "glm") +
+  #scale_x_continuous(limits = c(-0.149, 0.725)) +
+  # scale_color_viridis_c(option="inferno")+
+  labs(
+    y = "Number of Separate Visits",
+    x = "Average NDVI",
+    title = "Null Model -- all visited cells"
+  ) +
+  theme_minimal() +
+  theme(plot.title = element_text(hjust = .5))
+
+filter(biweekly_null, nsv.43200 > 1) %>%
+  ggplot(aes(y = nsv.43200, x = nullNDVI_id * .0001)) +
+  stat_bin2d(aes(fill = log(stat(count))), binwidth = c(.05, 1)) +
+  # stat_bin2d(aes(color = ..count..), size = 1, bins = 30) +
+  scale_fill_gradient("Count",
+                      low = "#c6dbef", high = "#08306b",
+                      labels = function(x) {
+                        round(exp(1)^x)
+                      }
+  ) +
+  # geom_jitter(size = .6, alpha = .3) +
+  #geom_smooth(method = "glm", color = "red") +
+  # scale_color_viridis_c(option="inferno")+
+  labs(
+    y = "Number of Separate Visits",
+    x = "Average NDVI",
+    title = "Null Model -- ID specific"
+  ) +
+  theme_minimal() +
+  theme(plot.title = element_text(hjust = .5))
+
+
+filter(biweekly_null, nsv.43200 > 1) %>%
+  ggplot(aes(y = nsv.43200, x = nullNDVI_id * .0001)) +
+   geom_jitter(size = .6, alpha = .3) +
+  stat_smooth(method = "glm", formula = y ~ poly(x, 2), size = 1) +
+  geom_smooth(method = "glm") +
+  # scale_color_viridis_c(option="inferno")+
+  labs(
+    y = "Number of Separate Visits",
+    x = "Average NDVI",
+    title = "Null Model - per individual"
+  ) +
+  theme_minimal() +
+  theme(plot.title = element_text(hjust = .5))
+
+lm <- lm(`nsv.43200` ~ (meanNDVI*.0001),  data = biweekly_null)
+quad <- lm(nsv.43200 ~ poly(meanNDVI*.0001, 2),  data = biweekly_null)
+
+anova(lm, quad)
+
+filter(biweekly_nullrev, nsv.43200 > 1) %>%
+  ggplot(aes(y = nsv.43200, x = meanNDVI * .0001)) +
+  geom_jitter(size = .6, alpha = .3) +
+  stat_smooth(method = "lm", formula = y ~ poly(x, 2), size = 1)
+  #geom_smooth(method = "glm") +
+  #geom_smooth(aes(y = nsv.43200, x = nullNDVI * .0001), method = "glm", col = "red") +
+  # scale_color_viridis_c(option="inferno")+
+  labs(
+    y = "Number of Separate Visits",
+    x = "Average NDVI",
+    title = "Null Model - per individual per interval"
+  ) +
+  theme_minimal() +
+  theme(plot.title = element_text(hjust = .5))
+
+
+
+biweekly_nullrev %>% 
+  filter(nsv.43200 > 1) %>% 
+  summarise(avg = mean(meanNDVI, na.rm = T), sd = sd(meanNDVI, na.rm = T),
+            null = mean(nullNDVI, na.rm = T), nsd = sd(nullNDVI, na.rm = T)) %>% 
+clipr::write_clip()
+
+biweekly_nullid %>% ungroup() %>%  
+  filter(nsv.43200 > 1) %>% 
+  summarise(avg = mean(meanNDVI, na.rm = T), sd = sd(meanNDVI, na.rm = T),
+            null = mean(nullNDVI, na.rm = T), nsd = sd(nullNDVI, na.rm = T)) %>% 
+  clipr::write_clip()
+
+
+### If you want to find an intermediate value you have to test the fit of a
+### quadratic y=ax^2 + bx + c and see if this fit better from the AIC point of 
+### view than your linear fit.  
+
+mod_biweekly <- biweekly_null %>% 
+  mutate(meanNDVI = .0001*meanNDVI, 
+         nullNDVI_id = .0001*nullNDVI_id
+         )
+
+mod1 <- glm(nsv.43200 ~ nullNDVI_id, data = mod_biweekly, family = "poisson") 
+mod2 <- glm(nsv.43200 ~ meanNDVI, data = mod_biweekly, family = "poisson") 
+mod3 <- glm(nsv.43200 ~ poly(nullNDVI_id,2), data = mod_biweekly, family = "poisson") 
+mod4 <- glm(nsv.43200 ~ poly(meanNDVI, 2), data = mod_biweekly, family = "poisson") 
+
+anova(mod1, mod2, mod3, mod4)
+AIC(mod1, mod2, mod3, mod4)
+
+Cairo::CairoPNG(filename = "quad.png", width = 1200, height = 1000,
+                res = 180)
+mod_biweekly %>%
+  ggplot(aes(y = nsv.43200, x = meanNDVI * .0001)) +
+  geom_jitter(size = .6, alpha = .3) +
+  stat_smooth(method = "lm", formula = y ~ poly(x, 2), size = 1) +
+  labs(
+  y = "Number of Separate Visits",
+  x = "Average NDVI",
+  title = "Average NDVI values in visited cells, fitted with a Quadratic Model"
+) +
+  theme_minimal() +
+  theme(plot.title = element_text(hjust = .5))
+dev.off()
+
+
+### Replicate Fig 7 - Side by Side
+p <- filter(biweekly_null, nsv.43200 > 1) %>%
+  ggplot(aes(y = nsv.43200, x = meanNDVI * .0001)) +
+  stat_bin2d(aes(fill = log(stat(count))), binwidth = c(.05, 1)) +
+  # stat_bin2d(aes(color = ..count..), size = 1, bins = 30) +
+  scale_fill_gradient("Count",
+                      low = "#c6dbef", high = "#08306b",
+                      labels = function(x) {
+                        round(exp(1)^x)
+                      }
+  ) +
+  # geom_jitter(size = .6, alpha = .3) +
+  # scale_color_viridis_c(option="inferno")+
+  labs(
+    y = "Number of Separate Visits",
+    x = "Average NDVI",
+    title = "Relationship between 16-day revisitation and grid cell NDVI"
+  
+    ) +
+  theme_minimal() +
+  theme(plot.title = element_text(hjust = .5), 
+        legend.position = "None")
+
+p2 <- filter(biweekly_null, nsv.43200 > 1) %>%
+  ggplot(aes(y = nsv.43200, x = nullNDVI_id * .0001)) +
+  stat_bin2d(aes(fill = log(stat(count))), binwidth = c(.05, 1)) +
+  # stat_bin2d(aes(color = ..count..), size = 1, bins = 30) +
+  scale_fill_gradient("Count",
+                      low = "#c6dbef", high = "#08306b",
+                      labels = function(x) {
+                        round(exp(1)^x)
+                      }
+  ) +
+  scale_x_continuous(breaks = c(0,.25, .50, .75)) +
+  # geom_jitter(size = .6, alpha = .3) +
+  # scale_color_viridis_c(option="inferno")+
+  labs(
+    y = "",
+    x = "Average NDVI",
+    title = "Relationship between 16-day revisitation and grid cell NDVI", 
+    subtitle = "NULL MODEL") +
+  theme_minimal() +
+  theme(plot.title = element_text(hjust = .5), 
+        plot.subtitle = element_text(hjust =.5))
+
+library(patchwork)
+
+Cairo::CairoPNG(filename = "Fig7_SidebySide.png", width = 1200, height = 600,
+                res = 100)
+p + p2
+
+dev.off()
diff --git a/SupplementaryAnalyses/DisplacementPerHour.R b/SupplementaryAnalyses/DisplacementPerHour.R
new file mode 100644
index 0000000..7174a9b
--- /dev/null
+++ b/SupplementaryAnalyses/DisplacementPerHour.R
@@ -0,0 +1,12 @@
+Cairo::CairoPNG(filename = "stddisp.png", width = 1200, height = 1000,
+                res = 180)
+ggplot(lag4) +
+  geom_density(aes(x = dt4/time4, color = ToD, fill = ToD), alpha = .1, na.rm = T) +
+  theme_minimal() +
+  theme(legend.position = "top", plot.title = element_text(hjust = .5)) +
+  labs(title = "Displacement per hour distributions across time of day",
+       x = "24-hr displacement per hour (m/hr)") +
+  scale_color_manual("Time of Day", values = c("#26384f", "#57beff", "#4f8a99", "#00578b")) +
+  scale_fill_manual("Time of Day", values = c("#26384f", "#57beff", "#4f8a99", "#00578b"))
+dev.off()
+
diff --git a/SupplementaryAnalyses/HomerangeMethods.R b/SupplementaryAnalyses/HomerangeMethods.R
new file mode 100644
index 0000000..5254c8c
--- /dev/null
+++ b/SupplementaryAnalyses/HomerangeMethods.R
@@ -0,0 +1,138 @@
+##### LoCoh v. AKDE #####
+
+### To Compare the AKDE vs. k-LoCoH home range 
+
+# replacing the construct_UDs function in 16dayRecursion_HomeRange_Script
+# we can extract both types of home range at the same time. 
+proj4 = "+proj=utm +zone=33 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
+
+construct_UDs <- function(df) {
+  ## CTMM code
+  proj4 <- "+proj=utm +zone=33 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
+  telm <- stmove:::create_telemetry(df, proj4)
+  m.ouf <- ctmm::ctmm.guess(telm, interactive = FALSE) # automated model guess
+  M.OUF <- ctmm::ctmm.fit(telm, m.ouf) # this can take awhile...
+  UD <- ctmm::akde(telm, M.OUF, weights = F)
+  akde.isopolys <- ctmm::SpatialPolygonsDataFrame.UD(UD, level.UD = c(.9)) %>% st_as_sf()
+  
+  ## k-locoh Code
+  dropNA <- na.omit(df)
+  k <- round(sqrt(nrow(dropNA)))
+  lxy <- tlocoh::xyt.lxy(
+    xy = matrix(c(dropNA$x, dropNA$y), ncol = 2),
+    dt = dropNA$date, id = dropNA$id, tau.diff.max = 0, # disable filter to avoid errors
+    proj4string = sp::CRS(proj4),
+    status = F
+  )
+  lxy <- tlocoh::lxy.nn.add(lxy, s = 0, k = k, status = F)
+  lhs <- tlocoh::lxy.lhs(lxy,
+                         k = k, s = 0, iso.levels = c(0.90),
+                         iso.add = T, status = F
+  )
+  klocoh.isopolys <- tlocoh::isopleths(lhs)[[1]][1:3] %>% st_as_sf()
+
+  area_df <- tibble(
+    id = df$id[1],
+    date_id = paste0(year(df$date[1]), formatC(df$interval_start[1], width = 3, flag = "0")),
+    iso.level = 90,
+    method = c("akde", "klocoh")
+  ) %>%
+    mutate(
+      area = ifelse(method == "klocoh",
+                    klocoh.isopolys$area,
+                    st_area(akde.isopolys[str_detect(akde.isopolys$name, "ML"), ])
+      ),
+      edge = ifelse(method == "klocoh", klocoh.isopolys$edge.len,
+                    st_length(st_cast(akde.isopolys[str_detect(akde.isopolys$name, "ML"), ], "MULTILINESTRING"))
+      )
+    )
+  
+  return(list(
+    "AreaSummary" = area_df,
+    "AKDE" = akde.isopolys,
+    "klocoh" = klocoh.isopolys
+  ))
+}
+
+plan(multiprocess)
+polys <- future_map(dfs, safely(construct_UDs), .progress = T)
+
+join_constructs <- function(x){
+  akde <- x$result$AKDE[2,] %>% # ML 90%
+    dplyr::mutate(name = as.character(name),
+                  method = "akde")
+  klocoh <- x$result$klocoh %>% 
+    select(name = iso.level) %>% 
+    mutate(method = "klocoh")
+  
+  rbind(akde, klocoh) %>% 
+    dplyr::mutate(id = x$result$AreaSummary$id, 
+           date_id = x$result$AreaSummary$date_id, 
+           iso.level = x$result$AreaSummary$iso.level,
+           area = x$result$AreaSummary$area,
+           edge = x$result$AreaSummary$edge) %>% 
+    select(id, date_id, method, iso.level, area, edge)
+}
+
+polys_combined <- map(polys, join_constructs) %>% reduce(rbind)
+
+polys_by_interval <- polys_combined %>% split(.$date_id)
+
+polymeta <- names(polys) %>%
+  str_split("-", simplify = T) %>%
+  as_tibble() %>%
+  rename(id = V1, year = V2, date_id = V3) %>%
+  mutate(
+    interval = paste0( year,
+                       formatC(as.numeric(date_id), width = 3, flag = "0")),
+    file = paste0(
+      "~/Dropbox/Processed/", year,
+      formatC(as.numeric(date_id), width = 3, flag = "0"), ".tif"
+    ),
+    list_name = names(polys)
+  ) %>%
+  arrange(interval)
+
+unique(polymeta$interval) == names(polys_by_interval)
+
+extract_NDVI <- function(polys, file){
+  vlx <- velox(file)
+  polys <- st_cast(polys, "MULTIPOLYGON") # just to prevent multi-type frames
+  polys %>% 
+    st_set_geometry(NULL) %>% 
+    mutate(meanNDVI = 
+             vlx$extract(polys, fun = function(x) mean(x, na.rm = TRUE)),
+           sdNDVI =
+             vlx$extract(polys, fun = function(x) sd(x, na.rm = TRUE)),
+           maxNDVI = 
+             vlx$extract(polys, fun = function(x) max(x, na.rm = TRUE)),
+           minNDVI=
+             vlx$extract(polys, fun = function(x) min(x, na.rm = TRUE)),
+           medianNDVI = 
+             vlx$extract(polys, fun = function(x) median(x, na.rm = TRUE)))
+}
+
+extracted_UDs <- pmap_df(list(polys = polys_by_interval, 
+                              file = unique(polymeta$file)), 
+                         extract_NDVI) 
+
+scaled <- extracted_UDs %>% 
+  mutate(iso.level = as.factor(iso.level), 
+         date_id = as.character(date_id),
+         meanNDVI = meanNDVI *.0001,
+         logarea = log(area))%>% 
+  rename(met = method)
+
+
+scaled_k90 <- scaled %>% filter(met == "klocoh", iso.level == 90)
+scaled_a90 <- scaled %>% filter(met == "akde", iso.level == 90)
+
+
+lmm <- lme(logarea ~ meanNDVI, data = scaled_k90,
+           random = ~1|id, method = "ML",  na.action = na.omit)
+
+summary(lmm)
+
+lmm_a90 <- lme(logarea ~ meanNDVI, data = scaled_a90,
+           random = ~1|id, method = "ML",  na.action = na.omit)
+summary(lmm_a90)
diff --git a/SupplementaryAnalyses/SeasonalGroupEffects.R b/SupplementaryAnalyses/SeasonalGroupEffects.R
new file mode 100644
index 0000000..f5d1d56
--- /dev/null
+++ b/SupplementaryAnalyses/SeasonalGroupEffects.R
@@ -0,0 +1,161 @@
+##### Group effects, Seasonality
+
+
+### Seasonal Differences ###
+# wet season defined as november - april, according to rainfall data
+library(readxl)
+rainfall1981_2013 <- read_xlsx("Data/ENP_rainfall_daily_1981-2013_clean.xlsx",
+                               col_types = c(
+                                 "date", "numeric", "text", "numeric", # ignore the parsing errors
+                                 "numeric", "numeric", "numeric",
+                                 "numeric", "numeric", "numeric", "numeric", "numeric", "numeric"
+                               )
+)
+
+Cairo::CairoPNG(filename = "rainfall.png", width = 1200, height = 1000,
+                res = 180)
+rainfall1981_2013 %>%
+  gather(., area, rainfall, -month, -year, -month_cal, -day, -date) %>%
+  dplyr::group_by(month, year, area) %>%
+  dplyr::summarise(rainfall_month = sum(rainfall, na.rm = T)) %>%
+  group_by(month) %>%
+  summarise(mean_rainfall = mean(rainfall_month, na.rm = T)) %>%
+  ggplot(., aes(x = month, y = mean_rainfall)) + geom_bar(stat = "identity") +
+  ggtitle("Average Monthly Rainfall across all areas ENP, 1981-2013") +
+  theme_minimal() +
+  scale_x_continuous(breaks = c(1:12), 
+                     labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun",
+                                "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")) +
+  labs(y = "Average Monthly Rainfall (mm)") +
+  theme(plot.title = element_text(hjust = .5))
+dev.off()
+
+# wet == November - April 
+wet_chosen <- chosen9 %>% mutate(wet = ifelse(month(local_date) %in% c(1:4, 11:12), T, F))
+wet_lag <- lag4 %>% mutate(wet = ifelse(month(local_date) %in% c(1:4, 11:12), T, F))
+
+Cairo::CairoPNG(filename = "season9.png", width = 1200, height = 1000,
+                res = 180)
+ggplot(wet_chosen, aes(x = dt2, fill = wet)) +
+  geom_density(na.rm = T, alpha = .4) + 
+  labs(title = "24-hr Displacement Across Seasons - 9 Consistent Individuals", 
+       x = "24-hr displacement (m)") + 
+  scale_fill_discrete("Season", labels = c("dry", "wet")) +
+  facet_wrap(~ToD) +
+  theme_minimal() + 
+  theme(strip.background = element_rect(fill = "#F4F4F4"),
+        plot.title = element_text(hjust = .5))
+dev.off()
+
+wet_chosen %>% select(-n) %>% group_by(wet) %>% tally()
+
+Cairo::CairoPNG(filename = "season_TOD.png", width = 1200, height = 1000,
+                res = 180)
+ggplot(wet_lag, aes(x = dt4, fill = wet)) + 
+  geom_density(alpha = .4, na.rm = T) + facet_wrap(~ToD) +
+  theme_minimal() + 
+  labs(title = "24-hr Displacement Across Seasons", 
+       x = "24-hr displacement (m)") +
+  scale_fill_discrete("Season", labels = c("dry", "wet")) +
+  theme(strip.background = element_rect(fill = "#F4F4F4"),
+        plot.title = element_text(hjust = .5))
+dev.off()
+
+wet_lag %>% na.omit() %>% group_by(wet) %>% tally() 
+
+#### Demographics ####
+ids <- ToDfilter %>% pull(id) %>% unique()
+metadataDbbCss <- readxl::read_xlsx("../../../Dropbox/Rhino/Collar Info/Black rhino sat collars/Rhino Sat Collar Data/RhinoSAT_combined_editted.xlsx", sheet = 1)
+
+metadataDbb1718 <- readxl::read_xlsx("../../../Dropbox/Rhino/Rhino Movements/moredataetosha/ENPDbbmove2017_18.xlsx", sheet = 1)
+
+# bring it all together, sort out calf presence/pregnant. 
+old_collars <- metadataDbbCss %>% 
+  mutate(id = paste0(Type, `Collar ID`), 
+         id = ifelse(id == "SAT843", "SAT943", id)) %>%  # I have a hunch collar "943" is actually listed as collar "843" in this metadata. patching.. 
+  filter(id %in% ids) 
+
+Dbb1718_collars <- metadataDbb1718 %>% 
+  mutate(id = str_remove(str_remove(`Tag Number`, "IR-"), " ")) %>% 
+  filter(id %in% ids)
+
+nrow(old_collars) + nrow(Dbb1718_collars) == length(ids)
+ids %in% c(old_collars$id, Dbb1718_collars$id)
+
+collars_demographics <- old_collars %>% 
+  select(id, Sex, Composition, Age, Remarks) %>%
+  bind_rows(., Dbb1718_collars %>% 
+              select(id, Sex = Gender, Remarks = Notes)) %>% 
+  mutate(has_calf = case_when(
+    str_detect(Composition, "Size") ~ TRUE, 
+    str_detect(Composition, "calf") ~ TRUE, 
+    str_detect(Remarks,  "calf") ~ TRUE, 
+    TRUE ~ FALSE
+  ), 
+  pregnant = str_detect(Remarks, "pregnant"), 
+  pregnant = ifelse(is.na(pregnant), FALSE, pregnant), 
+  Sex = case_when(
+    str_detect(Sex, "^F") ~ "F", 
+    str_detect(Sex, "^M") ~ "M", 
+    TRUE ~ NA_character_
+  ), 
+  Condition = case_when(
+    has_calf ~  "with calf", 
+    pregnant ~ "with calf", 
+    !has_calf ~ "single")) 
+
+
+# When pregant females and females captured with calf are combine we have nearly even
+# groups. Note that one individual's, SAT2596, sex was not recorded and was dropped from
+# the analysis. 
+
+collars_demographics %>% filter(!is.na(Sex)) %>% 
+  group_by(Sex, Condition) %>% 
+  tally() %>% arrange(Sex) %>% knitr::kable("latex")
+  
+demo_lag <- lag4 %>% right_join(collars_demographics) %>% filter(!is.na(Sex))
+
+Cairo::CairoPNG(filename = "sexTOD.png", width = 1200, height = 1000,
+                res = 180)
+ggplot(demo_lag, aes(x = dt4, fill = ToD)) +
+  geom_density(na.rm = T, alpha = .4) + facet_wrap(~Condition +Sex) +
+  theme_minimal() + 
+  labs(x = "24-hour displacement") +
+  theme(strip.background = element_rect(fill = "#F4F4F4"),
+        plot.title = element_text(hjust = .5))
+dev.off()
+
+ggplot(filter(demo_lag, str_detect(ToD, "^d")), aes(x = dt4, fill = ToD)) +
+  geom_density(na.rm = T, alpha = .4) + facet_wrap(~Condition +Sex)
+
+
+
+##### Regional Differences
+
+# Kunene
+
+# Etosha --- across west/east/central 
+
+# Waterberg --  only 4 individuals. 
+library(stmove)
+means <- dist_map(st_set_geometry(all_rhinos, NULL), proj = "+proj=utm +zone=33 +south +datum=WGS84 +units=m +no_defs")
+park <- read_sf("../../../Dropbox/Research - PHD/Anthrax/ENP data for transmission study/ENP shapefiles/ENP_three_areas.shp") %>% 
+  st_transform(32733) %>% mutate(id = c("East", "West", "Central"))
+
+groups <- st_intersects(means, park, sparse = F) %>% 
+  as_tibble(.name_repair = "unique") %>% 
+  rename(East = "...1", West = "...2", Central = "...3") %>%
+  mutate(ids = means$id, 
+         Area = case_when(
+           East == TRUE ~ "East",
+           West == TRUE ~ "West", 
+           Central == TRUE  ~ "Central", 
+           TRUE ~ "OOP"
+         ))
+
+waterberg_ids <- paste0("SAT", c(469, 676:678, 680:683))
+groups[groups$ids %in% waterberg_ids, "Area"] <- "Waterberg"
+
+groups %>% group_by(Area) %>% tally()
+
+Rhino_meanlocs <- left_join(means, groups, by = c("id" = "ids"))
diff --git a/SupplementaryAnalyses/pairwise_ks_test.R b/SupplementaryAnalyses/pairwise_ks_test.R
new file mode 100644
index 0000000..baa62c9
--- /dev/null
+++ b/SupplementaryAnalyses/pairwise_ks_test.R
@@ -0,0 +1,41 @@
+pairwise_ks_test <- function(value, group, alternative = "two.sided"){
+  
+
+  
+  f <- function(x, y){ 
+    p <- ks.test(x, y, alternative = alternative, exact = 
+                   NULL)$p.value 
+    return(p) 
+  }
+  
+  res <- suppressWarnings(lapply(lst, function(x) lapply(lst, function(y) f(x, y)))) 
+  res <- unlist(res)
+  res <- res[res!=1] %>% unique()
+  res <- p.adjust(res, method = "bonferroni")
+  res <- matrix(res, nrow = length(lst), ncol = length(lst), byrow = T)
+  row.names(res) <- colnames(res) <- names(lst)
+  cat("Pairwise Kolmogorov-Smirnov Test with Correction: p-value Matrix","\n","\n")
+  return(res)
+}
+
+
+#dumbest way i can figure this out tonight... 
+expand.grid(1:4, 1:4) %>% filter(Var1 != Var2) -> x
+x[x[,2]>=x[,1],] 
+
+comb <- expand.grid(unique(levels(df$ToD)), unique(levels(df$ToD))) %>% 
+  filter(Var1 != Var2) %>% .[x[,2]>=x[,1],] %>% mutate_all(as.character())
+
+
+ToD_disp <- filter(lag4, !is.na(dt4)) %>%
+  mutate(disphr = dt4/time4) %>% 
+  select(ToD, disphr) %>% 
+  split(.$ToD, ) %>%
+  map(~pull(.x, disphr))
+
+results <- pmap_dbl(list(a = comb$Var1, b = comb$Var2), 
+     function(a,b){
+  suppressWarnings(ks.test(ToD_disp[[a]], ToD_disp[[b]])$p.value)
+}) %>% p.adjust(method = "bonferroni") 
+
+names(results) <- paste(comb$Var1, comb$Var2, sep =".")