.Rhistory

#'
#' @param process_df A dataframe containing the dyad data with columns "winner", "loser", and "interactions".
#'
#' @return A dataframe with the total number of interactions for each dyad and the winning percentage.
total_interaction_per_dyad <- function(process_df){
# Create a combined dyad_id
process_df$dyad_id <- apply(process_df[, c("winner", "loser")], 1, function(x) paste(sort(as.integer(x)), collapse = "-"))
# remove when winner and loser are the same, also remove the dyads that have 0 interactions
process_df2 <- process_df[as.character(process_df$winner) != as.character(process_df$loser), ]
# calculate the total number of interaction per unique dyad A>B and A<B are the same dyad
interact_per_dyad <- aggregate(process_df2$interactions, by = list(process_df2$dyad_id), FUN = sum)
colnames(interact_per_dyad) <- c("dyad_id", "total_interactions")
process_df3 <-merge(process_df2, interact_per_dyad)
process_df3 <- process_df3[which(process_df3$total_interactions >0),]
process_df3$win_pct <- process_df3$interactions/process_df3$total_interactions
process_df3 <- process_df3[order(process_df3$dyad_id),]
return(process_df3)
}
#' Process Dyads to Set Dominant Winner
#'
#' This function processes a dataframe to ensure that for each dyad, there's only one record
#' where the cow that wins more is set as the winner. If there's a tie in the number of wins,
#' the cow with the larger ID is set as the winner.
#'
#' for each dyad, they have an unique dyad_id.
#' as I need to record the winner and loser for each dyad, I set the winner cow to be
#' the cow that wins more over the other, when this dyad appear for the first time in a feeder occupancy level
#'
#' @param process_df3 A dataframe containing the dyad data with columns "winner", "loser", "win_pct", and "dyad_id".
#'
#' @return A dataframe with processed dyads where each dyad has only one record with the dominant winner set.
low_fo_dyad_set <- function(process_df3) {
# now there are 2 record for each dyad, because A>B and A<B are in 2 seperate rows.
# For each dyad, only keep 1 record, the record where the cow that wins more is placed on the winner
process_df4 <- process_df3[which(process_df3$win_pct >=0.5),]
process_df4_no_tie <- process_df4[which(process_df4$win_pct > 0.5),]
#for dyads where 2 individuals wins the same amount, assign the winner to be the cow with larger cow ID
process_df4_tie <- process_df4[which(process_df4$win_pct == 0.5),]
process_df4_tie$winner <- NULL
process_df4_tie$loser <- NULL
process_df4_tie <- unique(process_df4_tie)
split_ids <- strsplit(process_df4_tie$dyad_id, "-") # Split the dyad_id column on the hyphen
process_df4_tie$winner <- sapply(split_ids, function(x) x[1])  # Extract the winner (first cow ID before the hyphen)
process_df4_tie$loser <- sapply(split_ids, function(x) x[2])  # Extract the loser (second cow ID after the hyphen)
process_df4_tie <- process_df4_tie[, colnames(process_df4_no_tie)]
process_df5 <- rbind(process_df4_no_tie, process_df4_tie)
return(process_df5)
}
#' Process Dyads Based on Previous Appearances
#'
#' This function processes a dataframe to determine the sequence of winners and losers for dyads
#' based on their previous appearances in the lowest level of feeder occupancy. If a dyad hasn't
#' appeared in previous levels, it uses the method from the `low_fo_dyad_set` function to determine the sequence.
#'
#' @param prog_df A dataframe containing the dyad data with columns "dyad_id", "winner", and "loser".
#' @param interactions_by_dyad A dataframe containing previous interactions by dyad.
#'
#' @return A dataframe with processed dyads based on their previous appearances.
other_fo_dyad_set <- function(prog_df, interactions_by_dyad) {
# for the dyads that have showed up in previous levels of feeder occupancy
# keep the sequence of who is the winner and who is the loser
dyad_seq <- unique(interactions_by_dyad[, c("dyad_id", "winner", "loser")])
dyad_showed <- prog_df[which(prog_df$dyad_id %in% dyad_seq$dyad_id),]
dyad_showed_processed <- merge(prog_df, dyad_seq)
# for those did not show up in previous levels of feeder occupancy
# record the dyad using the same method as the loest feeder occupancy level
dyad_not_showed <- prog_df[which(!(prog_df$dyad_id %in% dyad_seq$dyad_id)),]
if (nrow(dyad_not_showed) > 0) { # if there are new dyad show up
dyad_not_showed_processed <- low_fo_dyad_set(dyad_not_showed)
# order column names
dyad_showed_processed <- dyad_showed_processed[, colnames(dyad_not_showed_processed)]
prog_df_processed <- rbind(dyad_showed_processed, dyad_not_showed_processed)
} else {
prog_df_processed <- dyad_showed_processed
}
return(prog_df_processed)
}
#' Find Dyads Present in Only One Level of Feeder Occupancy not the other levels
#'
#' This function identifies dyads that only appear in one specific level of feeder occupancy
#' and not in other levels. It returns these dyads along with the specific level they appear in,
#' as well as a count of how many such dyads are present in each level.
#'
#' @param data A dataframe containing the dyad data.
#' @param dyad_id_col The column name in the dataframe that represents the dyad ID.
#' @param occupancy_col The column name in the dataframe that represents the feeder occupancy level.
#'
#' @return A list containing two dataframes:
#'   - single_level_dyads: A dataframe with dyads that only appear in one specific level of feeder occupancy.
#'   - single_level_dyads_count: A dataframe with a count of how many such dyads are present in each level.
find_dyads_in_single_level <- function(data, dyad_id_col, occupancy_col) {
# Group the data by dyad_id and feeder_occupancy
grouped_data <- aggregate(data[[occupancy_col]], by = list(data[[dyad_id_col]], data[[occupancy_col]]), FUN = length)
# Count the number of unique levels of feeder_occupancy for each dyad_id
dyad_counts <- aggregate(grouped_data$x, by = list(grouped_data$Group.1), FUN = length)
# Filter the dyad_id that have counts equal to 1
single_level_dyads <- dyad_counts[dyad_counts$x == 1, "Group.1"]
# Create an empty dataframe to store results
result_df <- data.frame(dyad_id = character(), feeder_occupancy = numeric())
# Create an empty dataframe to store the count of single_level_dyads for each level
count_df <- data.frame(feeder_occupancy = numeric(), count = integer())
# Loop through each unique level of feeder_occupancy
for (level in unique(data[[occupancy_col]])) {
# Filter the dyad_id that only show up in the current level
dyads_in_current_level <- grouped_data[grouped_data$Group.1 %in% single_level_dyads & grouped_data$Group.2 == level, "Group.1"]
# Create a dataframe with the filtered dyad_id and their corresponding feeder_occupancy level
current_level_df <- data.frame(dyad_id = dyads_in_current_level, feeder_occupancy = level)
# Append the current level dataframe to the result dataframe
result_df <- rbind(result_df, current_level_df)
# Record the count of single_level_dyads for the current level
count_df <- rbind(count_df, data.frame(feeder_occupancy = level, count = length(dyads_in_current_level)))
}
return(list(single_level_dyads = result_df, single_level_dyads_count = count_df))
}
#' Calculate replacements by Dyad at Different Feeder Occupancy Levels
#'
#' This function calculates the number of replacements that occurred per dyad
#' at each of the 25 levels of feeder occupancy.
#'
#' @param repl_master A dataframe containing interaction data with columns 'feeder_occupancy', 'winner', and 'loser'.
#' @param occupancy_col A string specifying the column name for occupancy (either "feeder_occupancy" or "feeder_occupancy_grouped").
#'
#' @return A dataframe containing interactions by dyad at different feeder occupancy levels.
calculate_interactions_by_dyad <- function(repl_master, occupancy_col) {
fed_occupancy_list <- unique(repl_master[[occupancy_col]])
contingency_tables <- list()
interactions_by_dyad <- data.frame()
for (i in 1:length(fed_occupancy_list)) {
# Subset data by feeder occupancy level
occupancy_data <- subset(repl_master, repl_master[[occupancy_col]] == fed_occupancy_list[i])
# Create contingency table for winner-loser pairs with the same unique values
occupancy_dyad_table <- table(factor(occupancy_data$winner, levels = sort(unique(c(occupancy_data$winner, occupancy_data$loser)))),
factor(occupancy_data$loser, levels = sort(unique(c(occupancy_data$winner, occupancy_data$loser)))))
# Transform the contingency table into a dataframe
temp_df <- contingency_to_dataframe(occupancy_dyad_table, fed_occupancy_list[i], occupancy_col)
# create unique dyad_id, and calculate total number of interactions/dyad, and wining percentage
temp_df3 <- total_interaction_per_dyad(temp_df)
# when it's the lowest level of feeder occupancy, for each dyad, only keep 1 record,
# the record where the cow that wins more is placed on the winner
if (i == 1) {
temp_df5 <- low_fo_dyad_set(temp_df3)
} else {
temp_df5 <- other_fo_dyad_set(temp_df3, interactions_by_dyad)
}
# Append merged dyads to the final dataframe
interactions_by_dyad <- rbind(interactions_by_dyad, temp_df5)
}
interactions_by_dyad[[occupancy_col]] <- round(interactions_by_dyad[[occupancy_col]], digits = 2)
return(interactions_by_dyad)
}
#' Calculate Total Dyads for Each Level of Feeder Occupancy
#'
#' This function computes the total number of dyads for each unique level of feeder occupancy.
#'
#' @param interactions_by_dyad A dataframe containing interaction data for different dyads across various levels of feeder occupancy.
#'
#' @return A dataframe with two columns:
#'   - feeder_occupancy: The unique levels of feeder occupancy.
#'   - total_dyad: The total number of dyads for each level of feeder occupancy.
calculate_total_dyads <- function(interactions_by_dyad){
temp = interactions_by_dyad
temp$c = 1
total_dyad <- aggregate(temp$c, by = list(temp$feeder_occupancy), FUN = sum)
colnames(total_dyad) <- c("feeder_occupancy", "total_dyad")
return(total_dyad)
}
#' Calculate Percentage of 2-Way Dyads for Dyads with 2 or More Interactions
#'
#' This function calculates the percentage of 2-way dyads for dyads with 2 or more interactions.
#' It filters the dyads based on the total number of interactions and then calculates the percentage
#' of 2-way interactions for each level of feeder occupancy.
#'
#' @param interactions_by_dyad A dataframe containing interaction data for different dyads across various levels of feeder occupancy.
#'
#' @return A dataframe with columns:
#'   - feeder_occupancy: The unique levels of feeder occupancy.
#'   - dyads_mt2_interact: The total number of dyads with more than 2 interactions.
#'   - 2way_dyad: The total number of 2-way dyads.
#'   - 2way_pct: The percentage of 2-way dyads.
two_way_dyad_pct_calculation <- function(interactions_by_dyad) {
mt_2_dyad <- interactions_by_dyad[which(interactions_by_dyad$total_interactions >=2),]
mt_2_dyad_temp <- mt_2_dyad
mt_2_dyad_temp$c = 1
mt_2_dyad_total <- aggregate(mt_2_dyad_temp$c, by = list(mt_2_dyad_temp$feeder_occupancy), FUN = sum)
colnames(mt_2_dyad_total) <- c("feeder_occupancy", "dyads_mt2_interact")
mt_2_2way_dyad <- mt_2_dyad[which((abs(mt_2_dyad$win_pct) != 1 ) & (abs(mt_2_dyad$win_pct) != 0)),]
mt_2_2way_dyad$c = 1
mt_2_2way_dyad_total <- aggregate(mt_2_2way_dyad$c, by = list(mt_2_2way_dyad$feeder_occupancy), FUN = sum)
colnames(mt_2_2way_dyad_total) <- c("feeder_occupancy", "2way_dyad")
two_way_dyad_perct <- merge(mt_2_dyad_total, mt_2_2way_dyad_total)
two_way_dyad_perct$`2way_pct` <- two_way_dyad_perct$`2way_dyad`/two_way_dyad_perct$dyads_mt2_interact
return(two_way_dyad_perct)
}
#' Plot Percentage of Two-way Dyads by Feeder Occupancy
#'
#' This function creates a scatter plot to visualize the percentage of two-way dyads among dyads
#' with 2 or more interactions across different levels of feeder occupancy.
#'
#' @param dyad_summary2 A dataframe containing the summary of dyads with columns 'feeder_occupancy' and '2way_pct'.
#'
#' @return A ggplot object visualizing the percentage of two-way dyads by feeder occupancy.
#'
two_way_pct_by_feeder_occupancy_plot <- function(dyad_summary2) {
two_way_pct_plot <- ggplot(dyad_summary2, aes(x=feeder_occupancy, y = `2way_pct`)) +
geom_point(aes(y = `2way_pct`), size = 10, color = "royal blue") +
geom_smooth(method = "lm", se = FALSE, size= 2, color = "midnight blue", fullrange = TRUE) +
labs(y= "Percentage of Two-way Dyads", x = "Feeder Occupancy") +
theme_classic() +
theme(text = element_text(size = 55), axis.text.x = element_text(size = 50)) +
scale_y_continuous(expand=expansion(mult = c(0, .1)), limits = c(0, 0.63))
return(two_way_pct_plot)
}
#' Save Two-way Dyads Percentage Plot by Feeder Occupancy
#'
#' This function generates a scatter plot visualizing the percentage of two-way dyads among dyads
#' with 2 or more interactions across different levels of feeder occupancy and saves it to a specified directory.
#'
#' @param dyad_summary2 A dataframe containing the summary of dyads with columns 'feeder_occupancy' and '2way_pct'.
#' @param output_dir A string specifying the directory where the plot should be saved.
two_way_pct_plot <- function(dyad_summary2) {
two_way_pct_plot <- two_way_pct_by_feeder_occupancy_plot(dyad_summary2)
file_name = here("graphs/2way_pct_by_feeder_occupancy.png")
ggsave(file_name, plot = two_way_pct_plot, width = 15, height = 13, limitsize = FALSE)
}
repl_master <- replacement_sampled_master[, c("Actor_cow", "Reactor_cow", "end_density")]
colnames(df = repl_master) <- c("winner", "loser", "feeder_occupancy")
repl_master_grouped <- group_density_buckets(repl_master, "feeder_occupancy")
interactions_by_dyad_grouped <- calculate_interactions_by_dyad(repl_master_grouped, "feeder_occupancy_grouped")
View(interactions_by_dyad_grouped)
# find dyads that show up in all 5 levels of feeder occupancy
dyads_in_all_levels_grouped <- find_dyads_in_all_levels(interactions_by_dyad_grouped, "dyad_id", "feeder_occupancy_grouped")
interactions_by_dyad_same_across_all_levels_grouped <- interactions_by_dyad_grouped[which(interactions_by_dyad_grouped$dyad_id %in% dyads_in_all_levels_grouped),]
dyads_in_all_levels_grouped_num <- length(dyads_in_all_levels_grouped)
win_pct_changes <- calculate_win_pct_change_per_dyad(
interactions_by_dyad_same_across_all_levels_grouped,
"dyad_id",
"feeder_occupancy_grouped",
"win_pct"
)
# the changes in win pct is all relative to the lowest level of feeder occupancy after grouping (which is 0% - 27%)
win_pct_changes2 <- win_pct_changes[which(win_pct_changes$feeder_occupancy_grouped != 0.27),]
win_pct_changes2_0 <- win_pct_changes2[which(win_pct_changes2$win_pct_change ==0),]
win_pct_changes2_pos <- win_pct_changes2[which(win_pct_changes2$win_pct_change >0),]
win_pct_changes2_neg <- win_pct_changes2[which(win_pct_changes2$win_pct_change <0),]
# plot dyads with positive and negative changes
positive_win_pct_change_violin_boxplot(win_pct_changes2_pos)
negative_win_pct_change_violin_boxplot(win_pct_changes2_neg)
# Fit a mixed-effects model
win_pct_changes_pos <- lmerTest::lmer(win_pct_change ~ feeder_occupancy_grouped + (1 | dyad_id), data = win_pct_changes2_pos)
summary(win_pct_changes_pos)
win_pct_changes_neg <- lmerTest::lmer(win_pct_change ~ feeder_occupancy_grouped + (1 | dyad_id), data = win_pct_changes2_neg)
summary(win_pct_changes_neg)
# calculate the percentage of dyads with 0 change in win_pct across all levels
total_same_dyad <- length(unique(win_pct_changes2$dyad_id))
grouped_df <- dplyr::group_by(win_pct_changes2_0, feeder_occupancy_grouped)
unique_dyad_counts <- dplyr::summarise(grouped_df, n_unique_dyads = dplyr::n_distinct(dyad_id))
colnames(unique_dyad_counts) <- c("feeder_occupancy_grouped", "dyads_0_change")
unique_dyad_counts$dyads_0_change_pct <- unique_dyad_counts$dyads_0_change/total_same_dyad
View(unique_dyad_counts)
# calculate the percentage of dyads with negative changes in win pct in win_pct across all levels
total_same_dyad <- length(unique(win_pct_changes2$dyad_id))
grouped_df_neg <- dplyr::group_by(win_pct_changes2_neg, feeder_occupancy_grouped)
unique_dyad_counts_neg <- dplyr::summarise(grouped_df_neg, n_unique_dyads = dplyr::n_distinct(dyad_id))
colnames(unique_dyad_counts_neg) <- c("feeder_occupancy_grouped", "dyads_neg_change")
unique_dyad_counts_neg$dyads_neg_change_pct <- unique_dyad_counts_neg$dyads_neg_change / total_same_dyad
# calculate the percentage of dyads with positive changes in win pct in win_pct across all levels
total_same_dyad <- length(unique(win_pct_changes2$dyad_id))
grouped_df_pos <- dplyr::group_by(win_pct_changes2_pos, feeder_occupancy_grouped)
unique_dyad_counts_pos <- dplyr::summarise(grouped_df_pos, n_unique_dyads = dplyr::n_distinct(dyad_id))
colnames(unique_dyad_counts_pos) <- c("feeder_occupancy_grouped", "dyads_pos_change")
unique_dyad_counts_pos$dyads_pos_change_pct <- unique_dyad_counts_pos$dyads_pos_change / total_same_dyad
View(unique_dyad_counts_neg)
View(unique_dyad_counts_pos)
minute_thereshold <- 30
second_thereshold <- minute_thereshold * 60
percent_longer_than30 <- nrow(master_feeding[which(master_feeding$Duration > second_thereshold),])/nrow(master_feeding)
data_summary_dur <- summary(master_feeding$Duration)
Q1_dur <-data_summary_dur[["1st Qu."]]
Q3_dur <-data_summary_dur[["3rd Qu."]]
IQR_dur <-Q3_dur - Q1_dur
upper_whisker_dur <-  Q3_dur+(1.5)*IQR_dur
lower_whisker_dur <-  Q1_dur-(1.5)*IQR_dur
dur_ourliter_percent <- nrow(master_feeding[which(master_feeding$Duration > upper_whisker_dur),])/nrow(master_feeding)
# large feed intake
data_summary_intake <- summary(master_feeding$Intake)
Q1_intake <-data_summary_intake[["1st Qu."]]
Q3_intake <-data_summary_intake[["3rd Qu."]]
IQR_intake <-Q3_intake - Q1_intake
upper_whisker_intake <-  Q3_intake+(1.5)*IQR_intake
lower_whisker_intake <-  Q1_intake-(1.5)*IQR_intake
intake_ourliter_percent <- nrow(master_feeding[which(master_feeding$Intake > upper_whisker_intake),])/nrow(master_feeding)
# double detections
load(here("data/results/double_detection_1cow_2bins.rda"))
for (i in 1:length(double_detection_1cow_2bins)) {
if (i == 1) {
double_detection_master <- double_detection_1cow_2bins[[i]]
} else {
double_detection_master <- rbind(double_detection_master, double_detection_1cow_2bins[[i]])
}
}
# double detections
load(here("data/results/double_detection_1cow_2bins.rda"))
for (i in 1:length(double_bin_detection_list)) {
if (i == 1) {
double_detection_master <- double_bin_detection_list[[i]]
} else {
double_detection_master <- rbind(double_detection_master, double_bin_detection_list[[i]])
}
}
double_detection_percent <- nrow(double_detection_master)/nrow(master_feeding)
# calculate the average and SD of days each cow stayed in the trial
cow_day_cal <- unique(master_comb[, c("Cow", "date")])
cow_day_cal$count <- 1
cow_day_sum <- aggregate(cow_day_cal$count, by = list(cow_day_cal$Cow), FUN = sum)
colnames(cow_day_sum) <- c("Cow", "total_days")
cow_day_mean <- mean(cow_day_sum$total_days)
cow_day_sd <- sd(cow_day_sum$total_days)
### calculate the average time for morning feed delivery
load(here("data/results/feed_delivery.rda"))
feed_delivery$date <- ymd(feed_delivery$date, tz=time_zone)
load("C:/Users/skysheng/OneDrive - UBC/R package project and Git/competition_dominance_analysis/data/results/feed_delivery.rda")
View(time_interval_after_feed_added)
### calculate the average time for morning feed delivery
load(here("data/results/feed_delivery.rda"))
time_interval_after_feed_added$date <- ymd(time_interval_after_feed_added$date, tz=time_zone)
master_feed_replacement_all$date <- ymd(master_feed_replacement_all$date, tz=time_zone)
time_interval_after_feed_added <- time_interval_after_feed_added[which(time_interval_after_feed_added$date %in% master_feed_replacement_all$date),]
# Filter rows with time between 4 am and 8 am
time_interval_after_feed_added_filtered <- time_interval_after_feed_added %>% filter(hour(morning_feed_add_start) >= 4 & hour(morning_feed_add_start) < 8)
# Extract the time part from morning_feed_add_start
time_interval_after_feed_added_filtered <- time_interval_after_feed_added_filtered %>% mutate(morning_time_part = as.numeric(format(morning_feed_add_start, format = "%H")) * 60 + as.numeric(format(morning_feed_add_start, format = "%M")))
# Calculate the average time in minutes
average_time <- mean(time_interval_after_feed_added_filtered$morning_time_part)
# Convert the average time to hours and minutes format
average_hh_mm_am <- sprintf("%02d:%02d", average_time %/% 60, round(average_time %% 60))
# Print the average time expressed in hh:mm
average_hh_mm_am
### calculate the average time for afternoon feed delivery
# Define the desired time range (e.g., 1 pm to 5 pm)
start_hour <- 13
end_hour <- 17
# Filter rows with time within the desired range
time_interval_after_feed_added_filtered <- time_interval_after_feed_added %>% filter(hour(afternoon_feed_add_start) >= start_hour & hour(afternoon_feed_add_start) < end_hour)
# Extract the time part from afternoon_feed_add_start
time_interval_after_feed_added_filtered <- time_interval_after_feed_added_filtered %>% mutate(afternoon_time_part = as.numeric(format(afternoon_feed_add_start, format = "%H")) * 60 + as.numeric(format(afternoon_feed_add_start, format = "%M")))
# Calculate the average time in minutes
average_time <- mean(time_interval_after_feed_added_filtered$afternoon_time_part)
# Convert the average time to hours and minutes format
average_hh_mm_pm <- sprintf("%02d:%02d", average_time %/% 60, round(average_time %% 60))
# Print the average time expressed in hh:mm
average_hh_mm_pm
# calculate the number of days included at each level
load(here("data/results/replacement_sampled_master.rda"))
unique_day_level <- unique(replacement_sampled_master[, c("date", "start_density", "end_density")])
unique_day_level$count <- 1
unique_day_sum <- aggregate(unique_day_level$count, by= list(unique_day_level$end_density), FUN = sum)
colnames(unique_day_sum) <- c("end_occupancy", "total_days")
# Print the average time expressed in hh:mm
average_hh_mm_pm
### calculate the average time for morning feed delivery
load(here("data/results/feed_delivery.rda"))
time_interval_after_feed_added$date <- ymd(time_interval_after_feed_added$date, tz=time_zone)
master_feed_replacement_all$date <- ymd(master_feed_replacement_all$date, tz=time_zone)
time_interval_after_feed_added <- time_interval_after_feed_added[which(time_interval_after_feed_added$date %in% master_feed_replacement_all$date),]
# Filter rows with time between 4 am and 8 am
time_interval_after_feed_added_filtered <- time_interval_after_feed_added %>% filter(hour(morning_feed_add_start) >= 4 & hour(morning_feed_add_start) < 8)
# Extract the time part from morning_feed_add_start
time_interval_after_feed_added_filtered <- time_interval_after_feed_added_filtered %>% mutate(morning_time_part = as.numeric(format(morning_feed_add_start, format = "%H")) * 60 + as.numeric(format(morning_feed_add_start, format = "%M")))
# Calculate the average time in minutes
average_time <- mean(time_interval_after_feed_added_filtered$morning_time_part)
# Convert the average time to hours and minutes format
average_hh_mm_am <- sprintf("%02d:%02d", average_time %/% 60, round(average_time %% 60))
# Print the average time expressed in hh:mm
average_hh_mm_am
### calculate the average time for afternoon feed delivery
# Define the desired time range (e.g., 1 pm to 5 pm)
start_hour <- 13
end_hour <- 17
# Filter rows with time within the desired range
time_interval_after_feed_added_filtered <- time_interval_after_feed_added %>% filter(hour(afternoon_feed_add_start) >= start_hour & hour(afternoon_feed_add_start) < end_hour)
# Extract the time part from afternoon_feed_add_start
time_interval_after_feed_added_filtered <- time_interval_after_feed_added_filtered %>% mutate(afternoon_time_part = as.numeric(format(afternoon_feed_add_start, format = "%H")) * 60 + as.numeric(format(afternoon_feed_add_start, format = "%M")))
# Calculate the average time in minutes
average_time <- mean(time_interval_after_feed_added_filtered$afternoon_time_part)
# Convert the average time to hours and minutes format
average_hh_mm_pm <- sprintf("%02d:%02d", average_time %/% 60, round(average_time %% 60))
# Print the average time expressed in hh:mm
average_hh_mm_pm
load("C:/Users/skysheng/OneDrive - UBC/R package project and Git/competition_dominance_analysis/data/results/Bins_with_number_of_visits_daily.rda")
load("C:/Users/skysheng/OneDrive - UBC/R package project and Git/competition_dominance_analysis/data/results/Bins_with_number_of_visits_daily.rda")
View(bins_visit_num)
cache("bins_visit_num")
load("C:/Users/skysheng/Desktop/temp/Feeding and drinking analysis.rda")
cache("Insentec_final_summary")
save(bins_visit_num, file = (here::here(paste0("data/results/", "Bins_with_number_of_visits_daily.rdata"))))
save(Insentec_final_summary, file = (here::here(paste0("data/results/", "Feeding and drinking analysis.rdata"))))
load("C:/Users/skysheng/OneDrive - UBC/R package project and Git/competition_dominance_analysis/data/results/replacement_sampled_master.rda")
save(replacement_sampled_master, file = here("data/results/replacement_sampled_master.rdata"))
load("C:/Users/skysheng/Desktop/temp/number_of_bins_visited_by_each_cow.rda")
save(bin_num_visit_per_cow, file = (here::here(paste0("data/results/", "number_of_bins_visited_by_each_cow.rdata"))))
load("C:/Users/skysheng/Desktop/temp/number_of_visits_for_each_bin_for_each_cow.rda")
save(visit_per_bin_per_cow, file = (here::here(paste0("data/results/", "number_of_visits_for_each_bin_for_each_cow.rdata"))))
save(visit_per_bin_per_cow, file = (here::here(paste0("data/results/", "number_of_visits_for_each_bin_for_each_cow.rdata"))))
library(ProjectTemplate)
load.project()
View(warning_days)
save(warning_days, here("data/warning_days.rdata"))
here()
here("data")
save(warning_days, file=here("data/warning_days.rdata"))
load.project()
load.project()
View(warning_days)
View(daylight_saving_csv)
save(saylight_saving_csv, file=here("data/daylight_saving_csv.rdata"))
save(daylight_saving_csv, file=here("data/daylight_saving_csv.rdata"))
(30*29)/2
library('ProjectTemplate')
load.project()
total_dyad_possible_num <- total_dyad_long_term(all.comb2)
group_by <- 1
master_dyad_analysis_10mon_control <- dyad_relationship_long_term(replacement_sampled_master, res_occupancy_seq, group_by, total_dyad_possible_num)
master_unique_dyad_direction_10mon_control <- master_dyad_analysis_10mon_control[[1]]
dyads_unknown <- master_dyad_analysis_10mon_control[[2]]
plot_unknown(dyads_unknown, here("graphs/"))
# linear model fit
unknown_lm <- lm(percent_unkown ~ end_density, data = dyads_unknown)
unknown_lm_summary <- summary(unknown_lm)
################################################################################
############### Steepness changes as feeder occupancy increases ################
################################################################################
# plot the steepness of elo under different Feeder Occupancy
master_steepness2 <- master_steepness
elo_steepness_plot(master_steepness)
# fit a linear model for Elo Steepness X Feeder Occupancy
steepness_lm <- lm(master_steepness$steepness_mean ~ master_steepness$resource_occupancy)
steepness_lm_summary <- summary(steepness_lm)
steepness_lm_summary
repl_master <- replacement_sampled_master[, c("Actor_cow", "Reactor_cow", "end_density")]
colnames(repl_master) <- c("winner", "loser", "feeder_occupancy")
interactions_by_dyad <- calculate_interactions_by_dyad(repl_master, "feeder_occupancy")
# find dyads that show up in all levels of feeder occupancy
dyads_in_all_levels <- find_dyads_in_all_levels(interactions_by_dyad, "dyad_id", "feeder_occupancy")
# find dyads that only show up in certain levels
results <- find_dyads_in_single_level(interactions_by_dyad, "dyad_id", "feeder_occupancy")
single_level_dyads_df <- results$single_level_dyads
single_level_dyads_count_df <- results$single_level_dyads_count
colnames(single_level_dyads_count_df) <- c("feeder_occupancy", "unique_dyad_num")
total_dyad <- calculate_total_dyads(interactions_by_dyad)
single_level_dyads_count_df2 <- merge(single_level_dyads_count_df, total_dyad)
single_level_dyads_count_df2$unique_dyad_pct <- single_level_dyads_count_df2$unique_dyad_num/single_level_dyads_count_df2$total_dyad
# linear model fit for unique dyad per feeder occupancy level
single_level_dyads_lm <- lm(unique_dyad_pct ~ feeder_occupancy, data = single_level_dyads_count_df2)
single_level_dyads_lm_summary <- summary(single_level_dyads_lm)
################################################################################
########## Dyad Level Analysis: percentage of two way dyads among dyads ########
######################## with >= 2 interactions ################################
################### control for same number of replacements ####################
################################################################################
two_way_dyad_perct <- two_way_dyad_pct_calculation(interactions_by_dyad)
two_way_dyad_perct
# plot the changes of percentage of 2-way dyads (among dyads with >=2 interactions)
# as feeder occupancy increases
two_way_pct_plot(two_way_dyad_perct)
single_level_dyads_lm_summary
View(two_way_dyad_perct)
str(two_way_dyad_perct)
# linear model fit for two-way dyad per feeder occupancy level
two_way_dyads_lm <- lm(`2way_pct` ~ feeder_occupancy, data = two_way_dyad_perct)
two_way__dyads_lm_summary <- summary(two_way_dyads_lm)
two_way__dyads_lm_summary
repl_master <- replacement_sampled_master[, c("Actor_cow", "Reactor_cow", "end_density")]
colnames(df = repl_master) <- c("winner", "loser", "feeder_occupancy")
repl_master_grouped <- group_density_buckets(repl_master, "feeder_occupancy")
interactions_by_dyad_grouped <- calculate_interactions_by_dyad(repl_master_grouped, "feeder_occupancy_grouped")
# find dyads that show up in all 5 levels of feeder occupancy
dyads_in_all_levels_grouped <- find_dyads_in_all_levels(interactions_by_dyad_grouped, "dyad_id", "feeder_occupancy_grouped")
interactions_by_dyad_same_across_all_levels_grouped <- interactions_by_dyad_grouped[which(interactions_by_dyad_grouped$dyad_id %in% dyads_in_all_levels_grouped),]
dyads_in_all_levels_grouped_num <- length(dyads_in_all_levels_grouped)
win_pct_changes <- calculate_win_pct_change_per_dyad(
interactions_by_dyad_same_across_all_levels_grouped,
"dyad_id",
"feeder_occupancy_grouped",
"win_pct"
)
# the changes in win pct is all relative to the lowest level of feeder occupancy after grouping (which is 0% - 27%)
win_pct_changes2 <- win_pct_changes[which(win_pct_changes$feeder_occupancy_grouped != 0.27),]
win_pct_changes2_0 <- win_pct_changes2[which(win_pct_changes2$win_pct_change ==0),]
win_pct_changes2_pos <- win_pct_changes2[which(win_pct_changes2$win_pct_change >0),]
win_pct_changes2_neg <- win_pct_changes2[which(win_pct_changes2$win_pct_change <0),]
# plot dyads with positive and negative changes
positive_win_pct_change_violin_boxplot(win_pct_changes2_pos)
negative_win_pct_change_violin_boxplot(win_pct_changes2_neg)
# Fit a mixed-effects model
control <- lmerControl(optimizer = "bobyqa")
win_pct_changes_pos <- lmerTest::lmer(win_pct_change ~ feeder_occupancy_grouped + (feeder_occupancy_grouped | dyad_id), data = win_pct_changes2_pos, control = control)
summary(win_pct_changes_pos)
win_pct_changes_neg <- lmerTest::lmer(win_pct_change ~ feeder_occupancy_grouped + (feeder_occupancy_grouped | dyad_id), data = win_pct_changes2_neg, control = control)
summary(win_pct_changes_neg)
library('ProjectTemplate')
load.project()
##
total_dyad_possible_num <- total_dyad_long_term(all.comb2)
group_by <- 1
master_dyad_analysis_10mon_control <- dyad_relationship_long_term(replacement_sampled_master, res_occupancy_seq, group_by, total_dyad_possible_num)
master_unique_dyad_direction_10mon_control <- master_dyad_analysis_10mon_control[[1]]
dyads_unknown <- master_dyad_analysis_10mon_control[[2]]
plot_unknown(dyads_unknown, here("graphs/"))
# linear model fit
unknown_lm <- lm(percent_unkown ~ end_density, data = dyads_unknown)
unknown_lm_summary <- summary(unknown_lm)
unknown_lm_summary