paceR: Tools for working with the PACE Database

library(paceR)
Sys.setenv(TZ = 'UTC')

source("~/Dropbox/R/themes_fc.R")
source("~/Dropbox/R/helpers_fc.R")
source("~/Dropbox/R/plotting-functions.R")
library(tidyverse)
library(forcats)
library(lubridate)
library(RColorBrewer)

pace_db <- DBI::dbConnect(RMySQL::MySQL(), group = "PACE", user = "camposf", dbname = "monkey")
paceR_db <- DBI::dbConnect(RMySQL::MySQL(), group = "PACE", user = "camposf", dbname = "paceR")

census <- get_pace_tbl(paceR_db, "vCensusAnnual")

census <- census %>%
  mutate_each(funs(as.Date), DateStart, DateEnd, DateOf)

census <- filter(census, Project == "Santa Rosa")

census_group <- get_pace_tbl(pace_db, "tblGroup") %>%
  select(GroupID = ID, PrimateSpeciesID, NameLong)

species <- get_pace_tbl(pace_db, "tblPrimateSpecies") %>%
  select(PrimateSpeciesID = ID, Species = NameShort)

census <- census %>%
  inner_join(census_group) %>%
  inner_join(species)

`%ni%` = Negate(`%in%`)
bad_years <- c(1989, 1991, 1993:1998)

t1 <- census %>%
  filter(UseForTotalEstimate == 1) %>%
  group_by(Species, YearOfCensus) %>%
  summarise(total_n = sum(N)) %>%
  filter(total_n > 100)

# Add additional points manually
m <- tibble(Species = c("Mantled Howlers", "Capuchins"),
            YearOfCensus = c(1972, 1972),
            total_n = c(85, 297))
t1 <- bind_rows(t1, m)

# h_color <- brewer.pal(3, "Set1")[1]
h_color <- "firebrick"
c_color <- brewer.pal(3, "Set1")[2]

# Total population (howlers)
ggplot(filter(t1, Species == "Mantled Howlers" & YearOfCensus >= 1983),
       aes(x = YearOfCensus, y = total_n)) +
  # annotate(geom = "rect", xmin = 2007, xmax = 2013, ymin = -Inf, ymax = Inf,
  #          fill = "gray90") +
  geom_line(color = h_color) +
  geom_point(shape = 21, color = "white", fill = h_color, size = 3) +
  annotate(geom = "text", x = 1990, y = 300, label = "Mantled Howlers",
           hjust = 0, color = h_color) +
  coord_cartesian(y = c(0, 700)) +
  labs(x = "Year", y = "Number of Animals",
       title = "Total Population Size") +
  scale_x_continuous(breaks = seq(1980, 2015, by = 5)) +
  theme_journal_x2() +
  scale_color_manual(values = h_color,
                     name = "", guide = FALSE)

# Total population (both)
ggplot(filter(t1, YearOfCensus >= 1983),
       aes(x = YearOfCensus, y = total_n, color = Species, fill = Species)) +
  # annotate(geom = "rect", xmin = 2007, xmax = 2013, ymin = -Inf, ymax = Inf,
  #          fill = "gray90") +
  geom_line() +
  geom_point(shape = 21, color = "white", size = 3) +
  annotate(geom = "text", x = 1990, y = 300, label = "Mantled Howlers",
           hjust = 0, color = h_color) +
  annotate(geom = "text", x = 1983, y = 550, label = "Capuchins",
           hjust = 0, color = c_color) +
  coord_cartesian(y = c(0, 700)) +
  labs(x = "Year", y = "Number of Animals",
       title = "Total Population Size") +
  scale_x_continuous(breaks = seq(1980, 2015, by = 5)) +
  theme_journal_x2() +
  scale_color_manual(values = c(c_color, h_color),
                     name = "", guide = FALSE) +
  scale_fill_manual(values = c(c_color, h_color),
                    name = "", guide = FALSE) +
  geom_smooth(method = "loess", color = NA, span = 1)



