## Packages ####
if (!require("pacman"))
  install.packages("pacman")
pacman::p_load(lme4, emmeans, afex, psych, GPArotation, pbkrtest, MASS, piecewiseSEM, grid, readxl, ggpubr, tidyverse)
## Data ####
humanmacaque <-
  read_csv(
    "humanmacaque.csv"
    )

longhumanmacaque <-
  read_csv(
    "longhumanmacaque.csv"
    ) %>% 
  mutate(group = factor(group),
         specmir = factor(specmir))

## Principal Component Analysis - Human Infant Behaviour ####
chhbehav <-
  humanmacaque %>%
  filter(
    group == "Humans"
    ) %>%
  select(
    ca1,
    ca2,
    cb1,
    cc1,
    cneg
    ) %>%
  drop_na()

chhKMO <- chhbehav %>% KMO(.)
chhbartlett <- chhbehav %>% cortest.bartlett(., n = 100)
chhpa <- chhbehav %>% fa.parallel(., fa = "pc", n.iter = 5000)
chhpca <- chhbehav %>% principal(., chhpa$ncomp, rotate = "simplimax")

## Principal Component Analysis - Human Maternal Responses ####
mhresp <-
  humanmacaque %>%
  filter(
    group == "Humans"
    ) %>%
  select(
    m_e,
    a_e,
    nresp
    ) %>%
  drop_na()

mhrKMO <- mhresp %>% KMO(.)
mhrbartlett <- mhresp %>% cortest.bartlett(., n = 100)
mhrpa <- mhresp %>% fa.parallel(., fa = "pc", n.iter = 5000)
mhrpca <- mhresp %>% principal(., mhrpa$ncomp, rotate = "simplimax")

## Principal Component Analysis - Macaque Infant Behaviour ####
chmbehav <-
  humanmacaque %>%
  filter(
    group == "Rhesus Macaques"
    ) %>%
  select(
    ca1,
    ca2,
    ci3, 
    cneg
    ) %>% drop_na()

chmKMO <- chmbehav %>% KMO(.)
chmbartlett <- chmbehav %>% cortest.bartlett(., n = 100)
chmpa <- chmbehav %>% fa.parallel(., fa = "pc", n.iter = 5000)
chmpca <- chmbehav %>% principal(., chmpa$ncomp, rotate = "simplimax")

## Principal Component Analysis - Macaque Maternal Responses ####
mmresp <-
  humanmacaque %>%
  filter(
    group == "Rhesus Macaques"
    ) %>%
  select(
    m_e,
    a_e
    ) %>%
  drop_na()

mmrKMO <- mmresp %>% KMO(.)
mmrbartlett <- mmresp %>% cortest.bartlett(., n = 100)
mmrpa <- mmresp %>% fa.parallel(., fa = "pc", n.iter = 5000)
mmrpca <- mmresp %>% principal(., mmrpa$ncomp, rotate = "simplimax")

## Gaze ####

gaze_mixed <-
  humanmacaque %>%
  mixed(
    cd1 ~ age_allweeks * group + (1 | id),
    offset = log(humanmacaque %>% pull(seconds)),
    data = .,
    family = "poisson",
    method = "LRT")

gaze_humans <-
  humanmacaque %>%
  glmer(
    cd1 ~ age_allweeks * group + (1 | id) + offset(log(seconds)),
    data = .,
    family = "poisson") %>%
  summary()

gaze_macaques <-
  humanmacaque %>%
  mutate(
    group = factor(group) %>% relevel(., "Rhesus Macaques")
    ) %>%
  glmer(
    cd1 ~ age_allweeks * group + (1 | id) + offset(log(seconds)),
    data = .,
    family = "poisson") %>%
  summary()

