---
title: "history"
output: html_document
date: "2025-03-13"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
library(dplyr)
library(kableExtra)
library(ggplot2)
library(officer)
library(flextable)
library(openxlsx)
library(gridExtra)
library(ggtext)
```

```{r read_and_prepare_data, include=FALSE}
acet <- read.table(file='R_acetosa_observations.tsv', sep = '\t', header = TRUE, na.strings = c("NA", "", " "))
acet <- cbind(species="R. acetosa", acet)
acet$date <- gsub("X2021\\.", "", acet$date)
thys <- read.table(file='R_thyrsiflorus_observations.tsv', sep = '\t', header = TRUE, na.strings = c("NA", "", " "))
thys <- cbind(species="R. thyrsiflorus", thys)
thys$date <- gsub("X2021\\.", "", thys$date)

data <- bind_rows(acet, thys)
data <- data %>% 
  filter(!is.na(date))

ac_true <- c(F=18, M=17)
at_true <- c(F=17, M=9)
true=c("R. acetosa"=ac_true, "R. thyrsiflorus"=at_true)

```

```{r functions, include=FALSE}
save_plots <- function(p, name, spec, additional){
  ggsave(filename = paste(name, spec, "-", additional,".svg", sep = ""), plot = p, width = 10, height = 3.5, dpi = 300, bg = "transparent")
  ggsave(filename = paste(name, spec, "-", additional,".pdf", sep = ""), plot = p, width = 10, height = 3.5, bg = "transparent")
  ggsave(filename = paste(name, spec, "-", additional,".png", sep = ""), plot = p, width = 10, height = 3.5, dpi = 300, bg = "transparent")
}
plot_1 <- function(df, spec, species, prev, tekst, name = "", sex_prop) {
  df_RA <- df %>% 
    filter(species == spec) %>%
    mutate(date = as.Date(date, format = "%m.%d"))
  
  date_min <- min(df_RA$date) - prev
  date_max <- max(df_RA$date) 
  print(paste("******* min", min(df_RA$date), "Date min:", date_min, "\n"))
  

  formatted_caption <- paste("<b>",tekst,"</b>", "<b> </b>","<i> </i>", "<i>", species, "</i>",name, sep=" ")
 
  
  p <- ggplot(df_RA, aes(x = date, y = proportion)) +
    geom_line(color = "gray") +
    geom_hline(yintercept = as.numeric(sex_prop), linetype = "dashed", color = "darkgrey", size = 0.75) +
    geom_point(size = 3, shape = 21, fill = "gray", color = "black", stroke = 1) +
    geom_point(data = df_RA %>% filter(proportion > 0.5, p_value < 0.05),
               aes(x = date, y = proportion), color = "black", size = 3, stroke = 1) +
    geom_point(data = df_RA %>% filter(proportion < 0.5, p_value < 0.05),
               aes(x = date, y = proportion), shape = 21, fill = "white", color = "black", size = 3, stroke = 1) +
    scale_x_date(date_labels = "%m.%d", date_breaks = "1 day", breaks = df_RA$date, minor_breaks = NULL,
                 limits = c(date_min, date_max)) +
    scale_y_continuous(limits = c(-0.05, 1.05), expand = c(0,0)) +
    labs(caption = formatted_caption,
         x = "Date",
         y = "Frequency of males") +
    theme_minimal() +
    theme(
      axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
      plot.caption = element_markdown(
        hjust = 0, 
        size = 12, 
        color = "black", 
        margin = margin(t = 5, b = 0)
      ),
      plot.caption.position = "plot",
      plot.margin = margin(t = 8, r = 5, b = 10, l = 5, unit = "pt")
    )
  save_plots(p, name, spec,  "history")
  
  
  return(p)
}