t1 %>%
  spread(Species, total_n) %>%
  View()

census %>%
  filter(UseForTotalEstimate == 1) %>%
  group_by(Species, YearOfCensus, GroupID) %>%
  summarise(group_n = sum(N)) %>%
  ungroup() %>%
  group_by(Species, YearOfCensus) %>%
  summarise(total_n = sum(group_n),
            smallest = min(group_n),
            largest = max(group_n),
            mean = mean(group_n, na.rm = TRUE),
            n_groups = n()) %>%
  View()

d_alo <- census %>%
  filter(UseForDemography == 1 & Species == "Mantled Howlers") %>%
  unite(AgeSexClass, AgeClass, Sex)

d_alo$AgeSexClass <- fct_collapse(d_alo$AgeSexClass,
                                  "AdultFemale" = "A_F",
                                  "AdultMale" = "A_M",
                                  "Juvenile" = c("J_U", "LIM_U", "SA_M", "SIM_U"),
                                  "Infant" = c("I_U", "ID_U", "II_U"),
                                  "Unknown" = "U_U")
d_alo_summary <- d_alo %>%
  filter(YearOfCensus %ni% bad_years) %>%
  group_by(YearOfCensus, AgeSexClass) %>%
  summarise(n = sum(N)) %>%
  spread(AgeSexClass, n) %>%
  mutate(af_inf = AdultFemale / Infant,
         af_imm = AdultFemale / (Infant + Juvenile))



d_ceb <- census %>%
  filter(UseForDemography == 1 & Species == "Capuchins") %>%
  unite(AgeSexClass, AgeClass, Sex)

d_ceb$AgeSexClass <- fct_collapse(d_ceb$AgeSexClass,
                                  "AdultFemale" = "A_F",
                                  "AdultMale" = "A_M",
                                  "SubAdultMale" = "SA_M",
                                  "LImm" = c("LIM_U", "LIM_F", "LIM_M"),
                                  "SImm" = c("SIM_U", "SIM_F", "SIM_M"),
                                  "Juv" = c("J_U"),
                                  "Infant" = c("I_U", "ID_U", "II_U", "I_F", "I_M"),
                                  "Unknown" = "U_U",
                                  "Unborn" = "UB_U")
d_ceb_summary <- d_ceb %>%
  filter(YearOfCensus %ni% bad_years) %>%
  group_by(YearOfCensus, AgeSexClass) %>%
  summarise(n = sum(N))

d_ceb_summary %>%
  ungroup() %>%
  group_by(AgeSexClass) %>%
  mutate(diff = n - lag(n),
         prop_diff = 100 * (diff / lag(n))) %>%
  select(-diff, -prop_diff) %>%
  spread(AgeSexClass, n) %>%
  View()

d_ceb_plot <- filter(d_ceb_summary, AgeSexClass %ni% c("Unknown", "UB_U", "Juv"))

d_ceb_plot$AgeSexClass <- fct_collapse(d_ceb_plot$AgeSexClass,
                                     "Adult Male" = c("AdultMale",
                                                      "SubAdultMale"))

d_ceb_plot$AgeSexClass <- fct_recode(d_ceb_plot$AgeSexClass,
                                     "Adult Female" = "AdultFemale",
                                     "Large Immature" = "LImm",
                                     "Small Immature" = "SImm")

d_ceb_plot <- d_ceb_plot %>%
  ungroup() %>%
  group_by(YearOfCensus, AgeSexClass) %>%
  summarise(n = sum(n)) %>%
  ungroup() %>%
  group_by(YearOfCensus) %>%
  mutate(prop = n / sum(n, na.rm = TRUE))

d_ceb_plot$AgeSexClass <- factor(d_ceb_plot$AgeSexClass,
                                 levels = c("Adult Female",
                                            "Adult Male",
                                            "Large Immature",
                                            "Small Immature",
                                            "Infant"))