gaze_to_mother <-
  humanmacaque %>%
    transmute(
      group,
      perc = (cd1 * 100) / seconds,
      days =
        (age_allweeks *
           humanmacaque %>% filter(group == "Humans") %>% pull(days) %>% sd(., na.rm = T)) +
        humanmacaque %>% filter(group == "Humans") %>% pull(days) %>% mean(., na.rm = T)) %>%
    mutate(group = factor(ifelse(group == "Humans", "Humans", "Rhesus Macaque"))) %>%
    ggplot(data = .,
           aes(x=days, y=perc, color=group)) + geom_point() +
    scale_colour_hue(l=50) +
    geom_smooth(method=lm,
                se=FALSE,
                size = 2,
                fullrange=TRUE) +
    scale_y_continuous(name = "Gaze to Mother\n(as Percentage of Interaction Time)") +
    scale_x_continuous(name = "Age in Days (Humans)",
                       sec.axis =
                         sec_axis(~ ((. - mean(humanmacaque %>% filter(group == "Humans") %>% pull(days), na.rm = T)) /
                                       sd(humanmacaque %>% filter(group == "Humans") %>% pull(days), na.rm = T)) *
                                    sd(humanmacaque %>% filter(group == "Rhesus Macaques") %>% pull(days), na.rm = T) +
                                    mean(humanmacaque %>% filter(group == "Rhesus Macaques") %>% pull(days), na.rm = T),
                                  name = "Age in Days (Rhesus Macaques)",
                                  breaks = seq(0, 15, 2))) +
    theme(
      axis.title = element_text(face = "bold", size = 16),
      axis.text = element_text(face = "bold", colour = "black", size = 16),
      axis.ticks = element_blank(),
      legend.key = element_rect(fill = "white", colour = "grey80"),
      legend.text = element_text(face = "bold", size = 10),
      legend.title = element_blank(),
      plot.title = element_text(face = "bold", size = 18),
      panel.background = element_rect(fill = "transparent", colour = NA),
      panel.grid.major  = element_line(colour = "black", size = 0.5),
      panel.grid.minor  = element_line(colour = "grey70", size = 0.2),
      plot.background = element_rect(fill = "transparent", colour = NA),
      strip.background  = element_rect(fill = "transparent", colour = NA),
      strip.text = element_text(colour = "black", face = "bold", size = 24)
    )




## Social Expressiveness ####

social_mixed <-
  humanmacaque %>%
  filter(
    !is.na(age_lastweeks),
    !is.na(csb),
    callall > 0) %>%
  mixed(
    csb ~ group * age_lastweeks + (1 | id),
    offset =
      log(
        humanmacaque %>%
          filter(
            !is.na(age_lastweeks),
            !is.na(csb),
            callall > 0) %>%
          pull(callall)
        ),
    data = .,
    family = "poisson",
    method = "LRT")

social_min <-
  humanmacaque %>%
  filter(
    !is.na(age_lastweeks),
    !is.na(csb),
    callall > 0) %>%
  mutate(
    age_lastweeks = age_lastweeks - (humanmacaque %>% filter(!is.na(age_lastweeks)) %>% pull(age_lastweeks) %>% min(.))
  ) %>%
  glmer(
    csb ~ group * age_lastweeks + (1 | id) + offset(log(callall)),
    data = .,
    family = "poisson") %>%
  summary()

social_max <-
  humanmacaque %>%
  filter(
    !is.na(age_lastweeks),
    !is.na(csb),
    callall > 0) %>%
  mutate(
    age_lastweeks = age_lastweeks - (humanmacaque %>% filter(!is.na(age_lastweeks)) %>% pull(age_lastweeks) %>% max(.))
    ) %>%
  glmer(
    csb ~ group * age_lastweeks + (1 | id) + offset(log(callall)),
    data = .,
    family = "poisson") %>%
  summary()

social_human <-
  humanmacaque %>%
  filter(
    !is.na(age_lastweeks),
    !is.na(csb),
    callall > 0) %>%
  mutate(
    group = ifelse(group == "Humans", 0, 1)
    ) %>%
  glmer(
    csb ~ group * age_lastweeks + (1 | id) + offset(log(callall)),
    data = .,
    family = "poisson") %>%
  summary()

social_macaques <-
  humanmacaque %>%
  filter(
    !is.na(age_lastweeks),
    !is.na(csb),
    callall > 0) %>%
  mutate(
    group = ifelse(group == "Humans", 1, 0)
    ) %>%
  glmer(
    csb ~ group * age_lastweeks + (1 | id) + offset(log(callall)),
    data = .,
    family = "poisson") %>%
  summary()

infant_social_behaviour<-
  humanmacaque %>%
  filter(
    !is.na(age_lastweeks),
    !is.na(csb),
    callall > 0) %>%
    mutate(perc = (csb *100) / callall,
           group,
      days =
             age_lastweeks *
             sd(humanmacaque %>% filter(group == "Humans", !is.na(age_lastweeks)) %>% pull(days), na.rm = T) +
             mean(humanmacaque %>% filter(group == "Humans", !is.na(age_lastweeks)) %>% pull(days), na.rm = T)) %>%
    ggplot(
      data = .,
      aes(x = days, y = perc, color = group)) +
    geom_point() +
    geom_smooth(
      method = lm,
      se = FALSE,
      size = 2,
      fullrange = TRUE
    ) +
    scale_y_continuous(
      name = "Infant Social Behaviour\n(as Percentage of all Behaviours)") +
    scale_x_continuous(
      name = "Age in Days (Humans)",
      sec.axis =
        sec_axis(~ ((. - mean(humanmacaque %>% filter(group == "Humans", !is.na(age_lastweeks)) %>% pull(days), na.rm = T)) /
                      sd(humanmacaque %>% filter(group == "Humans", !is.na(age_lastweeks)) %>% pull(days), na.rm = T)) *
                   sd(humanmacaque %>% filter(group == "Rhesus Macaques") %>% pull(days), na.rm = T) +
                   mean(humanmacaque %>% filter(group == "Rhesus Macaques") %>% pull(days), na.rm = T),
                 name = "Age in Days (Rhesus Macaques)",
                 breaks = seq(0, 15, 2))) +
  theme(
    axis.title = element_text(face = "bold", size = 16),
    axis.text = element_text(face = "bold", colour = "black", size = 16),
    axis.ticks = element_blank(),
    legend.key = element_rect(fill = "white", colour = "grey80"),
    legend.text = element_text(face = "bold", size = 10),
    legend.title = element_blank(),
    plot.title = element_text(face = "bold", size = 18),
    panel.background = element_rect(fill = "transparent", colour = NA),
    panel.grid.major  = element_line(colour = "black", size = 0.5),
    panel.grid.minor  = element_line(colour = "grey70", size = 0.2),
    plot.background = element_rect(fill = "transparent", colour = NA),
    strip.background  = element_rect(fill = "transparent", colour = NA),
    strip.text = element_text(colour = "black", face = "bold", size = 24)
  )

