p_discolor_auditory_feedback/vocalisations/scripts/analysis.R

library(tidyverse)
library(lubridate)
library(nlme)

birthdates <- tribble(~calling_bat, ~calling_bat_birthdate,
                      "b2", "2017-01-26",
                      "b4", "2017-01-30",
                      "b7", "2017-02-08",
                      "b3", "2017-01-29",
                      "b5", "2017-02-02",
                      "b8", "2017-02-08")
birthdates <- birthdates %>%
    mutate(calling_bat_birthdate=ymd(calling_bat_birthdate))

load_pup_data <- function(filename) {
    data <- read_csv(filename, na="--",
                     col_types=cols(.default="?", session_id="f",
                                    calling_bat="f", calling_bat_deaf="l", before_deafening="l"))
    data <- inner_join(data, birthdates)

    # annotate with ages in weeks
    data <- data %>%
        mutate(age_in_weeks=
               floor(as.numeric(as.duration(data$calling_bat_birthdate %--% data$start_time),
                                "weeks")),
               bandwidth=f_max - f_min) %>%
        filter(!before_deafening)

    data
}

load_pup_activity_data <- function(filename) {
    data <- read_csv(filename,
                     col_types=cols(.default="?", session_id="f",
                                    calling_bat="f", calling_bat_deaf="l", before_deafening="l")) %>%
        rename(age_in_weeks=calling_bat_age_in_weeks)

    data
}

pup_data <- load_pup_data("../output/pup_calls_for_model.csv")
pup_activity_data <- load_pup_activity_data("../output/pup_activity_for_model.csv")

fit_pup_model <- function(data, response) {
    frm <- as.formula(paste(response, "~ calling_bat_deaf*age_in_weeks"))
    lme(frm, random=~1 + age_in_weeks | calling_bat,
        data=data, control=lmeControl(opt="optim"))
}

anova_p <- function(model) {
    result <- as.data.frame(anova(model))
    as.data.frame.list(c(p_group=result[["calling_bat_deaf", "p-value"]],
			 p_age=result[["age_in_weeks", "p-value"]],
			 p_interaction=result[["calling_bat_deaf:age_in_weeks", "p-value"]]))
}

parameters <- c("call_rms", "call_duration", "f_min", "f_max", "bandwidth",
                "f0_mean", "f0_min", "f0_max", "f0_start", "f0_end", "f0_slope",
                "spectral_centroid", "mean_aperiodicity")

pup_p <- bind_rows(lapply(parameters,
                            function(parameter)
                                anova_p(fit_pup_model(pup_data, parameter))),
                   anova_p(fit_pup_model(pup_activity_data, "calls_per_10s"))) %>%
    mutate(p_group_adjusted=p.adjust(p_group, "BH"),
	   p_age_adjusted=p.adjust(p_age, "BH"),
	   p_interaction_adjusted=p.adjust(p_age, "BH"),
           parameter=c(parameters, "calls_per_10s"))
write_csv(pup_p, "../output/pup_model_results.csv")