plot_1_diff <- function(df, spec, prev, tekst,  name="") {

  df_RA <- df %>% 
    filter(species == spec) %>%
    mutate(date = as.Date(date, format = "%m.%d"))
  
  date_min <- min(df_RA$date) - prev
  date_max <- max(df_RA$date) 
  formatted_caption <- paste0(tekst, " <i>", spec, "</i>")
  date_range_extended <- c(min(df_RA$date) - 1, max(df_RA$date) + 1)

  p <- ggplot(df_RA, aes(x = date, y = diff_prop)) +
    geom_line(color = "gray") +
    
    geom_point(size = 3, shape = 21, fill = "gray", color = "black", stroke = 1) +
    

    geom_point(data = df_RA %>% filter(diff_prop > 0, tt_p_val < 0.05),
               aes(x = date, y = diff_prop), color = "black", size = 3, stroke = 1) +
    

    geom_point(data = df_RA %>% filter(diff_prop < 0, tt_p_val < 0.05),
               aes(x = date, y = diff_prop), shape = 21, fill = "white", color = "black", size = 3, stroke = 1) +
    
    scale_x_date(date_labels = "%m.%d", date_breaks = "1 day", breaks = df_RA$date, minor_breaks = NULL,
                 limits = c(date_min, date_max)) +
    scale_y_continuous(limits = c(-0.1, 1.1), expand = c(0,0)) +
    ylim(-1, 1) +
    labs(title = paste(spec, "- difference in observable frequency of males in time"),
         x = "Date",
         y = "difference in proportion") +
    theme_minimal() +
    theme(
      axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),

      plot.caption = element_markdown(
        hjust = 0, 
        size = 10, 
        color = "black", 
        margin = margin(t = -10, b = 0)
      ),
      plot.caption.position = "plot",
      plot.margin = margin(t = 5, r = 5, b = 40, l = 5, unit = "pt")
    )
  save_plots(p, name, spec,  "history")

  
  return(p)
}

```

```{r binom_test, include=FALSE}
df <- data
df <- df %>%
  rowwise() %>%
  mutate(
    # proportions (male / (male + female))
    proportion = male / (male + female),
    true_prop = true[paste(species,".M",sep='')]/
                        (true[paste(species,".F",sep='')] + 
                                true[paste(species,".M",sep='')]),
    # binomial test H0: p = 0.5
    p_value = round(binom.test(x = male, n = male + female, p = 0.5)$p.value, 3),
    
    # conf. intervals (default 95%)
    conf_low = round(binom.test(x = male, n = male + female, p = 0.5)$conf.int[1], 3),
    conf_high = round(binom.test(x = male, n = male + female, p = 0.5)$conf.int[2], 3),
    diff_prop = round((male / (male + female)) -
                       (true[paste(species,".M",sep='')]/
                        (true[paste(species,".F",sep='')] + 
                                true[paste(species,".M",sep='')]) ), 2),
    tt_p_val = round(prop.test(x=c(male, true[paste(species,".M",sep='')]), 
                         n=c(male + female, (true[paste(species,".F",sep='')] + 
                                true[paste(species,".M",sep='')])))$p.value, 4),
    dp = round(proportion-true_prop, 2)
    
  ) %>%
  ungroup()

write.table(df, file = "history.tsv", 
            sep = "\t",           
            row.names = FALSE,    
            quote = FALSE,        
            na = "NA")

```

```{r plots, include=FALSE}

p1 <- plot_1(df, "R. acetosa", "R. acetosa",1, "a)","", 0.49)
p2 <- plot_1(df, "R. thyrsiflorus", "R. thyrsiflorus",8,"b)", "", 0.35)
plot_1_diff(df, "R. acetosa", 1,"a)" ,"diff_")
plot_1_diff(df, "R. thyrsiflorus", 8,"b)", "diff_")
```

```{r oba}
df <- data

df_summary <- df %>%
  group_by(date) %>%
  summarise(
    male = sum(male),
    female = sum(female)
  ) %>%
  mutate(species = "both") %>%
  select(species, date, male, female)


df_summary <- df_summary %>%
  mutate(date = format(as.Date(date, format = "%m.%d"), "%m.%d")) %>%
  arrange(date)


df_summary <- df_summary %>%
  rowwise() %>%
  mutate(

    proportion = round(male / (male + female), 2),

    p_value = round(binom.test(x = male, n = male + female, p = 0.5)$p.value, 3),

    conf_low = round(binom.test(x = male, n = male + female, p = 0.5)$conf.int[1], 3),
    conf_high = round(binom.test(x = male, n = male + female, p = 0.5)$conf.int[2], 3)
  ) %>%
  ungroup()

write.table(df_summary, file = "history-both.tsv", 
            sep = "\t",           
            row.names = FALSE,    
            quote = FALSE,        
            na = "NA")


p3 <- plot_1(df_summary, "both", "R. acetosa + R. thyrsiflorus", 1, "c)","", 0.42)

p_hist <- grid.arrange(p1, p2, p3, ncol = 1, heights = c(16, 16, 16))
#p_hist
name="all_history"
ggsave(filename = paste(name, ".svg", sep = ""), plot = p_hist, width = 10, height = 8, dpi = 300, bg = "transparent")
ggsave(filename = paste(name, ".pdf", sep = ""), plot = p_hist, width = 10, height = 8, bg = "transparent")
ggsave(filename = paste(name, ".png", sep = ""), plot = p_hist, width = 10, height = 8, dpi = 300, bg = "transparent")
```