## General Responsiveness ####

respall <-
  humanmacaque %>%
  filter(callall > 0) %>%
  mutate(child_behav_rate = callall / (seconds / 60)) %>% 
  mixed(
    r_eall ~ group * age_allweeks + child_behav_rate + (1 | id),
    offset =
      log(
        humanmacaque %>%
          filter(
            callall > 0
            ) %>%
          pull(callall)
        ),
    data = .,
    control=glmerControl(optimizer="bobyqa"),
    family = "poisson",
    method = "LRT")

general_responsiveness <-
  humanmacaque %>%
  filter(callall > 0) %>%
  mutate(
      fitted = fitted(respall$full_model)) %>%
    mutate(
      perc = (fitted * 100) / callall,
      days =
        age_allweeks *
        sd(humanmacaque %>% filter(group == "Humans") %>% pull(days), na.rm = T) +
        mean(humanmacaque %>% filter(group == "Humans") %>% pull(days), na.rm = T)) %>%
    select(group, perc, days) %>%
    ggplot(data = .,
           aes(x = days, y = perc, color = group)) +
    geom_point() +
    geom_smooth(
      method = lm,
      se = FALSE,
      size = 2,
      fullrange = TRUE
    ) +
    scale_y_continuous(name = "General Responsiveness to Infant Behaviours\n(as Estimated Percentage of all Infant Behaviours)") +
    scale_x_continuous(name = "Age in Days (Humans)",
                       sec.axis =
                         sec_axis(~ ((. - mean(humanmacaque %>% filter(group == "Humans") %>% pull(days), na.rm = T)) /
                                       sd(humanmacaque %>% filter(group == "Humans") %>% pull(days), na.rm = T)) *
                                    sd(humanmacaque %>% filter(group == "Rhesus Macaques") %>% pull(days), na.rm = T) +
                                    mean(humanmacaque %>% filter(group == "Rhesus Macaques") %>% pull(days), na.rm = T),
                                  name = "Age in Days (Rhesus Macaques)",
                                  breaks = seq(0, 15, 2))) +
  theme(
    axis.title = element_text(face = "bold", size = 16),
    axis.text = element_text(face = "bold", colour = "black", size = 16),
    axis.ticks = element_blank(),
    legend.key = element_rect(fill = "white", colour = "grey80"),
    legend.text = element_text(face = "bold", size = 10),
    legend.title = element_blank(),
    plot.title = element_text(face = "bold", size = 18),
    panel.background = element_rect(fill = "transparent", colour = NA),
    panel.grid.major  = element_line(colour = "black", size = 0.5),
    panel.grid.minor  = element_line(colour = "grey70", size = 0.2),
    plot.background = element_rect(fill = "transparent", colour = NA),
    strip.background  = element_rect(fill = "transparent", colour = NA),
    strip.text = element_text(colour = "black", face = "bold", size = 24)
  )

## Responsiveness to Social Behaviour ####
respsb <-
  humanmacaque %>%
  filter(
    !is.na(age_lastweeks),
    csb > 0) %>%
  mutate(social_behav_ratio = csb / callall) %>% 
  mixed(
    r_sb ~ group * age_lastweeks + social_behav_ratio + (1 | id),
    offset =
      log(
        humanmacaque %>%
          filter(
            !is.na(age_lastweeks),
            csb > 0) %>%
          pull(csb)
        ),
    data = .,
    family = "poisson",
    method = "LRT")

