# Crucial settings
rm(list = ls()) # clean-up the environment

# Load libraries and read and adjust data
library(tidyverse)
library(readxl)
# install.packages("remotes")
# remotes::install_github("jenzopr/silvermantest")
library(silvermantest) # for testing modality
library(circular)

# Download data

td <- read_excel("C:/Users/KWJ/Dropbox/Liak2019/NonBreeders/Time of the day in colony/Files used/Timeday.xlsx")
td <- td %>% mutate(group = paste0(Sex, "_", Breedingstatus, sep = ""))


format_cells <- function(df, rows ,cols, value = c("italics", "bold", "strikethrough")){

  # select the correct markup
  # one * for italics, two ** for bold
  map <- setNames(c("*", "**", "~~"), c("italics", "bold", "strikethrough"))
  markup <- map[value]  

  for (r in rows){
    for(c in cols){

      # Make sure values are not factors
      df[[c]] <- as.character( df[[c]])

      # Update formatting
      df[r, c] <- paste0(markup, df[r, c], markup)
    }
  }

  return(df)
}

Let’s first see overal patterns

# Visualise data - plot ---------------------------------------------------

# With sex
ggplot(data = td, mapping = aes(x = Hours, y = Time, col = Sex)) + 
  geom_point(alpha = 0.1) +
  geom_smooth() +
  facet_wrap(~Breedingstatus, scales = "free") +
  theme_bw() +
  labs(title = "Sex and breeding status separately")

# Pooled sex
ggplot(data = td, mapping = aes(x = Hours, y = Time)) + 
  geom_point(alpha = 0.1) +
  geom_smooth() +
  facet_wrap(~Breedingstatus, scales = "free") +
  theme_bw()   +
  labs(title = "Sexes pooled")

Modality

First, let’s set for modality (with Silverman tests results:) for each sex and breeding status separately

# Testing modality - all groups sex/status separated
gr <- unique(td$group)

# A loop for each status/sex group [i], 
# and k = 1:3 [j], where k is number of modes

slv_k_res <- list()

 
for (j in 1:3) {
  
  slv_res <- numeric()
  
    for (i in 1:length(gr)) {
    
    dt_temp <- td %>% filter(group == gr[i])
    
    slv_res_temp <- silverman.test(dt_temp$Time, k = j)
    slv_res[i] <- slv_res_temp@p_value
    
  }
  slv_k_res[[j]] <- slv_res
  
}

# Some tuning for the final table (content p values for the silverman test)
sil_output <- data.frame(Reduce(rbind, slv_k_res))
colnames(sil_output) <- gr
rownames(sil_output) <- c("k=1", "k=2","k=3")

sil_output %>%
 format_cells(2, 1, "bold") %>%
 knitr::kable()
male_B female_B male_NB female_NB
k=1 0.121121121121121 0.1701702 0.3633634 0.4004004
k=2 0.0800800800800801 0.0970971 0.0850851 0.2782783
k=3 0.158158158158158 0.3383383 0.4954955 0.2672673
There is just a slight tendency for breeding males (male_B) for two peaks of activity (p < 0.1 for k=2)

Now consider modality with sexed pooled

# Testing modality - all sexes combined
gr <- unique(td$Breedingstatus)

# A loop for each status/sex group [i], 
# and k = 1:3 [j], where k is number of modes

slv_k_res <- list()

for (j in 1:3) {
  
  slv_res <- numeric()
  
  for (i in 1:length(gr)) {
    
    dt_temp <- td %>% filter(Breedingstatus == gr[i])
    
    slv_res_temp <- silverman.test(dt_temp$Time, k = j)
    slv_res[i] <- slv_res_temp@p_value
    
  }
  
  slv_k_res[[j]] <- slv_res
  
}

# Some tuning for the final table (content p values for the silverman test)
sil_output <- data.frame(Reduce(rbind, slv_k_res))
colnames(sil_output) <- gr
rownames(sil_output) <- c("k=1", "k=2","k=3")

sil_output %>%
  format_cells(1, 1, "bold") %>%
  format_cells(2, 1, "bold") %>%
  format_cells(3, 2, "bold") %>%
  knitr::kable()
B NB
k=1 0.0810810810810811 0.520520520520521
k=2 0.0600600600600601 0.266266266266266
k=3 0.776776776776777 0.0520520520520521
Looks there is some tendency for one/two and three peaks for breeders and none breeders, respectively, though the stongest results stand for k = 3 for nonbreeders

Distributions

Comparing distribution for status/breeder groups, (with Kolmogorow-Smirnov pairwise tests).
To analyze distribution we first need to unify data (calculate means per hour, to have distribution along a time axis (not time+idvididuals))

First, each group treated separately:

# Adjusting data
td_means <- td %>% group_by(Breedingstatus, Sex, group, Hours) %>% summarise(time_mean = mean(Time)) 

# KS pairwise tests