census_years <- as.numeric(levels(factor(d_ceb_plot$YearOfCensus)))

ggplot(d_ceb_plot, aes(x = YearOfCensus, y = prop)) +
  geom_line(size = 1) +
  geom_point(shape = 21, fill = "black", color = "white", size = 2,
             stroke = 0.75) +
  facet_grid(AgeSexClass ~ .) +
  theme_journal() +
  scale_x_continuous(breaks = census_years) +
  labs(x = "Year", y = "Proportion of Population") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplot(d_ceb_plot, aes(x = YearOfCensus, y = prop, fill = AgeSexClass)) +
  geom_area(stat = "identity", position = "stack") +
  theme_journal() +
  scale_x_continuous(breaks = census_years) +
  labs(x = "Year", y = "Proportion of Population") +
  scale_fill_brewer(palette = "Set1", name = NULL) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

d_ceb_grp <- d_ceb %>%
  filter(YearOfCensus %ni% bad_years) %>%
  group_by(YearOfCensus) %>%
  mutate(n_grps = n_distinct(GroupID))

d_ceb_grp <- filter(d_ceb_grp,
                    AgeSexClass %ni% c("Unknown", "UB_U", "Juv"))

d_ceb_grp$AgeSexClass <- fct_collapse(d_ceb_grp$AgeSexClass,
                                      "Adult Male" = c("AdultMale",
                                                       "SubAdultMale"))

d_ceb_grp$AgeSexClass <- fct_recode(d_ceb_grp$AgeSexClass,
                                    "Adult Female" = "AdultFemale",
                                    "Large Immature" = "LImm",
                                    "Small Immature" = "SImm")

d_ceb_grp <- d_ceb_grp %>%
  group_by(YearOfCensus, AgeSexClass) %>%
  summarise(n_indivs = sum(N),
            n_grps = first(n_grps)) %>%
  ungroup() %>%
  mutate(n_per_grp = n_indivs / n_grps)

d_ceb_grp$AgeSexClass <- factor(d_ceb_grp$AgeSexClass,
                                levels = c("Adult Female",
                                           "Adult Male",
                                           "Subadult Male",
                                           "Large Immature",
                                           "Small Immature",
                                           "Infant"))

ggplot(d_ceb_grp, aes(x = YearOfCensus, y = n_per_grp)) +
  geom_line(size = 1) +
  geom_point(shape = 21, fill = "black", color = "white", size = 2,
             stroke = 0.75) +
  facet_grid(AgeSexClass ~ .) +
  theme_journal() +
  scale_x_continuous(breaks = census_years) +
  labs(x = "Year", y = "Individuals Per Group") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_y_continuous(breaks = 0:6) +
  expand_limits(y = 0)


outl <- filter(d_alo_summary, YearOfCensus == 2011)

# AF-Infant ratio
ggplot(d_alo_summary, aes(x = YearOfCensus, y = af_inf)) +
  annotate(geom = "rect", xmin = 2007, xmax = 2013, ymin = -Inf, ymax = Inf,
           fill = "gray90") +
  geom_line() +
  annotate(geom = "point", x = outl$YearOfCensus, y = outl$af_inf,
           size = 4, color = h_color) +
  geom_point(shape = 21, fill = "black", color = "white", size = 2) +
  geom_smooth(method = "loess", color = h_color, fill = h_color) +
  expand_limits(y = 0.5) +
  geom_hline(yintercept = 1, lty = 2, color = "gray50") +
  theme_journal_x2() +
  labs(x = "Year", y = "# Adult Females / # Infants",
       title = "Adult Female-to-Infant Ratio") +
  scale_x_continuous(breaks = seq(1980, 2015, by = 5))

