From babd747f583137faa7d57ff9209d69a54a882a5f Mon Sep 17 00:00:00 2001 From: Dana Seidel Date: Thu, 18 Jul 2019 16:58:29 -0700 Subject: [PATCH 1/9] typo fix in Figure 8 --- Figures_Script.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Figures_Script.R b/Figures_Script.R index 032e5fc..0eb928d 100644 --- a/Figures_Script.R +++ b/Figures_Script.R @@ -3,7 +3,7 @@ ### (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 @@ -151,5 +151,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)) From 59df261887b79c0ee2a967edbf14490637776758 Mon Sep 17 00:00:00 2001 From: Dana Seidel Date: Thu, 18 Jul 2019 16:58:54 -0700 Subject: [PATCH 2/9] Outline major revisions --- README.md | 5 +++-- SupplementaryAnalyses/BiweekelyNullModel.R | 6 ++++++ SupplementaryAnalyses/HomerangeMethods.R | 1 + SupplementaryAnalyses/SeasonalGroupEffects.R | 2 ++ 4 files changed, 12 insertions(+), 2 deletions(-) create mode 100644 SupplementaryAnalyses/BiweekelyNullModel.R create mode 100644 SupplementaryAnalyses/HomerangeMethods.R create mode 100644 SupplementaryAnalyses/SeasonalGroupEffects.R 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..45d46cb --- /dev/null +++ b/SupplementaryAnalyses/BiweekelyNullModel.R @@ -0,0 +1,6 @@ +##### 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? diff --git a/SupplementaryAnalyses/HomerangeMethods.R b/SupplementaryAnalyses/HomerangeMethods.R new file mode 100644 index 0000000..a44e4f6 --- /dev/null +++ b/SupplementaryAnalyses/HomerangeMethods.R @@ -0,0 +1 @@ +##### LoCoh v. AKDE ##### diff --git a/SupplementaryAnalyses/SeasonalGroupEffects.R b/SupplementaryAnalyses/SeasonalGroupEffects.R new file mode 100644 index 0000000..4d87a1f --- /dev/null +++ b/SupplementaryAnalyses/SeasonalGroupEffects.R @@ -0,0 +1,2 @@ +##### Group effects, Seasonality + From 44522628d3a40ffcf18ba101e88df167be0c9457 Mon Sep 17 00:00:00 2001 From: Dana Seidel Date: Tue, 23 Jul 2019 10:24:59 -0700 Subject: [PATCH 3/9] Adaption to Time of Day filter to correctly select closest mid-night fixes from empirical trajectories. Previous code was incorrectly calculating offset and grouping days. --- 24hrDisplacement_Script.R | 33 +++++++++++++++++++++++---------- 1 file changed, 23 insertions(+), 10 deletions(-) diff --git a/24hrDisplacement_Script.R b/24hrDisplacement_Script.R index 561a5e6..f26171c 100644 --- a/24hrDisplacement_Script.R +++ b/24hrDisplacement_Script.R @@ -7,7 +7,11 @@ source("DataParsing_Script.R") library(caret) ## Time of Day displacement analysis -ToDfilter <- all_rhinos %>% st_set_geometry(NULL) %>% +sampling_days <- c(unique(all_rhinos$local_date), ymd("2017-09-20")) +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,13 +27,25 @@ 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"))) @@ -41,9 +57,6 @@ ToD <- ToDfilter %>% split(.$id) %>% arrange(local_date, ToD)) lag4 <- ToD %>% - keep(., function(x) { - na.omit(unique(x$id)) != "SAT645" - }) %>% # drop 645, too short, throws error. map_df(., function(x) { mutate(x, dx = c(diff(x), NA), From 6374931ab908daa4e0cbcfb4a44cd70ff5c879c9 Mon Sep 17 00:00:00 2001 From: Dana Seidel Date: Tue, 23 Jul 2019 19:05:44 -0700 Subject: [PATCH 4/9] Calculate empirical time between points (reviewer 2 suggestion). Further adaptions to code to correctly ensure complete intervals and correct ordering before diff calculation. --- 24hrDisplacement_Script.R | 53 +++++++++++++++++++++++++++++++-------- 1 file changed, 43 insertions(+), 10 deletions(-) diff --git a/24hrDisplacement_Script.R b/24hrDisplacement_Script.R index f26171c..61be6f4 100644 --- a/24hrDisplacement_Script.R +++ b/24hrDisplacement_Script.R @@ -7,7 +7,8 @@ source("DataParsing_Script.R") library(caret) ## Time of Day displacement analysis -sampling_days <- c(unique(all_rhinos$local_date), ymd("2017-09-20")) +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"))) @@ -51,14 +52,19 @@ ToDfilter <- all_rhinos %>% st_set_geometry(NULL) %>% 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 %>% 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), @@ -67,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")) @@ -103,10 +136,10 @@ other_class <- filter(classes, Class == "Other") %>% pull(x) ks.test(dawn_class, other_class) ## Time of day points near water? -buffers <- st_buffer(water_utm, 500) +buffers <- st_buffer(water_utm, 250) 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_pts %>% group_by(ToD) %>% tally() From 05a76a1c1c127f5f74869de792eaf183cb087aa7 Mon Sep 17 00:00:00 2001 From: Dana Seidel Date: Tue, 23 Jul 2019 19:06:16 -0700 Subject: [PATCH 5/9] Add data to repository and gitignore, correct paths, to make revision testing easier. --- .gitignore | 1 + 16dayRecursion_HomeRange_Script.R | 4 ++-- DataParsing_Script.R | 6 +++--- 3 files changed, 6 insertions(+), 5 deletions(-) diff --git a/.gitignore b/.gitignore index d714431..4e6b125 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,4 @@ .Rhistory *.Rproj .DS_Store +Data/ diff --git a/16dayRecursion_HomeRange_Script.R b/16dayRecursion_HomeRange_Script.R index 248c631..0369a45 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) @@ -226,7 +226,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) 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 From 109728fc6c52bf63d353848c0af5217fa4167b89 Mon Sep 17 00:00:00 2001 From: Dana Seidel Date: Tue, 23 Jul 2019 22:53:22 -0700 Subject: [PATCH 6/9] minor label fix Fig 3. --- Figures_Script.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Figures_Script.R b/Figures_Script.R index 0eb928d..b9a5a57 100644 --- a/Figures_Script.R +++ b/Figures_Script.R @@ -67,7 +67,7 @@ chosen9 %>% theme_minimal() + labs( title = "Dawn-Dawn 24-displacement through time", - x = "Julian day", y = " 24-hr displacement (m)", + 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" ) + theme(plot.title = element_text(hjust = .5), legend.position = "none", panel.spacing = unit(1, "lines")) From 0447cc487ea7aa7eee55bb84cb6d0e0af4697c99 Mon Sep 17 00:00:00 2001 From: Dana Seidel Date: Tue, 23 Jul 2019 22:53:41 -0700 Subject: [PATCH 7/9] Add periodicity tests -- reviewer 2 suggestion --- 24hrDisplacement_Script.R | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/24hrDisplacement_Script.R b/24hrDisplacement_Script.R index 61be6f4..9f02acb 100644 --- a/24hrDisplacement_Script.R +++ b/24hrDisplacement_Script.R @@ -127,6 +127,15 @@ chosen9 <- con_12hr[ids] %>% group_by(id) %>% mutate(index = seq(1, n(), 1)) +### 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$dt2))}) +# 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)) From b80675d74a137d0ff542cb176402feca443b29c9 Mon Sep 17 00:00:00 2001 From: Dana Seidel Date: Tue, 13 Aug 2019 20:09:32 -0700 Subject: [PATCH 8/9] revisions --- .gitignore | 2 + 16dayRecursion_HomeRange_Script.R | 20 +-- 24hrDisplacement_Script.R | 83 +++++++--- AnnualRecursion_Script.R | 16 +- Figures_Script.R | 35 ++-- SupplementaryAnalyses/BiweekelyNullModel.R | 151 ++++++++++++++++++ SupplementaryAnalyses/DisplacementPerHour.R | 12 ++ SupplementaryAnalyses/HomerangeMethods.R | 137 ++++++++++++++++ SupplementaryAnalyses/SeasonalGroupEffects.R | 159 +++++++++++++++++++ SupplementaryAnalyses/pairwise_ks_test.R | 41 +++++ 10 files changed, 606 insertions(+), 50 deletions(-) create mode 100644 SupplementaryAnalyses/DisplacementPerHour.R create mode 100644 SupplementaryAnalyses/pairwise_ks_test.R diff --git a/.gitignore b/.gitignore index 4e6b125..1d73627 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,5 @@ *.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 0369a45..da29a35 100644 --- a/16dayRecursion_HomeRange_Script.R +++ b/16dayRecursion_HomeRange_Script.R @@ -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 ######## @@ -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 9f02acb..1f1a4fd 100644 --- a/24hrDisplacement_Script.R +++ b/24hrDisplacement_Script.R @@ -115,17 +115,30 @@ 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) %>% @@ -133,22 +146,54 @@ chosen9 %>% split(.$id) %>% # a couple have 105-155 relocations but most < 100 chosen9 %>% split(.$id) %>% - map(function(x){ptest::ptestg(na.omit(x$dt2))}) + 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],] -## Time of day points near water? -buffers <- st_buffer(water_utm, 250) -water_pts <- st_intersection(buffers, - st_as_sf(lag4, - coords = c("x", "y"), - crs = 32733, na.fail = F)) -water_pts %>% group_by(ToD) %>% tally() +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? +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/Figures_Script.R b/Figures_Script.R index b9a5a57..793cdb7 100644 --- a/Figures_Script.R +++ b/Figures_Script.R @@ -42,6 +42,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 +60,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-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 +95,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 stardardized 24-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 +150,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) %>% diff --git a/SupplementaryAnalyses/BiweekelyNullModel.R b/SupplementaryAnalyses/BiweekelyNullModel.R index 45d46cb..cd92bd8 100644 --- a/SupplementaryAnalyses/BiweekelyNullModel.R +++ b/SupplementaryAnalyses/BiweekelyNullModel.R @@ -4,3 +4,154 @@ ### 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(meanNDVI, 2), data = mod_biweekly, family = "poisson") + + +anova(mod1, mod2, mod3) +AIC(mod1, mod2, mod3) + +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() 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 index a44e4f6..5254c8c 100644 --- a/SupplementaryAnalyses/HomerangeMethods.R +++ b/SupplementaryAnalyses/HomerangeMethods.R @@ -1 +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 index 4d87a1f..f5d1d56 100644 --- a/SupplementaryAnalyses/SeasonalGroupEffects.R +++ b/SupplementaryAnalyses/SeasonalGroupEffects.R @@ -1,2 +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 =".") From 90588c69717d07aa5fc642c3ed1f76b9496b2f45 Mon Sep 17 00:00:00 2001 From: Dana Seidel Date: Sun, 27 Oct 2019 13:45:27 -0700 Subject: [PATCH 9/9] final revisions --- Figures_Script.R | 5 +- SupplementaryAnalyses/BiweekelyNullModel.R | 62 ++++++++++++++++++++-- 2 files changed, 60 insertions(+), 7 deletions(-) diff --git a/Figures_Script.R b/Figures_Script.R index 793cdb7..a86cfb5 100644 --- a/Figures_Script.R +++ b/Figures_Script.R @@ -8,7 +8,6 @@ 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 @@ -68,7 +67,7 @@ chosen9 %>% facet_wrap(~id_f, scales = "free_x", ncol = 3) + theme_minimal() + labs( - title = "Dawn-Dawn standardized 24-displacement through time", + 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" ) + @@ -104,7 +103,7 @@ chosen9 %>% facet_wrap(~id_f, ncol = 3) + theme_minimal() + labs( - title = "Dawn-Dawn stardardized 24-displacement through time", + 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"), diff --git a/SupplementaryAnalyses/BiweekelyNullModel.R b/SupplementaryAnalyses/BiweekelyNullModel.R index cd92bd8..c1898fc 100644 --- a/SupplementaryAnalyses/BiweekelyNullModel.R +++ b/SupplementaryAnalyses/BiweekelyNullModel.R @@ -135,11 +135,11 @@ mod_biweekly <- biweekly_null %>% 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(meanNDVI, 2), 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) -AIC(mod1, mod2, mod3) +anova(mod1, mod2, mod3, mod4) +AIC(mod1, mod2, mod3, mod4) Cairo::CairoPNG(filename = "quad.png", width = 1200, height = 1000, res = 180) @@ -155,3 +155,57 @@ mod_biweekly %>% 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()