# remotes::install_github("happyrabbit/DataScienceR")
# library(DataScienceR)
# value <- as.vector(td_means$time_mean)
# group <-as.factor(td_means$group)
# pairwise_ks_test(value,group,warning = -1) # doesn't work for me, and I have not idea why...

# an alternative (to do the same)

gr_means <- unique(td_means$group)
loop_vecA <- c(rep(gr_means[1], 3),
               rep(gr_means[2], 2),
               gr_means[3])

loop_vecB <- c(gr_means[2],
               rep(c(gr_means[3],gr_means[4]), 2),
               gr_means[4])

ks_res <- numeric()

for(i in 1:length(loop_vecA)) {

  dfA <- td_means %>% filter(group == loop_vecA[i])
  dfB <- td_means %>% filter(group == loop_vecB[i])
  
  ks_res[i] <- ks.test(dfA$time_mean, dfB$time_mean)$p.value  

  }

ks_res_df <- data.frame(loop_vecA, loop_vecB, ks_res)
ks_res_df$ks_res <- round(ks_res_df$ks_res, 5)

# sil_output %>%
#   format_cells(1, 1, "bold") %>%
#   format_cells(2, 1, "bold") %>%
#   format_cells(3, 2, "bold") %>%
  knitr::kable(ks_res_df)
loop_vecA loop_vecB ks_res
female_B male_B 0.86743
female_B female_NB 0.00000
female_B male_NB 0.00001
male_B female_NB 0.00000
male_B male_NB 0.00001
female_NB male_NB 0.29213
Well, seems like breeding status but not the sex makes the difference

Let’s then try to run KS test for pooled sexes.

td_means2 <- td %>% group_by(Breedingstatus, Hours) %>% summarise(time_mean = mean(Time)) 

df_N <- td_means %>% filter(Breedingstatus == "B")
df_NB <- td_means %>% filter(Breedingstatus == "NB")

ks.test(df_N$time_mean, df_NB$time_mean)
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  df_N$time_mean and df_NB$time_mean
## D = 0.76984, p-value = 1.312e-11
## alternative hypothesis: two-sided
Looks like the difference between distribution of breeders and nonbreeders is very significant

To plot this:

td_means2 %>% 
  ggplot(aes(x = Hours, y = time_mean)) + 
  geom_point() +
  geom_smooth() +
  facet_wrap(~Breedingstatus, scales = "free")

Note the scales are different for the two panels, and breeders spent just a little time in the colony; what makes a crowd in so call night hours are non-breeders

Circular statistics

# Data convertion:

td_circ <- td %>% 
  mutate(cHr = circular(Hours,
                        type = "angles", 
                        units = "hours",
                        template = "clock24",
                        modulo = "2pi",
                        zero = 0,
                        rotation = "clock")) # conversion of the hour for circular

Testing for uniformity

Plot and formal test with Rayleight test; groups (status/sex) treated separately.

ggplot(td_circ, aes(x = cHr, fill = Sex)) + 
  geom_density(alpha = 0.2) + theme_bw() + 
  coord_polar() + 
  facet_wrap(~ Breedingstatus)

c_gr <- unique(td_circ$group)
ray_p_val <- numeric()

for(i in 1:length(c_gr)) {

  c_df_temp <- td_circ %>% filter(group == c_gr[i])
  ray_p_val[i] <- rayleigh.test(c_df_temp$cHr)$p.value
  
}

ray_res <- data.frame(c_gr, ray_p_val)

knitr::kable(ray_res)
c_gr ray_p_val
male_B 0.0513570
female_B 0.2785155
male_NB 0.0000018
female_NB 0.0000045
Not-uniform distribution detected mostly for non-breeders, both sexes

Comparing districutions

Pairwise comparisons with Watson two-sample test

c_gr <- unique(td_circ$group)
wats_stat <- numeric()

loop_vecA <- c(rep(c_gr[1], 3),
               rep(c_gr[2], 2),
               c_gr[3])

loop_vecB <- c(c_gr[2],
               rep(c(c_gr[3],c_gr[4]), 2),
               c_gr[4])

for(i in 1:length(loop_vecA)) {
 
  td_circ_A <- td_circ %>% filter(group == loop_vecA[i])
  td_circ_B <- td_circ %>% filter(group == loop_vecB[i])
  
  wats_st <- watson.two.test(td_circ_A$cHr, td_circ_B$cHr,alpha = 0.05)
  wats_stat[i] <- wats_st$statistic
   
}

wats_res_df <- data.frame(loop_vecA, loop_vecB, wats_stat)
crit_val <- 0.187 

wats_res_df$Sigf <- if_else(wats_res_df$wats_stat>=crit_val, "0.05", "ns.")


knitr::kable(wats_res_df)
loop_vecA loop_vecB wats_stat Sigf
male_B female_B 0.1848639 ns.
male_B male_NB 0.3991350 0.05
male_B female_NB 0.5517867 0.05
female_B male_NB 0.2188754 0.05
female_B female_NB 0.1960029 0.05
male_NB female_NB 0.1048484 ns.
Significant difference between birds of different status but not between sexes within each breedingstatus group