g_alo <- d_alo %>%
  filter(YearOfCensus %ni% bad_years) %>%
  group_by(YearOfCensus, GroupID) %>%
  summarise(group_n = sum(N)) %>%
  summarise(group_mean = mean(group_n),
            group_med = median(group_n),
            group_max = max(group_n),
            group_min = min(group_n))

temp <- gather(g_alo, var, value, -YearOfCensus)

# Group Size
ggplot(filter(g_alo, !YearOfCensus %in% c(1996, 1997, 1998)),
       aes(x = YearOfCensus)) +
  annotate(geom = "rect", xmin = 2007, xmax = 2013, ymin = -Inf, ymax = Inf,
           fill = "gray90") +
  geom_line(aes(y = group_med), color = h_color) +
  geom_line(aes(y = group_max), color = h_color, alpha = 0.7, size = 0.3) +
  geom_line(aes(y = group_min), color = h_color, alpha = 0.7, size = 0.3) +
  geom_point(aes(y = group_med), shape = 21, size = 3,
             color = "white", fill = h_color) +
  geom_point(aes(y = group_max), shape = 21, size = 1.5,
             color = "white", fill = h_color) +
  geom_point(aes(y = group_min), shape = 21, size = 1.5,
             color = "white", fill = h_color) +
  geom_ribbon(aes(ymin = group_min, ymax = group_max),
              fill = h_color, alpha = 0.3) +
  annotate(geom = "text", x = 1995, y = 43, label = "Max", color = h_color) +
  annotate(geom = "text", x = 1995, y = 0, label = "Min", color = h_color) +
  annotate(geom = "text", x = 1995, y = 17, label = "Median", color = h_color) +
  theme_journal_x2() +
  labs(x = "Year", y = "# Animals",
       title = "Group Size") +
  scale_x_continuous(breaks = seq(1980, 2015, by = 5))




temp <- d_alo %>%
  group_by(YearOfCensus) %>%
  summarise(total_n = sum(N))

d_alo$AgeSexClass <- factor(d_alo$AgeSexClass,
                            levels = c("AdultFemale", "AdultMale", "Juvenile", "Infant"))

d_alo_as <- d_alo %>%
  filter(YearOfCensus %ni% bad_years & AgeSexClass != "Unknown") %>%
  group_by(YearOfCensus, AgeSexClass) %>%
  summarise(n = sum(N)) %>%
  inner_join(temp) %>%
  mutate(as_prop = n / total_n)

d_alo_as_med <- d_alo_as %>%
  ungroup() %>%
  group_by(AgeSexClass) %>%
  summarise(class_med = median(as_prop))

ggplot() +
  annotate(geom = "rect", xmin = 2007, xmax = 2013, ymin = -Inf, ymax = Inf,
           fill = "gray90") +
  geom_hline(data = d_alo_as_med, aes(yintercept = class_med),
             color = "gray80") +
  geom_line(data = d_alo_as,
            aes(x = YearOfCensus, y = as_prop),
            color = h_color) +
  geom_point(data = d_alo_as,
             aes(x = YearOfCensus, y = as_prop),
             shape = 21, fill = h_color, color = "white", size = 3) +
  facet_wrap(~AgeSexClass, ncol = 1) +
  scale_x_continuous(breaks = seq(1980, 2015, by = 5)) +
  theme_journal_x2() +
  labs(x = "Year", y = "Proportion of Population",
       title = "Population Age/Sex Composition")

# Diseases
dis <- read_csv("~/Desktop/AAPA/diseases.csv")