responses_to_social <-
  humanmacaque %>%
  filter(
    !is.na(age_lastweeks),
    csb > 0) %>%
  mutate(
      fitted = fitted(respsb$full_model)) %>%
    mutate(
      perc = (fitted * 100) / csb,
      days =
        age_lastweeks *
        sd(humanmacaque %>% filter(group == "Humans", !is.na(age_lastweeks)) %>% pull(days), na.rm = T) +
        mean(humanmacaque %>% filter(group == "Humans", !is.na(age_lastweeks)) %>% pull(days), na.rm = T)) %>%
    select(group, perc, days) %>%
    ggplot(data = .,
           aes(x = days, y = perc, color = group)) +
    geom_point() +
    geom_smooth(
      method = lm,
      se = FALSE,
      size = 2,
      fullrange = TRUE
    ) +
    scale_y_continuous(name = "Responsiveness to Infant Social Behaviours\n(as Estimated Percentage of Infant Social Behaviours)") +
    scale_x_continuous(name = "Age in Days (Humans)",
                       sec.axis =
                         sec_axis(~ ((. - mean(humanmacaque %>% filter(group == "Humans", !is.na(age_lastweeks)) %>% pull(days), na.rm = T)) /
                                       sd(humanmacaque %>% filter(group == "Humans", !is.na(age_lastweeks)) %>% pull(days), na.rm = T)) *
                                    sd(humanmacaque %>% filter(group == "Rhesus Macaques") %>% pull(days), na.rm = T) +
                                    mean(humanmacaque %>% filter(group == "Rhesus Macaques") %>% pull(days), na.rm = T),
                                  name = "Age in Days (Rhesus Macaques)",
                                  breaks = seq(0, 15, 2))) +
  theme(
    axis.title = element_text(face = "bold", size = 16),
    axis.text = element_text(face = "bold", colour = "black", size = 16),
    axis.ticks = element_blank(),
    legend.key = element_rect(fill = "white", colour = "grey80"),
    legend.text = element_text(face = "bold", size = 10),
    legend.title = element_blank(),
    plot.title = element_text(face = "bold", size = 18),
    panel.background = element_rect(fill = "transparent", colour = NA),
    panel.grid.major  = element_line(colour = "black", size = 0.5),
    panel.grid.minor  = element_line(colour = "grey70", size = 0.2),
    plot.background = element_rect(fill = "transparent", colour = NA),
    strip.background  = element_rect(fill = "transparent", colour = NA),
    strip.text = element_text(colour = "black", face = "bold", size = 24)
  )

## Mirroring to Social Behaviour ####
mirsb <-
  humanmacaque %>%
  filter(
    !is.na(age_lastweeks),
    csb > 0
    ) %>%
  mutate(
    social_behav_ratio = csb / callall,
  ) %>% 
  mixed(
    m_sb ~ group * age_lastweeks + social_behav_ratio + (1 | id),
    offset =
      log(
        humanmacaque %>%
          filter(
            !is.na(age_lastweeks),
            csb > 0
            ) %>%
          pull(csb)),
    data = .,
    family = "poisson",
    method = "LRT")

mirroring_to_social <-
  humanmacaque %>%
  filter(
    !is.na(age_lastweeks),
    csb > 0
  ) %>%
  mutate(
    fitted = fitted(mirsb$full_model)) %>%
  mutate(
    perc = (fitted * 100) / csb,
    days =
      age_lastweeks *
      sd(humanmacaque %>% filter(group == "Humans", !is.na(age_lastweeks)) %>% pull(days), na.rm = T) +
      mean(humanmacaque %>% filter(group == "Humans", !is.na(age_lastweeks)) %>% pull(days), na.rm = T)) %>%
    select(group, perc, days) %>%
    ggplot(data = .,
           aes(x = days, y = perc, color = group)) +
    geom_point() +
    geom_smooth(
      method = lm,
      se = FALSE,
      size = 2,
      fullrange = TRUE
    ) +
    scale_y_continuous(name = "Mirroring to Infant Social Behaviours\n(as Estimated Percentage of Infant Social Behaviours)") +
    scale_x_continuous(name = "Age in Days (Humans)",
                       sec.axis =
                         sec_axis(~ ((. - mean(humanmacaque %>% filter(group == "Humans", !is.na(age_lastweeks)) %>% pull(days), na.rm = T)) /
                                       sd(humanmacaque %>% filter(group == "Humans", !is.na(age_lastweeks)) %>% pull(days), na.rm = T)) *
                                    sd(humanmacaque %>% filter(group == "Rhesus Macaques") %>% pull(days), na.rm = T) +
                                    mean(humanmacaque %>% filter(group == "Rhesus Macaques") %>% pull(days), na.rm = T),
                                  name = "Age in Days (Rhesus Macaques)",
                                  breaks = seq(0, 15, 2))) +
  theme(
    axis.title = element_text(face = "bold", size = 16),
    axis.text = element_text(face = "bold", colour = "black", size = 16),
    axis.ticks = element_blank(),
    legend.key = element_rect(fill = "white", colour = "grey80"),
    legend.text = element_text(face = "bold", size = 10),
    legend.title = element_blank(),
    plot.title = element_text(face = "bold", size = 18),
    panel.background = element_rect(fill = "transparent", colour = NA),
    panel.grid.major  = element_line(colour = "black", size = 0.5),
    panel.grid.minor  = element_line(colour = "grey70", size = 0.2),
    plot.background = element_rect(fill = "transparent", colour = NA),
    strip.background  = element_rect(fill = "transparent", colour = NA),
    strip.text = element_text(colour = "black", face = "bold", size = 24)
  )

