# 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")
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
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
# 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
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
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