ggplot(dis, aes(x = YearOf, y = Cases)) +
  annotate(geom = "rect", xmin = 2007, xmax = 2013, ymin = -Inf, ymax = Inf,
           fill = "gray90") +
  annotate(geom = "text", x = 2000, y = 18000, label = "Dengue",
           color = brewer.pal(3, "Set1")[2]) +
  annotate(geom = "text", x = 2011, y = 3000, label = "Chikungunya",
           color = brewer.pal(3, "Set1")[1]) +
  annotate(geom = "text", x = 2016, y = 500, label = "Zika",
           color = brewer.pal(3, "Set1")[3]) +
  geom_line(aes(color = Disease)) +
  geom_point(aes(fill = Disease),
             shape = 21, color = "white", size = 3) +
  theme_journal_x2() +
  scale_y_continuous(trans = sqrt_sign_trans()) +
  expand_limits(y = -1) +
  scale_color_brewer(palette = "Set1", guide = FALSE) +
  scale_fill_brewer(palette = "Set1", guide = FALSE) +
  theme(plot.caption = element_text(size = 7)) +
  labs(x = "Year", y = "# Cases",
       title = substitute('Cases of'~italic("Flavivirus")~'diseases'),
       subtitle = "Costa Rica",
       caption = "Data: PAHO / WHO")



# Rain
weather <- get_pace_tbl(pace_db, "tblWeather", collect = FALSE) %>%
  filter(ProjectID == 1) %>%
  select(DateOf, TemperatureMax, TemperatureMin, Rainfall) %>%
  collect()

rain <- weather %>%
  mutate(DateOf = ymd(DateOf),
         MonthOf = month(DateOf),
         YearOf = year(DateOf),
         DayOfYear = yday(DateOf)) %>%
  filter(!is.na(DateOf)) %>%
  select(-contains("Temperature"))

rain$MonthOf <- factor(rain$MonthOf, labels = month.abb)

rain_wide <- rain %>%
  group_by(YearOf, MonthOf) %>%
  summarise(Rainfall = sum(Rainfall)) %>%
  spread(YearOf, Rainfall) %>%
  as.data.frame()

rownames(rain_wide) <- rain_wide$MonthOf
rain_wide <- select(rain_wide, -MonthOf)

rain_monthly <- rain %>%
  group_by(YearOf, MonthOf) %>%
  summarise(RainMonthly = sum(Rainfall)) %>%
  ungroup() %>%
  group_by(MonthOf) %>%
  summarise(AvgRainMonthly = mean(RainMonthly, na.rm = TRUE))

rain_yearly <- rain %>%
  group_by(YearOf) %>%
  summarise(RainYearly = sum(Rainfall, na.rm = TRUE))

my_years <- c(rep("Growth", times = length(which(rain_yearly$YearOf <= 1992))),
              rep("Stable", times = length(which(rain_yearly$YearOf > 1992 & rain_yearly$YearOf < 2007))),
              rep("Decline", times = length(which(rain_yearly$YearOf >= 2007 & rain_yearly$YearOf < 2015))),
              rep("Recovery", times = length(which(rain_yearly$YearOf >= 2015))))

my_years <- factor(my_years, levels = c("Growth", "Stable", "Decline", "Recovery"))

library(superheat)


# ---- rainfall_heatmap ---------------------------------------------------

png("~/Desktop/AAPA/superheat2.png", width = 1024 * 1.5, height = 768 * 1.5, res = 150)
superheat(rain_wide,
          scale = F,

          # Aesthetics
          title = "Monthly and Annual Rainfall (mm)",
          title.size = 8,
          heat.pal = brewer.pal(9, "Blues"),
          legend.breaks = seq(0, 1000, by = 250),
          left.label.col = "white",
          bottom.label.col = "white",
          left.label.text.size = 4,
          bottom.label.text.size = 4,

          # bottom.label.text.angle = 90,
          # grid.hline.col = "white",
          # grid.vline.col = "white",
          # grid.vline.size = 0.25,
          # grid.hline.size = 0.25,

          membership.cols = my_years,
          grid.hline.col = "white",
          grid.vline.col = "white",
          grid.vline.size = 2,
          grid.hline.size = 0.25,

          # Right marginal plot
          yr = rain_monthly$AvgRainMonthly,
          yr.axis.name = "Mean Monthly Rainfall",
          yr.plot.type = "bar",
          yr.bar.col = "black",
          yr.obs.col = rep("beige", nrow(rain_monthly)),

          # Top marginal plot
          yt = rain_yearly$RainYearly,
          yt.axis.name = "Total Annual Rainfall",
          yt.plot.type = "bar",
          yt.bar.col = "black",
          yt.obs.col = rep("beige", nrow(rain_yearly)),
          yt.plot.size = 0.6)