## Affirmative Marking to Social Behaviour ####
assb <-
  humanmacaque %>%
  filter(
    !is.na(age_lastweeks),
    csb > 0
    ) %>%
  mutate(social_behav_ratio = csb / callall) %>% 
  mixed(
    a_sb ~ group * age_lastweeks + social_behav_ratio + (1 | id),
    offset =
      log(
        humanmacaque %>%
          filter(
            !is.na(age_lastweeks),
            csb > 0
            ) %>%
          pull(csb)),
    data = .,
    family = "poisson",
    method = "LRT")

marking_to_social <-
  humanmacaque %>%
  filter(
    !is.na(age_lastweeks),
    csb > 0
  ) %>%
  mutate(
      fitted = fitted(assb$full_model)) %>%
    mutate(
      perc = (fitted * 100) / csb,
      days =
        age_lastweeks *
        sd(humanmacaque %>% filter(group == "Humans", !is.na(age_lastweeks)) %>% pull(days), na.rm = T) +
        mean(humanmacaque %>% filter(group == "Humans", !is.na(age_lastweeks)) %>% pull(days), na.rm = T)) %>%
    select(group, perc, days) %>%
    ggplot(data = .,
           aes(x = days, y = perc, color = group)) +
    geom_point() +
    geom_smooth(
      method = lm,
      se = FALSE,
      size = 2,
      fullrange = TRUE
    ) +
    scale_y_continuous(name = "Marking to Infant Social Behaviours\n(as Estimated Percentage of Infant Social Behaviours)") +
    scale_x_continuous(name = "Age in Days (Humans)",
                       sec.axis =
                         sec_axis(~ ((. - mean(humanmacaque %>% filter(group == "Humans", !is.na(age_lastweeks)) %>% pull(days), na.rm = T)) /
                                       sd(humanmacaque %>% filter(group == "Humans", !is.na(age_lastweeks)) %>% pull(days), na.rm = T)) *
                                    sd(humanmacaque %>% filter(group == "Rhesus Macaques") %>% pull(days), na.rm = T) +
                                    mean(humanmacaque %>% filter(group == "Rhesus Macaques") %>% pull(days), na.rm = T),
                                  name = "Age in Days (Rhesus Macaques)",
                                  breaks = seq(0, 15, 2))) +
    theme(
      axis.title = element_text(face = "bold", size = 16),
      axis.text = element_text(face = "bold", colour = "black", size = 16),
      axis.ticks = element_blank(),
      legend.key = element_rect(fill = "white", colour = "grey80"),
      legend.text = element_text(face = "bold", size = 10),
      legend.title = element_blank(),
      plot.title = element_text(face = "bold", size = 18),
      panel.background = element_rect(fill = "transparent", colour = NA),
      panel.grid.major  = element_line(colour = "black", size = 0.5),
      panel.grid.minor  = element_line(colour = "grey70", size = 0.2),
      plot.background = element_rect(fill = "transparent", colour = NA),
      strip.background  = element_rect(fill = "transparent", colour = NA),
      strip.text = element_text(colour = "black", face = "bold", size = 24)
    )

## Specific Mirroring Responses ####
pshmir <-
  longhumanmacaque %>%
  filter(group == "Humans") %>% 
  glmer(
    psbin ~ specmir + age_lastweeks + baserate + (1| id),
    data = .,
    family = "binomial",
    control = glmerControl(optimizer = "bobyqa")
    )

pshmir_mixed <-
  longhumanmacaque %>%
  filter(group == "Humans") %>% 
  mixed(
    psbin ~ specmir + age_lastweeks + baserate + (1| id),
    data = .,
    family = "binomial",
    control = glmerControl(optimizer = "bobyqa"),
    method = "LRT"
    )

pshmir_emmeans <-
  emmeans(
    pshmir,
    pairwise ~ specmir,
    type = "response",
    adjust = "fdr"
    )

vocmir <-
  longhumanmacaque %>%
  filter(
    group == "Humans",
    specmir != "Enriched Mir."
    ) %>%
  mutate(
    specmir = factor(specmir, levels = c("Direct Mir.", "Modified Mir."))
    ) %>%
  glmer(
    vocbin ~ specmir + age_lastweeks + baserate + (1| id),
    data = .,
    family = "binomial",
    control = glmerControl(optimizer = "bobyqa")
    )

vocmir_mixed <-
  longhumanmacaque %>%
  filter(
    group == "Humans",
    specmir != "Enriched Mir."
    ) %>%
  mutate(specmir = factor(specmir, levels = c("Direct Mir.", "Modified Mir."))) %>% 
  mixed(
    vocbin ~ specmir + age_lastweeks + baserate + (1| id),
    data = .,
    family = "binomial",
    control = glmerControl(optimizer = "bobyqa"),
    method = "LRT")
  
vocmir_emmeans <-
  emmeans(
    vocmir,
    pairwise ~ specmir,
    type = "response",
    adjust = "fdr")

comhmir <-
  longhumanmacaque %>%
  filter(group == "Humans") %>% 
  glmer(
    combin ~ specmir + age_lastweeks + baserateaffcom + (1| id),
    data = .,
    family = "binomial",
    control = glmerControl(optimizer = "bobyqa")
    )

comhmirmix <-
  longhumanmacaque %>%
  filter(group == "Humans") %>% 
  mixed(
    combin ~ specmir + age_lastweeks + baserateaffcom + (1| id),
    data = .,
    family = "binomial",
    control = glmerControl(optimizer = "bobyqa"),
    method = "LRT")

comhmir_emmeans <-
  emmeans(
    comhmir,
    pairwise ~ specmir,
    type = "response",
    adjust = "fdr")

smilemir <-
  longhumanmacaque %>%
  filter(group == "Humans") %>% 
  glmer(
    smilebin ~ specmir + age_lastweeks + baserate + (1| id),
    data = .,
    family = "binomial",
    control = glmerControl(optimizer = "bobyqa")
    )

smilemirmix <-
  longhumanmacaque %>%
  filter(group == "Humans") %>% 
  mixed(
    smilebin ~ specmir + age_lastweeks + baserate + (1| id),
    data = .,
    family = "binomial",
    control = glmerControl(optimizer = "bobyqa"),
    method = "LRT")

smilemir_emmeans <-
  emmeans(
    smilemir,
    pairwise ~ specmir,
    type = "response",
    adjust = "fdr")

psmmir <-
  longhumanmacaque %>%
  filter(group != "Humans") %>% 
  glmer(
    psbin ~ specmir + age_lastweeks + baserate + (1| id),
    data = .,
    family = "binomial",
    control = glmerControl(optimizer = "bobyqa")
    )

psmmirmix <-
  longhumanmacaque %>%
  filter(group != "Humans") %>% 
  mixed(
    psbin ~ specmir + age_lastweeks + baserate + (1| id), data = .,
    family = "binomial",
    control = glmerControl(optimizer = "bobyqa"),
    method = "LRT")

psmmir_emmeans <-
  emmeans(
    psmmir,
    pairwise ~ specmir,
    type = "response",
    adjust = "fdr")

lsmir <-
  longhumanmacaque %>%
  filter(group != "Humans") %>% 
  glmer(
    lsbin ~ specmir + age_lastweeks + baserate + (1| id),
    data = .,
    family = "binomial",
    control = glmerControl(optimizer = "bobyqa")
    )

lsmirmix <-
  longhumanmacaque %>%
  filter(group != "Humans") %>% 
  mixed(
    lsbin ~ specmir + age_lastweeks + baserate + (1| id),
    data = .,
    family = "binomial",
    control = glmerControl(optimizer = "bobyqa"),
    method = "LRT")

lsmir_emmeans <-
  emmeans(
    lsmir,
    pairwise ~ specmir,
    type = "response",
    adjust = "fdr")

affmir <-
  longhumanmacaque %>% 
  glmer(
    affbin ~ specmir * group + age_lastweeks + baserateaffcom + (1| id/age_lastweeks),
    data = .,
    family = "binomial",
    control = glmerControl(optimizer = "bobyqa")
    )

affmirmix <-
  longhumanmacaque %>% 
  mixed(
    affbin ~ specmir * group + age_lastweeks + baserateaffcom + (1| id/age_lastweeks),
    data = .,
    family = "binomial",
    control = glmerControl(optimizer = "bobyqa"),
    method = "LRT")

affmir_emmeans <-
  emmeans(
    affmir,
    pairwise ~ specmir,
    type = "response",
    adjust = "fdr")

commir <-
  longhumanmacaque %>% 
  glmer(
    combin ~ specmir * group + age_lastweeks + baserateaffcom + (1| id/age_lastweeks),
    data = .,
    family = "binomial",
    control = glmerControl(optimizer = "bobyqa")
    )