dev.off()


# ---- next ---------------------------------------------------------------

rain_monthly2 <- rain %>%
  group_by(YearOf, MonthOf) %>%
  summarise(RainMonthly = sum(Rainfall)) %>%
  inner_join(rain_monthly) %>%
  mutate(RainAnom = RainMonthly - AvgRainMonthly)

rain_monthly2$season <- fct_collapse(rain_monthly2$MonthOf,
                                     EarlyDry = c("Dec", "Jan", "Feb"),
                                     LateDry = c("Mar", "Apr", "May"),
                                     EarlyWet = c("Jun", "Jul", "Aug"),
                                     LateWet = c("Sep", "Oct", "Nov"))

rain_season <- rain_monthly2 %>%
  group_by(YearOf, season) %>%
  summarise(RainSeason = sum(RainMonthly)) %>%
  ungroup() %>%
  group_by(season) %>%
  summarise(AvgRainSeason = mean(RainSeason, na.rm = TRUE))

rain_monthly3 <- rain_monthly2 %>%
  group_by(YearOf, season) %>%
  summarise(RainSeason = sum(RainMonthly)) %>%
  inner_join(rain_season) %>%
  mutate(RainAnom = RainSeason - AvgRainSeason)


lim <- max(abs(rain_monthly2$RainAnom), na.rm = TRUE)

ggplot(filter(rain_monthly2, MonthOf == "Oct"),
       aes(x = YearOf, y = RainAnom)) +
  annotate(geom = "rect", xmin = 2007, xmax = 2013, ymin = -Inf, ymax = Inf,
           fill = "gray90") +
  geom_line() +
  geom_point(aes(fill = RainAnom), shape = 21, color = "black", size = 3) +
  facet_wrap(~MonthOf) +
  geom_hline(yintercept = 0, lty = 2) +
  theme_journal_x2() +
  scale_y_continuous(trans = sqrt_sign_trans()) +
  labs(x = "Year", y = "Rain Anomaly (mm)",
       title = "October Rainfall Anomalies (mm)") +
  scale_x_continuous(breaks = seq(1980, 2015, by = 5)) +
  scale_fill_gradientn(colours = brewer.pal(11, "BrBG"),
                       limits = c(-lim, lim), name = "") +
  theme(legend.key.height = unit(0.3, "cm"),
        legend.key.width = unit(1.5, "cm"))


lim <- max(abs(rain_monthly3$RainAnom), na.rm = TRUE)

ggplot(rain_monthly3, aes(x = YearOf, y = RainAnom)) +
  annotate(geom = "rect", xmin = 2007, xmax = 2013, ymin = -Inf, ymax = Inf,
           fill = "gray90") +
  geom_line() +
  geom_point(aes(fill = RainAnom), shape = 21, color = "black", size = 3) +
  facet_wrap(~season) +
  geom_hline(yintercept = 0, lty = 2) +
  theme_journal_x2() +
  scale_y_continuous(trans = sqrt_sign_trans()) +
  labs(x = "Year", y = "Rain Anomaly (mm)",
       title = "Seasonal Rainfall Anomalies (mm)") +
  scale_x_continuous(breaks = seq(1980, 2015, by = 5)) +
  scale_fill_gradientn(colours = brewer.pal(11, "BrBG"),
                       limits = c(-lim, lim), name = "") +
  theme(legend.key.height = unit(0.3, "cm"),
        legend.key.width = unit(1.5, "cm"))