commirmix <-
  longhumanmacaque %>% 
  mixed(
    combin ~ specmir * group + age_lastweeks + baserateaffcom + (1| id/age_lastweeks),
    data = .,
    family = "binomial",
    control = glmerControl(optimizer = "bobyqa"),
    method = "LRT")

commir_emmeans <-
  emmeans(
    commir,
    pairwise ~ specmir,
    type = "response",
    adjust = "fdr")

specific_mirroring <-
ggarrange(
         {ggplot(
           cbind(
             chisquaremir2 %>%
               group_by(group,
                        chbehav,
                        specmir) %>%
               summarise(n = n()) %>%
               complete(specmir, fill = list(n = 0)) %>%
               filter(specmir %in% c(
                 "Direct Mir.", "Enriched Mir.", "Modified Mir."
               )) %>% bind_cols(order=c(1, 2, 3, 7, 8, 9, 4, 5, 6, 10, 11, 12, 13, 14, 15)) %>% arrange(order),
             tot = c(rep(
               data.frame(
                 chisquaremir2 %>%
                   filter(
                     group == "Humans",
                     specmir %in% c("Direct Mir.", "Enriched Mir.", "Modified Mir.")
                   ) %>%
                   group_by(specmir) %>%
                   summarise(n = n())
               )$n,
               3
             ), rep(
               data.frame(
                 chisquaremir2 %>%
                   filter(
                     group == "Rhesus Macaques",
                     specmir %in% c("Direct Mir.", "Enriched Mir.", "Modified Mir.")
                   ) %>%
                   group_by(specmir) %>%
                   summarise(n = n())
               )$n,
               2
             ))
           ) %>%
             ungroup() %>% 
             mutate(perc = (n * 100) / tot,
                    chbehav = factor(ifelse(chbehav == "Proto-Communicative", "Proto-Communicative\nMouth Gestures\n(P)", ifelse(chbehav == "Vocalisations", "Vocalisations\n(V)", as.character(chbehav))), levels = c("Proto-Communicative\nMouth Gestures\n(P)", "Vocalisations\n(V)", "Smiles", "Lip-Smacking"))) %>%
             select(group,
                    specmir,
                    chbehav,
                    perc) %>% filter(group == "Humans", chbehav != "Smiles"),
           aes(x = chbehav, y = perc, fill = specmir)
         ) +
             geom_bar(stat = "identity", position = "dodge", color = "black") +
             scale_fill_manual(
               values = c("chartreuse4", "yellow3", "orangered2")
             ) +
             ggtitle("Individual Human Communicative Behaviours") +
             theme(
               axis.title.x = element_blank(),
               axis.title.y = element_text(size = 16, face = "bold"),
               axis.text = element_text(
                 colour = "black",
                 size = 12,
                 face = "bold"
               ),
               axis.ticks = element_line(colour = "black"),
               legend.key = element_rect(fill = "white", colour = "grey80"),
               legend.title = element_blank(),
               plot.margin = margin(30, 5.5, 5.5, 5.5, "pt"),
               legend.background = element_rect(fill = "transparent", colour = NA),
               legend.text = element_text(face = "bold", size = 14),
               plot.title = element_text(hjust = 0.5, size = 14, face = "bold")
             ) +
             scale_y_continuous(limits = c(0, 100),
                                name = "Maternal Mirroring (%)")},
         
         {ggplot(cbind(
           chisquaremir2 %>%
             group_by(group,
                      chbehav,
                      specmir) %>%
             summarise(n = n()) %>%
             complete(specmir, fill = list(n = 0)) %>%
             filter(specmir %in% c(
               "Direct Mir.", "Enriched Mir.", "Modified Mir."
             )) %>% bind_cols(order=c(1, 2, 3, 7, 8, 9, 4, 5, 6, 10, 11, 12, 13, 14, 15)) %>% arrange(order),
           tot = c(rep(
             data.frame(
               chisquaremir2 %>%
                 filter(
                   group == "Humans",
                   specmir %in% c("Direct Mir.", "Enriched Mir.", "Modified Mir.")
                 ) %>%
                 group_by(specmir) %>%
                 summarise(n = n())
             )$n,
             3
           ), rep(
             data.frame(
               chisquaremir2 %>%
                 filter(
                   group == "Rhesus Macaques",
                   specmir %in% c("Direct Mir.", "Enriched Mir.", "Modified Mir.")
                 ) %>%
                 group_by(specmir) %>%
                 summarise(n = n())
             )$n,
             2
           ))
         ) %>%
           ungroup() %>%
           mutate(perc = (n * 100) / tot,
                  chbehav = factor(chbehav, levels = c("Proto-Communicative", "Vocalisations", "Smiles", "Lip-Smacking"))) %>%
           select(group,
                  specmir,
                  chbehav,
                  perc) %>% filter(group == "Humans") %>% 
           mutate(group = factor(ifelse(chbehav == "Smiles", "Affiliative\n(Smiles)", "Communicative\n(P+V)"), levels = c("Communicative\n(P+V)", "Affiliative\n(Smiles)")),
                  perc2 = c(65.5, 54.9, 10.4, 91.4, 54.9, 97.9, 8.62, 45.1, 2.08),
                  Number = c("P", "P", "P", "V", NA, "V", NA, NA, NA),
                  order = c(4, 5, 6, 1, 2, 3, 7, 8, 9)) %>% 
           arrange(order) %>% 
           select(group, specmir, perc2, Number),
         aes(x=group, y=perc2, fill=specmir, label = Number)) +
             geom_bar(stat="identity", position="dodge", color = "black") +
             geom_text(position=position_dodge(width=0.9), vjust=1.6, size = 4, fontface = 2) +
             scale_fill_manual(
               values = c("chartreuse4", "yellow3", "orangered2")
             ) +
             ggtitle("Humans") +
             theme(
               axis.title.x = element_blank(),
               axis.title.y = element_text(size = 16, face = "bold"),
               axis.text = element_text(
                 colour = "black",
                 size = 12,
                 face = "bold"
               ),
               axis.ticks = element_line(colour = "black"),
               legend.key = element_rect(fill = "white", colour = "grey80"),
               legend.title = element_blank(),
               plot.margin = margin(30, 5.5, 5.5, 5.5, "pt"),
               legend.position = "none",
               legend.background = element_rect(fill = "transparent", colour = NA),
               legend.text = element_text(face = "bold", size = 14),
               plot.title = element_text(hjust = 0.5, size = 14, face = "bold")
             ) +
             scale_y_continuous(name = "Maternal Mirroring (%)") +
             coord_cartesian(ylim=c(0,100))},
         
         {ggplot(cbind(
           chisquaremir2 %>%
             group_by(group,
                      chbehav,
                      specmir) %>%
             summarise(n = n()) %>%
             complete(specmir, fill = list(n = 0)) %>%
             filter(specmir %in% c(
               "Direct Mir.", "Enriched Mir.", "Modified Mir."
             )) %>% bind_cols(order=c(1, 2, 3, 7, 8, 9, 4, 5, 6, 10, 11, 12, 13, 14, 15)) %>% arrange(order),
           tot = c(rep(
             data.frame(
               chisquaremir2 %>%
                 filter(
                   group == "Humans",
                   specmir %in% c("Direct Mir.", "Enriched Mir.", "Modified Mir.")
                 ) %>%
                 group_by(specmir) %>%
                 summarise(n = n())
             )$n,
             3
           ), rep(
             data.frame(
               chisquaremir2 %>%
                 filter(
                   group == "Rhesus Macaques",
                   specmir %in% c("Direct Mir.", "Enriched Mir.", "Modified Mir.")
                 ) %>%
                 group_by(specmir) %>%
                 summarise(n = n())
             )$n,
             2
           ))
         ) %>%
           ungroup() %>%
           mutate(perc = (n * 100) / tot,
                  chbehav = factor(chbehav, levels = c("Proto-Communicative", "Vocalisations", "Smiles", "Lip-Smacking"))) %>%
           select(group,
                  specmir,
                  chbehav,
                  perc) %>% filter(group == "Rhesus Macaques") %>% 
           mutate(chbehav = factor(ifelse(chbehav == "Lip-Smacking", "Affiliative\n(Lip-Smacking)", "Communicative\n(P)"), levels = c("Communicative\n(P)", "Affiliative\n(Lip-Smacking)"))),
         aes(x = chbehav, y = perc, fill = specmir)
         ) +
             geom_bar(stat = "identity", position = "dodge", color = "black") +
             scale_fill_manual(
               values = c("chartreuse4", "yellow3", "orangered2")
             ) +
             ggtitle("Rhesus Macaques") +
             theme(
               axis.title.x = element_blank(),
               axis.title.y = element_text(size = 16, face = "bold"),
               axis.text = element_text(
                 colour = "black",
                 size = 12,
                 face = "bold"
               ),
               axis.ticks = element_line(colour = "black"),
               legend.key = element_rect(fill = "white", colour = "grey80"),
               legend.title = element_blank(),
               plot.margin = margin(30, 5.5, 5.5, 5.5, "pt"),
               legend.background = element_rect(fill = "transparent", colour = NA),
               legend.text = element_text(face = "bold", size = 14),
               plot.title = element_text(hjust = 0.5, size = 14, face = "bold")
             ) +
             scale_y_continuous(limits = c(0, 100),
                                name = "Maternal Mirroring (%)")},
         ncol = 1, nrow = 3, common.legend = TRUE, labels = c("a)", "b)", "c)"), font.label = list(size = 20, face = "bold"))