Plot theme & colors

theme_set(obj$theme)
obj$colors_sex <- c('#eeb054', '#74ccd4')
obj$colors_rank <- c('#F26A21', '#BBC0A4', '#BADAEE', '#4FB2DD')

# Female Control, Male Control, Female CSM, Male CMS
obj$four_color = c('#eeb054', '#74ccd4', '#cd6632', '#317f8f')

Overview

We provide the combined outcomes of all tests in a single object that is loaded here as obj. Below is the structure of this object to help anyone access the data. The obj is available on the repository under data/obj.Rds.

`master`

This contains the master list of Social Box data for each individual on each tracking day. Here are some of the variables in this data table:

obj$master %>% select(RealMouseNumber, GroupNumber, Color, Sex, Day, Suspect)

`bodyweights`

This table contains the bodyweights & time of assessment for all individuals:

obj$bodyweights %>% as_tibble()

`organs`

This table contains several physiological parameters collected throughout the experiment, including adrenal weights, thymus weights, fur scores, and corticosterone values. The bodyweights at the time of sacrifice are also included and used for adjusting bodyweight-dependent variables.

obj$organs %>% as_tibble()

`tests`

This list contains separate tables containing the results of the classical tests that were run here.

obj$tests$EPM - Elevated Plus Maze results (& EPM in 1 min intervals for data exploration)
obj$tests$Nest - Nest building
obj$tests$OFT - Open Field Test
obj$tests$SP - Sucrose preference
obj$tests$TST - Tail suspension test
obj$tests$Splash - Splash test

`tests_combined`

obj$tests_combined_adj - combined table with all tests included after each readout has been adjusted for Batch effects and centered on the mean of the Female Control condition

Outlier removal

CORT - a single outlier with a value of 453.4 was removed
EPM - one failed video/tracking test was removed
Nest building - the “percent intact” value for one animal could not be read back from the document

obj$master <- obj$master %>% filter(!Suspect)
obj$profiles <- obj$profiles %>% filter(!Suspect)

Figure 1 - Baseline Dominance

Baseline Cumulative DS

Cumulative David’s Score based on chases

res <- lapply(1:4, function(d) {
  p <- explore_hierarchies(obj$master %>% filter(Day %in% 1:d), 
                           idVars = c('Sex', 'Condition'), 
                           verbose = F, DS_only = T, daily = F) %>% 
    select(GroupNumber, Sex, DS) %>% 
    unnest(cols = DS) %>% 
    mutate(MouseID = as.numeric(MouseID)) %>% 
    rename(DS_cum = DS)
  p$Day <- d
  return(p)
})
res <- do.call('rbind', res)
res <- left_join(res, obj$mice[,c('MouseID', 'GroupNumber', 'Sex', 'RealMouseNumber')])

# Rank day 4
p2 <- res %>% 
  filter(Day == 4) %>% 
  group_by(GroupNumber, Sex) %>% 
  mutate(
    Rank = rank(-DS_cum), 
    Rank = plyr::mapvalues(Rank, 1:4, c('alpha', 'beta', 'gamma', 'delta')) ,
    Rank = factor(Rank, levels = c('alpha', 'beta', 'gamma', 'delta'))
  )
res <- left_join(res, p2[,c('RealMouseNumber', 'Rank')])

r2 <- res %>% group_by(Sex, Day, Rank) %>% 
  dplyr::summarize(
    M = mean(DS_cum, na.rm = T),
    se = sd(DS_cum, na.rm = T)/sqrt(n()),
  ) 
f1_1 <- ggplot(r2, aes(x = as.factor(Day), 
                       y = M, 
                       fill = Rank, 
                       color = Rank)) + 
  geom_line(data = res, 
            aes(group = RealMouseNumber, y = DS_cum), 
            alpha = 0.4, size = 0.4) + 
  geom_point(data = res, 
             aes(group = RealMouseNumber, y = DS_cum), 
             alpha = 0.4, size = 0.4) + 
  geom_errorbar(aes(ymin = M-se, ymax = M+se), width = 0) +
  geom_line(aes(group = Rank), size = 0.8) + 
  labs(y = "Baseline David's Score", x = 'Day') +
  geom_point(size = 2.5, pch = 21, color = 'grey20') + 
  scale_fill_manual(values = obj$colors_rank, name = 'Final rank', 
                    labels = c('\u03B1', '\u03B2', '\u03B3', '\u03B4')) + 
  scale_color_manual(values = obj$colors_rank, name = 'Final rank', 
                     labels = c('\u03B1', '\u03B2', '\u03B3', '\u03B4')) + 
  scale_x_discrete(expand = c(0.1, 0.1)) +
  facet_grid(~Sex)  + theme(legend.position = 'top')
f1_1

Correlations with SB behaviors

Cumulative David’s Scores from days 1:4
SB behaviors summarized as means of days 2:4

h_cum <- explore_hierarchies(obj$master %>% filter(Day %in% 1:4), 
                             idVars = c('Sex', 'Condition'), 
                             DS_only = T, verbose = F, daily = F) %>%
  select(GroupNumber, Sex, Condition, DS) %>% 
  unnest(cols = DS) %>% 
  mutate(MouseID = as.numeric(as.character(MouseID))) %>% 
  left_join(obj$mice[,c('RealMouseNumber', 'Sex', 'GroupNumber', 'MouseID')])

h_cum <- h_cum %>% 
  group_by(GroupNumber, Sex) %>% 
  mutate(
    Rank = rank(-DS), 
    Rank = plyr::mapvalues(Rank, 1:4, c('alpha', 'beta', 'gamma', 'delta')) ,
    Rank = factor(Rank, levels = c('alpha', 'beta', 'gamma', 'delta'))
  )

h_daily <- explore_hierarchies(obj$master %>% filter(Day %in% 1:4), 
                               idVars = c('Sex', 'Condition'), 
                               DS_only = T, verbose = F, daily = T) %>%
  select(GroupNumber, Sex, Condition, Day, DS) %>% 
  unnest(cols = DS) %>% 
  mutate(MouseID = as.numeric(as.character(MouseID))) %>% 
  left_join(obj$mice[,c('RealMouseNumber', 'Sex', 'GroupNumber', 'MouseID')])


p <- subset(obj$profiles, !Suspect & Day %in% 2:4) %>% 
  select(RealMouseNumber, all_props) %>% 
  group_by(RealMouseNumber) %>% 
  summarise_at(.vars = all_props, mean, na.rm = T)


p <- subset(obj$profiles, !Suspect & Day %in% 2:4) %>% 
  select(RealMouseNumber, Batch, all_props) %>% 
  pivot_longer(cols = all_of(all_props), names_to = "prop") %>% 
  filter(!is.na(value)) %>%
  group_by(prop, Batch) %>%
  mutate(
    value = RNOmni::rankNorm(value, k = 3/8)
  ) %>% 
  group_by(RealMouseNumber, prop) %>% 
  summarize(value = mean(value)) %>% 
  pivot_wider(names_from = 'prop')

d <- left_join(h_cum, p)

idVars <- c('RealMouseNumber', 'GroupNumber', 'Sex', 'DS')
props <- c('DistanceOutside', 'GridEntropy6x6', 'FractionOfTimeOutside')

p <- d %>% 
  ungroup() %>% 
  select(idVars, props) %>% 
  pivot_longer(cols = -idVars, names_to = "prop", values_to = "value") %>% 
  mutate(prop = plyr::mapvalues(prop, all_props, all_props_pretty))

f1_2 <- ggplot(p, aes(x = DS, y = value, fill = Sex)) +
  geom_smooth(method = 'lm', se = F, lty = 2, aes(col = Sex)) +
  geom_point(pch = 21, size = 2) +
  labs(x = "Baseline David's Score", y = '') +
  facet_grid(prop~Sex, scales = 'free_y', switch = "y") +
  scale_fill_manual(values = obj$colors_sex) +
  scale_color_manual(values = obj$colors_sex) +
  theme(legend.position = 'none', strip.placement = "outside")
f1_2

Stats

props <- all_props
p <- d %>% 
  ungroup() %>% 
  select(idVars, props) %>% 
  pivot_longer(cols = -idVars, names_to = "prop", values_to = "value") %>% 
  group_by(prop, Sex) %>% 
  mutate(N = n()) %>%
  group_by(prop, Sex, N) %>% 
  nest() %>% 
  mutate(
    M = map(data, function(d) {
      tidy(cor.test(d$value, d$DS, method = 'spearman'))
    })
  ) %>% 
  select(-data) %>% unnest() %>% 
  ungroup() %>% mutate(
    prop = fct_reorder(prop, abs(estimate))
  ) %>% 
  group_by(Sex) %>% 
  mutate(qval = p.adjust(p.value, method = 'BH')) %>% 
  arrange(p.value) %>% 
  mutate(sig_qval = qval <= 0.05)

Social box features with significant q-values:

knitr::kable(table(p$prop[p$sig_qval], p$Sex[p$sig_qval]))

	Female	Male
FractionOfTimeInFeederOutside	0	0
ProximateVsDistantWater	0	0
FractionOfTimeInTheOpenOutside	0	0
FractionOfTimeInSmallNestOutside	0	0
GridMI6x6	0	0
FoodOrWaterPerTimeOutside	0	0
DistanceFromWallsInOpen	0	0
MeanContactDuration	0	0
AngularVelocity	0	0
TangentialVelocity	0	0
EntropyOutside	0	0
MedianContactDuration	0	0
DistanceFromNest	0	0
ForagingCorrelation	0	0
HighPlacePerTimeOutside	0	0
FractionOfTimeInContactOutside	0	0
ContactRateOutside	0	0
FractionOfTimeOnRampOutside	0	0
MeanSpeedOutside	0	0
FractionOfTimeInLabyrinthOutside	0	0
ContactRate	0	0
ProximateVsDistantFood	0	0
ApproachRateOutside	0	0
MedianSpeedOutside	0	0
BeingApproachedRate	0	0
VisitsOutsideRate	0	0
FractionOfTimeAtHighPlace	0	0
DiffBetweenApproachesAndChases	0	0
FractionOfTimeNearFoodOrWater	0	1
FractionOfTimeInWaterOutside	0	0
FractionOfTimeAloneOutside	0	1
ApproachRate	0	1
Entropy	0	1
GridEntropy6x6	0	1
NumberOfApproachs	0	1
FractionOfTimeOutside	0	1
BeingApproachedRateOutside	0	1
NumberOfApproachesPerMiceOut	0	1
DistanceOutside	0	1
FractionOfBeingApproachedPerContact	0	1
FractionOfApproachesPerContact	0	1
FractionOfApproachEscapeBehaviorPerAggression	1	1
AggressiveChaseRateOutside	1	1
AggressiveChaseRate	0	1
AggressiveEscapeRate	1	1
NAEscapeRate	1	1
FractionOfChasesPerContact	0	1
NAChaseRateOutside	1	1
FollowRateOutside	1	1
BeingFollowedRate	1	1
NAChaseRate	1	1
FollowRate	1	1
AggressiveEscapeRateOutside	1	1
NAEscapeRateOutside	1	1
FractionOfEscapesPerContact	1	1
BeingFollowedRateOutside	1	1
FractionOfNAEscapesPerContact	1	1
FractionOfBeingFollowedPerContact	1	1
FractionOfFollowsPerContact	1	1
FractionOfNAChasesPerContact	1	1

Distance Outside:

p[p$prop == 'DistanceOutside',]

Fraction of Time Outside:

p[p$prop == 'FractionOfTimeOutside',]

Grid Entropy [6x6]

p[p$prop == 'GridEntropy6x6',]

Sum of daily changes

# Generate a distribution of expected rank changes
set.seed(42)
sampleGroupRanks <- function (i) {sample(1:i, size = i, replace = F)}

# Fake data
nMice <- length(unique(obj$master$RealMouseNumber))
nDays <- 4
nIter <- 10

ref <- lapply(1:nIter, function(i) {
  dat <- data.frame(
    Day = rep(1:nDays, nMice),
    MouseNumber = rep(1:nMice, each = nDays), 
    GroupNumber = rep(1:(nMice/4), each = 4*nDays)
  )
  dat <- dat %>% 
    group_by(GroupNumber, Day) %>% 
    nest() %>%
    mutate(
      Rank = map(data, function(d) sampleGroupRanks(nrow(d)))
    ) %>% 
    unnest(cols = c(data, Rank))
  
  d <- dat %>% 
    group_by(GroupNumber, MouseNumber) %>% 
    arrange(GroupNumber, MouseNumber, Day) %>% 
    mutate(
      lagRank = lag(Rank, 1),
      RankChange = abs(Rank - lagRank)
    ) %>% filter(!is.na(lagRank)) 
  ref <- d %>% 
    group_by(GroupNumber, MouseNumber, RankChange) %>%
    tally() %>%
    group_by(GroupNumber, MouseNumber) %>% 
    mutate(
      s = sum(n),
      SumRankChangesFraction = n / sum(n) 
    ) %>% 
    group_by(RankChange) %>% 
    tally() %>% 
    mutate(
      SumRankChangesFraction = n / sum(n)
    )
  ref$iter = i
  return(ref)
})

ref <- do.call('rbind', ref)
ref <- ref %>% 
  group_by(RankChange) %>% 
  dplyr::summarise(
    M = mean(SumRankChangesFraction, na.rm = T), 
    Qt = qt(0.975,n()-1) * sd(SumRankChangesFraction)/sqrt(n()),
    CI_95 = M + Qt,
    CI_05 = M - Qt
  )

# Real ranks throughout
h_daily <- explore_hierarchies(subset(obj$master, Day %in% 1:4), 
                               idVars = c('Sex', 'Condition'), 
                               DS_only = T,
                               verbose = F, 
                               daily = T) %>% 
  select(GroupNumber, Sex, Condition, Day, DS) %>% 
  unnest(cols = DS) %>% 
  mutate(MouseID = as.numeric(MouseID)) %>% 
  group_by(GroupNumber, Sex, Condition, Day) %>% 
  mutate(DS_rank = rank(-DS)) %>% 
  left_join(obj$mice[,c('RealMouseNumber', 'GroupNumber', 'Sex', 'MouseID')])

p <- left_join(h_daily, h_cum[,c('RealMouseNumber', 'Rank')])

p <- p %>% 
  group_by(GroupNumber, Sex, Condition, MouseID, Rank) %>% 
  arrange(GroupNumber, Sex, MouseID, Day) %>% 
  mutate(
    lagRank = lag(DS_rank),
    RankChange = abs(DS_rank - lagRank)
  ) %>% 
  filter(!is.na(RankChange))

pp <- p %>% 
  group_by(Sex, Rank, RankChange) %>%
  tally() %>%
  group_by(Sex, Rank) %>%
  mutate(
    s = sum(n),
    SumRankChangesFraction = n / sum(n, na.rm = T) 
  ) 

ggplot(pp, aes(x = RankChange, y = SumRankChangesFraction)) +
  geom_ribbon(data = ref, aes(ymax = CI_95, ymin = CI_05, y = M), 
              fill = 'grey', alpha = 0.5) +
  geom_line(data = ref, color = 'grey', lty = 1, aes(y = M)) +
  geom_line(aes(col = Rank)) +
  geom_point(aes(fill = Rank), pch = 21, color = 'grey20', size = 2.5) +
  facet_grid(.~Sex) + 
  labs(y = 'Probability') +
  scale_color_manual(values = obj$colors_rank) +
  scale_fill_manual(values = obj$colors_rank) + 
  theme(legend.position = 'top')

x <- left_join(pp, ref)
b <- seq(-2, 2, 1)
l <- c(-4, -2, 0, 2, 4)
ggplot(x, aes(x = as.factor(RankChange), y = log2(SumRankChangesFraction / M))) +
  geom_hline(yintercept = 0, lty = 2, color = 'grey') +
  geom_line(aes(col = Rank, group = Rank), size = 0.5, 
            position = position_dodge(0.4)) +
  geom_point(aes(fill = Rank, col = Rank), pch = 21,
             color = 'grey20',
             size = 2.5, 
             position = position_dodge(0.4)) +
  labs(y = 'Fold change in probability\n(true / expected)', x = 'Absolute rank change') +
  facet_grid(.~Sex) + 
  scale_y_continuous(breaks = b, labels = l) +
  scale_fill_manual(values = obj$colors_rank) + 
  scale_color_manual(values = obj$colors_rank) + 
  theme(legend.position = 'none')

Change/no change analysis

p <- left_join(h_daily, h_cum[,c('RealMouseNumber', 'Rank')]) %>% 
  group_by(GroupNumber, Sex, Condition, MouseID, Rank) %>% 
  arrange(GroupNumber, Sex, MouseID, Day) %>% 
  mutate(
    lagRank = lag(DS_rank),
    RankChange = abs(DS_rank - lagRank)
  ) %>% 
  filter(!is.na(RankChange))

pp <- p %>% 
  dplyr::mutate(
    NoChange = RankChange == 0
  )

ppp <- pp %>% 
  group_by(Rank, Sex) %>%
  dplyr::summarise(
    N = n(), 
    S = sum(NoChange)
  ) %>%
  dplyr::mutate(
    P = S / N
  )

ppp$sig <- NA
ppp$sig[ppp$Rank == 'alpha'] <- '***'
ppp$sig[ppp$Rank == 'delta' & ppp$Sex == 'Male'] <- '***'
ppp$sig[ppp$Rank %in% c('beta', 'gamma') & ppp$Sex == 'Male'] <- '*'

ppp$Rank <- plyr::mapvalues(ppp$Rank, c('alpha', 'beta', 'gamma', 'delta'), 
                            c( '\u03B1', '\u03B2', '\u03B3', '\u03B4'))
f1_3 <- ggplot(ppp, aes(x = Rank, y = P / .25, fill = Rank)) +
  geom_bar(stat = 'identity', width = 0.7, alpha = 1, aes(color = Rank)) +
  geom_hline(yintercept = 1, lty = 1, color = 'grey30') +
  geom_text(aes(label = Rank, y = 1), nudge_y = -0.08, parse = F, fontface = 'italic', color = 'grey30') +
  geom_text(aes(label = N), color = 'grey20', nudge_y = 0.08, parse = F, size = 3, fontface = 'italic', color = 'grey30') +
  geom_text(aes(label = sig), color = 'grey20', nudge_y = 0.15, parse = F, size = 4, fontface = 'bold', color = 'grey30') +
  facet_grid(.~Sex) + 
  labs(y = 'Rank maintenance odds', x = '') +
  scale_color_manual(values = obj$colors_rank) +
  scale_fill_manual(values = obj$colors_rank) +
  scale_y_continuous(trans = 'log2', breaks = c(0.66, 1, 1.5, 2, 2.5, 3)
                     # , limits = c(0.66, 3)
  ) +
  theme(legend.position = 'none', 
        axis.text.x =  element_blank())

f1_3

Stats:

res <- lapply(1:nrow(ppp), function(i) {
  d = ppp[i,]
  tidy(binom.test(d$S, n = d$N, p = .25, alternative = 'greater'))
})
res <- do.call('rbind', res)
ppp <- cbind(as.data.frame(ppp), as.data.frame(res)) %>% arrange(Sex, Rank)
as.data.frame(ppp)

DS baseline v stress

baseline <- h_cum
stress <- explore_hierarchies(subset(obj$master, Day == 5), 
                              idVars = c('Sex', 'Condition'), 
                              verbose = F, 
                              DS_only = T,
                              daily = T) %>%
  select(GroupNumber, Sex, DS) %>% 
  unnest(cols = DS) %>% 
  rename(DS_stress = DS) %>% 
  mutate(MouseID = as.numeric(MouseID))
h <- left_join(baseline, stress)

Baseline v. Stress DS

f1_4 <- ggplot(h, aes(x = DS, y = DS_stress, fill = Sex)) +
  geom_smooth(method = 'lm', se = F, lty = 2, aes(col = Sex)) +
  geom_point(pch = 21, size = 2) +
  labs(y = "David's Score\nfollowing acute stress", x = "Baseline David's Score") +
  facet_grid(.~Sex) +
  scale_fill_manual(values = obj$colors_sex) +
  scale_color_manual(values = obj$colors_sex) +
  theme(legend.position = 'none')
f1_4

Stats

h %>% group_by(Sex) %>% 
  filter(!(is.na(DS) | is.na(DS_stress))) %>% 
  dplyr::summarize(n = n())

Males:

hh <- subset(h, Sex == 'Male')
data.frame(tidy(cor.test(hh$DS, hh$DS_stress)), N = nrow(hh))

Females:

hh <- subset(h, Sex == 'Female')
data.frame(tidy(cor.test(hh$DS, hh$DS_stress)), N = nrow(hh))

Chases after stress

baseline <- explore_hierarchies(subset(obj$master, 
                                       Day %in% 2:4), 
                                idVars = c('Sex', 'Condition'), 
                                DS_only = T,
                                verbose = F, 
                                daily = T) %>% 
  select(GroupNumber, Day, Sex, Condition, data, mat) %>% 
  unnest(cols = data) %>% 
  group_by(GroupNumber, Sex, Condition, MouseID, Day) %>% 
  summarize(wins = sum(wins, na.rm = T)) %>% # Collapse over OtherID
  group_by(GroupNumber, Sex, Condition, MouseID) %>% 
  summarize(wins = mean(wins, na.rm = T)) # Collapse over Day

stress <- explore_hierarchies(subset(obj$master, Day == 5), 
                              idVars = c('Sex', 'Condition'), 
                              DS_only = T,
                              verbose = F, 
                              daily = T) %>%
  select(GroupNumber, Sex, data) %>% 
  unnest(cols = data) %>% 
  rename(wins_stress = wins) %>% 
  group_by(GroupNumber, Sex, MouseID) %>% 
  summarize(wins_stress = sum(wins_stress, na.rm = T))  # Collapse over OtherID

h2 <- left_join(baseline, stress)

Chases from baseline to stress

Mean number of chases daily: baseline (days 2:4) vs. stress

pp <- h2 %>% 
  group_by(GroupNumber, MouseID, Sex, Condition) %>% 
  summarise_at(.vars = c('wins', 'wins_stress'), mean, na.rm = T) %>% 
  gather(Stress, Chases, -GroupNumber, -Sex, -Condition, -MouseID) %>% 
  mutate(Stress = plyr::mapvalues(Stress, c('wins', 'wins_stress'), c('Baseline', 'Acute Stress')), 
         Stress = factor(Stress, levels = c('Baseline', 'Acute Stress')))

Md <- pp %>% group_by(Sex, Stress) %>% 
  dplyr::summarise(
    Md = median(Chases, na.rm = T),
    Mn = mean(Chases, na.rm = T),
    SE = sd(Chases, na.rm = T) / sqrt(n()),
    MAD = stats::mad(Chases, na.rm = T)
  )

f1_5 <- ggplot(pp,
               aes(x = Chases, fill = Sex, color = Sex)) +
  geom_density(alpha = 0.7) +
  geom_rug(aes(color = Sex), size = 1, show.legend = F) +
  geom_errorbarh(data = Md[Md$Sex == 'Male',], inherit.aes = F, 
                 aes(xmin = Md, xmax = Md + MAD, y = -0.1, col = Sex), height = 0.1, show.legend = F) +
  geom_point(data = Md[Md$Sex == 'Male',], 
             pch = 21, size = 2.5, aes(x = Md, y = -0.1), color = 'grey20', show.legend = F) +
  geom_errorbarh(data = Md[Md$Sex == 'Female',], inherit.aes = F, 
                 aes(xmin = Md, xmax = Md + MAD, y = -0.2, col = Sex), height = 0.1, show.legend = F) +
  geom_point(data = Md[Md$Sex == 'Female',], 
             pch = 21, size = 2.5, aes(x = Md, y = -0.2), color = 'grey20', show.legend = F) +
  scale_fill_manual(values = obj$colors_sex) +
  scale_color_manual(values = obj$colors_sex) +
  facet_grid(Stress~.) + 
  scale_x_log10() + labs(y = '') +
  theme(legend.position = c(0.1, 0.9), 
        legend.title = element_blank(),
        panel.grid.major.x = element_line(color = 'grey', size = 0.2, linetype = 2),
        # panel.grid.minor.x = element_line(color = 'grey', size = 0.1),
        panel.border = element_blank(), 
        axis.text.y = element_blank())
f1_5

Stats

Repeated-measures ANOVA on log-transformed chases

s <- h2 %>% gather(Stage, Chases, -GroupNumber, -Sex, -Condition, -MouseID)
s$Stage <- plyr::mapvalues(s$Stage, c('wins', 'wins_stress'), c('Baseline', 'Stress'))
s <- left_join(s, obj$mice[,c('GroupNumber', 'Sex', 'MouseID', 'RealMouseNumber')])
summary(aov(log(Chases + 1) ~ Sex * Stage + Error(as.factor(RealMouseNumber)), data = s))

## 
## Error: as.factor(RealMouseNumber)
##           Df Sum Sq Mean Sq F value Pr(>F)
## Sex        1   0.02  0.0177   0.011  0.916
## Residuals 86 136.06  1.5821               
## 
## Error: Within
##           Df Sum Sq Mean Sq F value      Pr(>F)    
## Stage      1  5.963   5.963  29.040 0.000000611 ***
## Sex:Stage  1  0.216   0.216   1.053       0.308    
## Residuals 86 17.660   0.205                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Figure panel

col1 <- cowplot::plot_grid(plotlist = list(f1_1, f1_2), 
                           labels = c('a', 'b'), 
                           label_fontface = 'plain', 
                           ncol = 1, align = 'v', axis = 'lr', rel_heights = c(1, 2))
col2 <- cowplot::plot_grid(plotlist = list(f1_3, f1_4, f1_5), 
                           labels = c('c', 'd', 'e'), 
                           label_fontface = 'plain', 
                           ncol = 1, align = 'v', axis = 'lr', rel_heights = c(1, 0.8, 1.2))
cowplot::plot_grid(plotlist = list(col1, col2), ncol = 2, align = 'h')

if(exportFigures) {
  ggsave('../manuscript/figures/1_Panel.pdf', width = 9.5, height = 9.3, dpi = 1200)
  ggsave('../manuscript/figures/1_Panel.png', width = 9.5, height = 9.3, dpi = 1200)
}

Figure 2 - CMS outcomes

Physiological Outcomes

dat <- obj$tests_combined_adj_long %>% 
  filter(!has_NA) %>%
  mutate(
    Sex = factor(Sex, levels = c('Female', 'Male')), 
    SexCondition = paste0(Sex, '_', Condition), 
    SexCondition = factor(SexCondition, 
                          levels = c('Female_Control', 'Male_Control', 'Female_CMS', 'Male_CMS')
    )
  )

Bodyweight

Stats:

Summary:

var_int <- 'ORG_BodyweightChange'
d <- dat %>% filter(prop == var_int)
d %>% group_by(Sex, Condition) %>% 
  dplyr::summarize(
    n = n(), 
    n_groups = length(unique(GroupNumber)),
    M = mean(value, na.rm = T), 
    sd = sd(value, na.rm = T)
  )

Shapiro-Wilk test

d %>% group_by(Sex, Condition) %>% nest() %>% 
  mutate(SW = map(data, function(d) shapiro.test(d$value)), 
         tSW = map(SW, broom::tidy)) %>% 
  select(-data, -SW) %>% 
  unnest(cols = c(tSW))

Normality is not violated in any group

Levene’s Test

knitr::kable(car::leveneTest(value ~ Sex * Condition, d))

	Df	F value	Pr(>F)
group	3	1.23259	0.3032162
	82	NA	NA

Groupwise variances are not heteroscedastic

Males

knitr::kable(car::leveneTest(value ~ Condition, d[d$Sex == 'Male',]))

	Df	F value	Pr(>F)
group	1	1.736373	0.1957002
	37	NA	NA

Females

knitr::kable(car::leveneTest(value ~ Condition, d[d$Sex == 'Female',]))

	Df	F value	Pr(>F)
group	1	0.4261087	0.5172257
	45	NA	NA

Test

tidy(aov(value ~ Condition * Sex, data = d)) %>% as.data.frame()

summary(aov(value ~ Condition * Sex, data = d))

##               Df Sum Sq Mean Sq F value             Pr(>F)    
## Condition      1   69.4    69.4   7.394            0.00798 ** 
## Sex            1  755.3   755.3  80.453 0.0000000000000827 ***
## Condition:Sex  1   21.7    21.7   2.315            0.13201    
## Residuals     82  769.8     9.4                               
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise

dd <- d %>% 
  group_by(Sex) %>% 
  nest() %>% 
  dplyr::mutate(
    t = map(data, function(d) {
      # x = tidy(t.test(d$value ~ d$Condition))
      t.test(d$value ~ d$Condition)
      # x = tidy(kruskal.test(value ~ Condition, d))
    }
    )
  ) %>% 
  dplyr::select(-data)

Females:

dd$t[dd$Sex == "Female"][[1]]

## 
##  Welch Two Sample t-test
## 
## data:  d$value by d$Condition
## t = 3.9447, df = 43.784, p-value = 0.0002849
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  1.557491 4.812360
## sample estimates:
##      mean in group Control          mean in group CMS 
##  0.00000000000000004108419 -3.18492564471728734076805

Males:

dd$t[dd$Sex == "Male"][[1]]

## 
##  Welch Two Sample t-test
## 
## data:  d$value by d$Condition
## t = 1.1064, df = 36.937, p-value = 0.2757
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.9536554  3.2477764
## sample estimates:
## mean in group Control     mean in group CMS 
##              4.841554              3.694494

dd <- data.frame(Sex = c('Female', 'Male'), Sig = c('***', ' '), y = c(10, 10))
f2_1 <- ggplot(d, aes(x = Condition, y = value, fill = SexCondition)) +
  geom_hline(yintercept = 0, color = 'grey', lty = 2) +
  geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, aes(color = SexCondition)) +
  geom_point(pch = 21, size = 2, 
             position = position_jitterdodge(jitter.width = 0, 
                                             jitter.height = 0, dodge.width = 0.5)) + 
  scale_fill_manual(values = obj$four_color) + 
  scale_color_manual(values = obj$four_color) + 
  # scale_y_continuous(breaks = seq(-30, 30, 10)) +
  labs(x = '', y = 'Bodyweight change') +
  facet_grid(.~Sex) + theme(legend.position = 'none') + 
  geom_text(data = dd, inherit.aes = F, aes(x = 1.5, y = y, label = Sig), 
            fontface = 'bold', color = 'grey30')
f2_1

Coat state

d <- dat %>% 
  group_by(Sex, Condition, SexCondition) %>% 
  filter(prop == 'ORG_FurScore_Cumulative')

dd <- d %>%
  summarize(
    Md = median(value, na.rm = T),
    MAD = stats::mad(value, na.rm = T)
  )

sig_sum <- data.frame(Sex = c('Female', 'Male'), Sig = c('***', '***'), y = c(4, 4))

f2_2 <- ggplot(d, aes(x = Condition, y = value, fill = SexCondition)) +
  geom_dotplot(binaxis='y', stackdir='center', dotsize = 1, position = position_dodge(0.7)) +
  scale_fill_manual(values = obj$four_color) + 
  scale_color_manual(values = obj$four_color) + 
  labs(x = '', y = 'Coat state') +
  facet_grid(.~Sex) + theme(legend.position = 'none') +
  geom_text(data = sig_sum, inherit.aes = F, aes(x = 1.5, y = y, label = Sig), fontface = 'bold', color = 'grey30')
f2_2

Stats:

Summary:

d %>% group_by(Sex, Condition) %>% 
  dplyr::summarize(
    n = n(), 
    n_groups = length(unique(GroupNumber)),
    M = mean(value, na.rm = T), 
    sd = sd(value, na.rm = T)
  )

Shapiro-Wilk test

d %>% group_by(Sex, Condition) %>% 
  nest() %>% 
  mutate(SW = map(data, function(d) shapiro.test(d$value)), 
         tSW = map(SW, broom::tidy)) %>% select(-data, -SW) %>% 
  unnest(cols = c(tSW))

Normality is violated in both male groups as well as female controls

Levene’s Test

knitr::kable(car::leveneTest(value ~ Sex * Condition, d))

	Df	F value	Pr(>F)
group	3	1.346189	0.2651497
	82	NA	NA

Groupwise variances are not heteroscedastic

Transform data:

d <- d %>% 
  ungroup() %>% 
  dplyr::mutate(value_trans = rank(value)) 

ggplot(d, aes(x = Condition, y = value_trans, fill = SexCondition)) +
  geom_dotplot(binaxis='y', stackdir='center', dotsize = 1, position = position_dodge(0.7)) +
  scale_fill_manual(values = obj$four_color) + 
  scale_color_manual(values = obj$four_color) + 
  labs(x = '', y = 'Coat state') +
  facet_grid(.~Sex) + theme(legend.position = 'none')

d %>% group_by(Sex, Condition) %>% nest() %>% 
  mutate(SW = map(data, function(d) shapiro.test(d$value_trans)), 
         tSW = map(SW, broom::tidy)) %>% 
  select(-data, -SW) %>% 
  unnest(cols = c(tSW))

No transformation helped for the female controls

Test

Condition:

tidy(kruskal.test(value ~ Condition, d))

Sex:

tidy(kruskal.test(value ~ Sex, d))

Condition per sex

d %>% group_by(Sex) %>% 
  nest() %>% 
  mutate(
    t = map(data, ~tidy(kruskal.test(value ~ Condition, .x)))
  ) %>% 
  select(-data) %>% 
  unnest(cols = c(t))

Adrenals

d <- dat %>% 
  group_by(Sex, Condition, SexCondition) %>% 
  filter(prop == 'ORG_Adrenals_adj')
dd <- data.frame(Sex = c('Female', 'Male'), Sig = c('', '**'), y = c(0, 0))
f2_3 <- ggplot(d, aes(x = Condition, y = value, fill = SexCondition)) +
  geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, aes(color = SexCondition)) +
  geom_point(pch = 21, size = 2, 
             position = position_jitterdodge(jitter.width = 0, 
                                             jitter.height = 0, 
                                             dodge.width = 0.5)) + 
  scale_fill_manual(values = obj$four_color) + 
  scale_color_manual(values = obj$four_color) + 
  labs(x = '', y = 'Adrenal weight') +
  facet_grid(.~Sex) + theme(legend.position = 'none') +
  geom_text(data = dd, inherit.aes = F, aes(x = 1.5, y = y, label = Sig), fontface = 'bold', color = 'grey30')
f2_3

Stats:

Summary:

d %>% group_by(Sex, Condition) %>% 
  dplyr::summarize(
    n = n(), 
    n_groups = length(unique(GroupNumber)),
    M = mean(value, na.rm = T), 
    sd = sd(value, na.rm = T)
  )

Shapiro-Wilk test

d %>% group_by(Sex, Condition) %>% nest() %>% 
  mutate(SW = map(data, function(d) shapiro.test(d$value)), 
         tSW = map(SW, broom::tidy)) %>% select(-data, -SW) %>% unnest(cols = c(tSW))

Normality is not violated in any group

Levene’s Test

knitr::kable(car::leveneTest(value ~ Sex * Condition, d))

	Df	F value	Pr(>F)
group	3	0.6155812	0.6069253
	79	NA	NA

Groupwise variances are not heteroscedastic

Test

tidy(aov(value ~ Condition * Sex, data = d)) %>% as.data.frame()

Without interaction:

tidy(aov(value ~ Condition + Sex, data = d))  %>% as.data.frame()

Pairwise

d %>% group_by(Sex) %>% 
  nest() %>% 
  mutate(
    t = map(data, function(d) {
      x = tidy(t.test(d$value ~ d$Condition))
    })
  ) %>% 
  select(-data) %>% 
  unnest(cols = c(t))

Pretty labels

dat <- obj$tests_combined_adj
props <- colnames(dat)[regexpr('_', colnames(dat)) > 0]
props <- props[order(props)]

props_pretty = c(
  'EPM | Center time', 
  'EPM | Closed arm distance', 
  'EPM | Closed arm entries', 
  'EPM | Closed arm mean visit',  
  'EPM | Closed arm time', 
  'EPM | Distance', 
  'EPM | Open arm distance from center', 
  'EPM | Open arm distance', 
  'EPM | Open arm entries', 
  'EPM | Open arm mean visit', 
  'EPM | Open arm time', 
  'EPM | Closed arm preference', 
  
  'NBT | Intact material', 
  'NBT | Nest score', 
  
  'OFT | Mean distance from center', 
  'OFT | Distance',
  'OFT | Inner zone mean speed', 
  'OFT | Inner zone distance', 
  'OFT | Inner zone visits',
  'OFT | Inner zone mean visit duration', 
  
  'OFT | Max. speed', 
  'OFT | Meander',
  'OFT | Outer zone mean speed', 
  'OFT | Outer zone distance', 
  'OFT | Outer zone mean visit duration',
  
  'OFT | Outer zone preference', 
  'OFT | Immobility time',
  
  'PHYS | Adrenal weight',
  'PHYS | Bodyweight change', 
  'PHYS | Corticosterone', 
  'PHYS | Coat state', 
  
  'SPT | Median preference', 
  'SPT | Total preference', 
  
  'SPL | Grooming time',
  'SPL | Latency to groom', 
  
  'TST | Immobility episodes',
  'TST | Immobility time'
)

Figure panel

cowplot::plot_grid(plotlist = list(f2_1, f2_2, f2_3), ncol = 3, align = 'hv',
                   labels = c('b', 'c', 'd', 'e'), label_fontface = 'plain')

if(exportFigures) {
  ggsave('../manuscript/figures/2_Panel.pdf', width = 8.8, height = 2.5, dpi = 1200)
  ggsave('../manuscript/figures/2_Panel.png', width = 8.8, height = 2.5, dpi = 1200)
}

Figure 3 - PCA

PCA

Adjusted for batch effects

dat <- obj$tests_combined_adj
pcaDat <- dat[,colnames(dat) %in% props]
pheno <- dat[,!colnames(dat) %in% props]

pca <- prcomp(pcaDat, center = T, scale. = T)
PCs <- colnames(as.data.frame(pca$x[,1:3]))

p <- cbind(as.data.frame(pheno), as.data.frame(pca$x[,1:3])) %>%
  mutate(
    Sex = factor(Sex, levels = c('Female', 'Male')), 
    SexCondition = paste0(Sex, '_', Condition), 
    SexCondition = factor(SexCondition, 
                          levels = c('Female_Control', 'Male_Control', 'Female_CMS', 'Male_CMS'))
  )
obj$analyses$PCA_OutcomeVars_BatchAdj <- p

Variance explained

eigs <- pca$sdev^2
e <- data.frame(
  PC = colnames(pca$x),
  Proportion = eigs/sum(eigs),
  Cumulative = cumsum(eigs)/sum(eigs))
e <- e[1:5,] %>% 
  mutate(PC = factor(PC, levels = colnames(pca$x)))
f3_1 <- ggplot(e, aes(x = PC, y = Proportion*100)) + 
  geom_line(aes(group = 1), color = 'grey') +
  geom_point(size = 2, pch = 21, fill = obj$colors[2]) +
  scale_y_continuous(limits = c(5, 25), breaks = seq(0, 25, 5)) +
  labs(x = '', y = '% Variance Explained')
f3_1

as_tibble(e)

PCs vs sex & condition

Check Normality

idVars <- c(colnames(pheno), 'SexCondition')
PCs <- c('PC1', 'PC2', 'PC3')

pp <- obj$analyses$PCA_OutcomeVars_BatchAdj %>%
  select(idVars, PCs) %>%
  pivot_longer(all_of(PCs), names_to = "PC") %>%
  group_by(PC, Sex, Condition) %>% nest() %>%
  mutate(
    SW = map(data, ~broom::tidy(shapiro.test(.x$value))), 
  ) %>% 
  unnest(cols = c(SW))

pp

Normality is violated in Male Controls for PC1

Check heteroskedasticity

obj$analyses$PCA_OutcomeVars_BatchAdj %>%
  select(idVars, PCs) %>%
  pivot_longer(all_of(PCs), names_to = "PC") %>%
  group_by(PC) %>% nest() %>%
  mutate(
    t = map(data, ~broom::tidy(car::leveneTest(value ~ Sex * Condition, .x)))
  ) %>%
  unnest(cols = c(t)) %>% 
  filter(term == 'group')

All variances are homogenous

idVars <- c(colnames(pheno), 'SexCondition')
PCs <- c('PC1', 'PC2', 'PC3')

pp <- obj$analyses$PCA_OutcomeVars_BatchAdj %>%
  dplyr::select(idVars, PCs) %>%
  gather(PC, value, -idVars) %>%
  group_by(PC) %>% nest() %>%
  mutate(
    m = map2(data, PC, function(d, p) {
      summary(aov(value ~ Condition * Sex, d))
    })
  ) %>% dplyr::select(-data) 

plots <- obj$analyses$PCA_OutcomeVars_BatchAdj %>%
  dplyr::select(idVars, PCs) %>%
  gather(PC, value, -idVars) %>%
  group_by(PC) %>% nest() %>% mutate(
    plot = map2(data, PC, function(d, p) {
      ggplot(d, aes(x = Condition, y = value, fill = SexCondition)) +
        geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, aes(color = SexCondition),
                     position = position_dodge(0.6)) +
        geom_point(pch = 21, position = position_dodge(0.6), size = 2) +
        labs(x = '', y = p) + facet_grid(~Sex) +
        scale_fill_manual(values = obj$four_color) +
        scale_color_manual(values = obj$four_color) +
        theme(legend.position = 'none')
    })
  )

PC1:

pp$m[[1]]

##               Df Sum Sq Mean Sq F value Pr(>F)
## Condition      1    5.0   4.959   0.608  0.438
## Sex            1    0.1   0.113   0.014  0.907
## Condition:Sex  1    7.4   7.437   0.912  0.342
## Residuals     82  669.0   8.158

PC2:

pp$m[[2]]

##               Df Sum Sq Mean Sq F value        Pr(>F)    
## Condition      1   0.25    0.25   0.077         0.782    
## Sex            1 156.62  156.62  48.269 0.00000000081 ***
## Condition:Sex  1   4.02    4.02   1.238         0.269    
## Residuals     82 266.06    3.24                          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

PC3:

pp$m[[3]]

##               Df Sum Sq Mean Sq F value Pr(>F)  
## Condition      1   23.2  23.231   5.920 0.0171 *
## Sex            1    7.4   7.371   1.879 0.1742  
## Condition:Sex  1    0.4   0.423   0.108 0.7435  
## Residuals     82  321.8   3.924                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cowplot::plot_grid(plotlist = plots$plot, ncol = 2)

f3_2 <- plots$plot[plots$PC == 'PC1'][[1]] + ylab('PC1 score')

PC1

baseline <- h_cum
pp <- left_join(baseline, obj$analyses$PCA_OutcomeVars_BatchAdj)
pp$Condition <- factor(pp$Condition, levels = c('Control', 'CMS'))
f3_3 <- ggplot(pp, aes(x = DS, y = PC1)) +
  geom_smooth(method = 'lm', se = F, lty = 2, aes(col = SexCondition)) +
  geom_point(pch = 21, size = 2, aes(fill = SexCondition)) +
  labs(y = "PC1 Score", x = "Baseline David's Score") +
  facet_grid(Condition~Sex) +
  scale_fill_manual(values = obj$four_color) +
  scale_color_manual(values = obj$four_color) +
  theme(legend.position = 'none')
f3_3

Stats

summary(lm(PC1 ~ DS * Sex, pp[pp$Condition == 'CMS',]))

## 
## Call:
## lm(formula = PC1 ~ DS * Sex, data = pp[pp$Condition == "CMS", 
##     ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.2597 -1.7029 -0.3459  1.8846  6.4257 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   -3.294      2.273  -1.449   0.1545  
## DS             2.145      1.464   1.464   0.1504  
## SexMale        7.165      2.796   2.562   0.0140 *
## DS:SexMale    -4.348      1.773  -2.453   0.0183 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.856 on 43 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.1491, Adjusted R-squared:  0.08975 
## F-statistic: 2.512 on 3 and 43 DF,  p-value: 0.07119

summary(aov(PC1 ~ DS * Sex, pp[pp$Condition == 'CMS',]))

##             Df Sum Sq Mean Sq F value Pr(>F)  
## DS           1    7.9    7.90   0.969 0.3305  
## Sex          1    4.5    4.49   0.551 0.4619  
## DS:Sex       1   49.1   49.06   6.016 0.0183 *
## Residuals   43  350.7    8.15                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 1 observation deleted due to missingness

Correlations with behaviors

p <- obj$analyses$PCA_OutcomeVars_BatchAdj
PC <- 'PC1'
p$PC <- p[[PC]]

cors <- tibble(
  feature = colnames(pcaDat)
) %>% 
  mutate(
    data = map(feature, ~ pcaDat[[.x]]),
    cor_test = map(data, ~ cor.test(.x, p$PC, method = 'spearman')), 
    rho = map_dbl(cor_test, ~.x$estimate), 
    pval = map_dbl(cor_test, ~.x$p.value)
  )
cors$feature <- factor(cors$feature, levels = cors$feature[order(abs(cors$rho), decreasing = F)])
cors <-  cors[order(abs(cors$rho), decreasing = T), ]
cors$Sign <-  ifelse(cors$rho < 0, yes = 'negative', no = 'positive')
cors$Sign <-  factor(cors$Sign, levels = c('positive', 'negative'))
cors$Sig <- p.adjust(cors$pval, method = 'bonferroni') <= 0.05
cors$feature_pretty <-  plyr::mapvalues(cors$feature, props, props_pretty)

f3_4 <-  ggplot(cors, aes(x = abs(rho), 
                          y = feature_pretty, 
                          fill = Sign)) +
  geom_bar(stat = 'identity', width = 0.15) +  
  geom_point(size = 2.5, pch = 19, alpha = 1, aes(col = Sign)) + 
  geom_point(data = subset(cors, Sig), size = 2.5, pch = 1, color = 'black', stroke = 1.2) +
  scale_fill_manual(values = c(hm.palette(9)[8], hm.palette(9)[1]), name = '', guide = F) +
  scale_color_manual(values = c(hm.palette(9)[8], hm.palette(9)[1]), name = '') +
  labs(y = NULL, x = 'Absolute Correlation') + 
  scale_x_continuous(expand = c(0, 0, 0, 0.1)) +
  scale_y_discrete(expand = c(0, 0.5, 0, 2)) +
  theme(plot.title = element_text(hjust = 0.5), 
        panel.border = element_blank(),
        legend.position = c(.65, .05))  
f3_4

Loadings

p <- pca$rotation[,1] %>%
  as.data.frame()
colnames(p) <- 'rho'
p$feature <- rownames(p)
p$feature <- factor(p$feature, levels = p$feature[order(abs(p$rho), decreasing = F)])
p <-  p[order(abs(p$rho), decreasing = T), ]

p$Sign <-  ifelse(p$rho < 0, yes = 'neg', no = 'pos')
p$Sign <-  factor(p$Sign, levels = c('pos', 'neg'))
p$feature_pretty <-  plyr::mapvalues(p$feature, props, props_pretty)

ggplot(p, aes(y = abs(rho), x = feature_pretty, col = Sign, fill = Sign)) +
  geom_bar(stat = 'identity', width = 0.01) +  
  geom_point(size = 2, pch = 21, alpha = 1) + 
  scale_fill_manual(values = c(hm.palette(9)[8], hm.palette(9)[1]), name = '') +
  scale_color_manual(values = c(hm.palette(9)[8], hm.palette(9)[1]), name = '') +
  xlab('') + ylab("Loading") + coord_flip() + 
  theme(plot.title = element_text(hjust = 0.5), legend.position = 'top')

Selected features

dat <- obj$tests_combined_adj
idVars <- c('Color', 'Sex', 'Condition', 'Batch', 'RealMouseNumber', 
            'GroupNumber', 'MouseID', 'GroupID', 'uGroup', 'GroupType')
props <- colnames(dat)[!colnames(dat) %in% idVars]
props <- props[order(props)]
baseline <- h_cum %>% rename(Dominance = DS)
d <- left_join(dat, baseline)

idVars <- c(idVars, 'Dominance')
plots <- d  %>% 
  select(idVars, props) %>% 
  gather(prop, value, -idVars) %>% 
  group_by(prop) 
plots$prop <- plyr::mapvalues(plots$prop, props, props_pretty)
plots$Sex <- factor(plots$Sex, levels = c('Female', 'Male'))
plots$Condition <- factor(plots$Condition, levels = c('Control', 'CMS'))
plots$SexCondition <- paste0(plots$Sex, '_', dat$Condition)
plots$SexCondition <- factor(plots$SexCondition, 
                             levels = c('Female_Control', 'Male_Control', 'Female_CMS', 'Male_CMS'))

selected_props <- c('OFT | Distance', 'EPM | Open arm distance')

res <- lapply(selected_props, function(i) {
  pp <- subset(plots, prop == i)
  l <- c(min(pp$value), max(pp$value))
  p1 <- ggplot(pp, aes(x = Condition, y = value, fill = SexCondition)) +
    geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, aes(color = SexCondition)) +
    geom_point(pch = 21, size = 2, 
               position = position_jitterdodge(jitter.width = 0, 
                                               jitter.height = 0, dodge.width = 0.5)) + 
    scale_fill_manual(values = obj$four_color) + 
    scale_color_manual(values = obj$four_color) + 
    facet_grid(Sex~.) + labs(x = '', y = i) +
    scale_y_continuous(limits = l) +
    theme(legend.position = 'none') + 
    theme(strip.text.y = element_blank())
  
  p2 <- ggplot(pp[pp$Condition == 'CMS',], aes(x = Dominance, y = value, fill = Sex, color = Sex)) +
    geom_smooth(method = 'lm', se = F, lty = 2) +
    geom_point(pch = 21, size = 2, color = 'grey20') +
    scale_fill_manual(values = obj$four_color[3:4]) +
    scale_color_manual(values = obj$four_color[3:4]) +
    scale_y_continuous(limits = l) +
    labs(y = '', x = "Baseline David's Score") +
    facet_grid(Sex~Condition) +
    theme(legend.position = 'none', 
          axis.text.y = element_blank(),
          plot.margin = margin(0,0,0,-10, "pt")
    )
  
  cowplot::plot_grid(plotlist = list(p1, p2), axis = 'tb',
                     align = 'h', rel_widths = c(1,1.5))
})
f3_5 <- cowplot::plot_grid(plotlist = res, ncol = 2, align = 'v', rel_widths = c(1,1),
                           labels = c('e', 'f'), label_fontface = 'plain')
f3_5

Figure panel

row1 <- cowplot::plot_grid(plotlist = list(f3_1, f3_2), ncol = 2, axis = 't', align = 'h',
                           labels = c('a', 'b'), label_fontface = 'plain', rel_widths = c(0.7,1))

row2 <- cowplot::plot_grid(plotlist = list(f3_3), ncol = 1, align = 'hv',
                           labels = c('c'), label_fontface = 'plain')

col1 <- cowplot::plot_grid(plotlist = list(row1, row2), ncol = 1, axis = 't', align = 'h',
                           rel_widths = c(0.5,1), rel_heights = c(1, 2.1))

col2 <- cowplot::plot_grid(plotlist = list(f3_4), ncol = 1, align = 'hv',
                           labels = c('d'), label_fontface = 'plain')

l1 <- cowplot::plot_grid(plotlist = list(col1, col2), align = 'hv', axis = 't',
                         ncol = 2, rel_widths = c(1, 0.9),
                         label_fontface = 'plain')

cowplot::plot_grid(plotlist = list(l1, f3_5), align = 'hv', ncol = 1, rel_heights = c(1, 0.7))

if(exportFigures) {
  ggsave('../manuscript/figures/3_Panel.pdf', width = 8.5, height = 11, dpi = 1200)
  ggsave('../manuscript/figures/3_Panel.png', width = 8.5, height = 11, dpi = 1200)
}

Figure 4

Supplement

DS vs SB behaviors

Cumulative David’s score based on chases
Mean of SB variables over days 2,3,4

p <- subset(obj$profiles, !Suspect & Day %in% 2:4) %>% 
  select(RealMouseNumber, Batch, all_props) %>% 
  pivot_longer(cols = all_of(all_props), names_to = "prop") %>% 
  filter(!is.na(value)) %>%
  filter(prop != 'NAChaseRate') %>%
  group_by(prop, Batch) %>%
  mutate(
    value = RNOmni::rankNorm(value, k = 3/8)
  ) %>% 
  group_by(RealMouseNumber, prop) %>% 
  summarize(value = mean(value)) %>% 
  pivot_wider(names_from = 'prop')

d <- left_join(h_cum, p)
idVars <- c('RealMouseNumber', 'GroupNumber', 'Sex', 'DS')
props <- all_props[!all_props == 'NAChaseRate']

p <- d %>% 
  select(idVars, props) %>% 
  pivot_longer(cols = !idVars, names_to = "prop", values_to = "value") %>% 
  group_by(prop, Sex) %>%
  nest() %>% 
  mutate(
    M = map(data, ~tidy(cor.test(.x$value, .x$DS, method = 'spearman'))) 
  ) %>% 
  select(-data) %>% 
  unnest(cols = M) %>%
  ungroup() %>% # added becuase prop was a grouping variable and thus not modifiable
  mutate(
    log10pval = -log10(p.value),
    prop = fct_reorder(prop, abs(estimate)), 
    sign = ifelse(estimate > 0, '+', '-'), 
    sig_nominal = p.value <= 0.05
  ) %>% 
  group_by(Sex) %>% 
  mutate(qval = p.adjust(p.value, method = 'BH')) %>% 
  arrange(qval)

plots <- d  %>% 
  select(idVars, props) %>% 
  gather(prop, value, -idVars) %>% 
  group_by(prop) %>% 
  nest() %>% 
  mutate(
    plot = pmap(list(data, prop), 
                function(d, p) {
                  ggplot(d, aes(x = DS, y = value, fill = Sex)) +
                    geom_smooth(method = 'lm', se = F, lty = 2, aes(col = Sex)) +
                    geom_point(pch = 21, size = 2) +
                    labs(y = p, x = "Baseline David's Score") +
                    facet_grid(.~Sex) +
                    scale_fill_manual(values = obj$colors_sex) +
                    scale_color_manual(values = obj$colors_sex) +
                    theme(legend.position = 'none')
                })) %>% dplyr::select(-data)

p <- left_join(p, plots)

displaying only behaviors that showed a q-val significant correlation with DS

p <- subset(p, !is.na(estimate))
p %>% group_by(prop) %>% 
  summarize(is_sig = sum(qval <= 0.05) > 0) %>% 
  group_by(is_sig) %>% tally()

p[p$prop == 'DistanceOutside',]

p[p$prop == 'FractionOfTimeOutside',]

p[p$prop == 'GridEntropy6x6',]

p$rho <- abs(p$estimate)
p$rho[p$qval > 0.05] <- NA

p <- p %>% 
  filter(qval < 0.05)

p$sign <- factor(p$sign, levels = c('+', '-'))
p$prop_pretty <- plyr::mapvalues(p$prop, all_props, all_props_pretty)

ord <- p$prop_pretty[p$Sex == 'Male'][order(p$estimate[p$Sex == 'Male'], decreasing = F)]
p$prop_pretty <- factor(p$prop_pretty, levels = ord)
p <- subset(p, !is.na(prop_pretty))
p$Sex <- factor(p$Sex, levels = c('Male', 'Female'))
ggplot(p[!is.na(p$sign),], aes(y = Sex, x =  prop_pretty, col = sign, fill = sign)) +
  geom_point(pch = 22, aes(size = rho)) +
  scale_fill_manual(values = c(hm.palette(9)[8], hm.palette(9)[1]), name = 'Direction') +
  scale_color_manual(values = c(hm.palette(9)[8], hm.palette(9)[1]), name = 'Direction') +
  scale_size_continuous(name = 'Correlation') +
  labs(x = '', y = '') + coord_flip() + 
  theme(plot.title = element_text(hjust = 0.5), legend.position = 'right', 
        panel.border = element_blank())

if(exportFigures) {
  ggsave('../manuscript/figures/Supp_DavidsScore_vs_SocialBox.pdf', height = 8, width = 7.2, dpi = 1200)
  ggsave('../manuscript/figures/Supp_DavidsScore_vs_SocialBox.png', height = 8, width = 7.2, dpi = 1200)
}

Hierarchy Panel

Hierarchy properties based on chases
Daily & mean of days 1:4 and day 5 separately

domProps <- c('Steepness', 'Despotism', 'DirectionalConsistency', 'Landau')
domPropsPretty <- c('Steepness', 'Despotism', 'Directional Consistency', "Landau's modified h'")

h <- explore_hierarchies(subset(obj$master, Day %in% 1:5), 
                         idVars = c('Sex', 'Condition'), 
                         DS_only = F,
                         verbose = F, daily = T)
h$uGroup <- paste0(h$GroupNumber, h$Sex)
plots <- lapply(1:(length(domProps)-1), function(i) {
  prop <- domProps[i]
  h$prop <- h[[prop]]
  lims <- c(min(h$prop, na.rm = T), max(h$prop, na.rm = T))
  p1 <- ggplot(h, aes(x = as.factor(Day), y = prop, fill = Sex)) +
    geom_line(aes(color = Sex, group = uGroup), lty = 1, size = 0.2, position = position_dodge(0)) +
    geom_point(aes(group = uGroup, col = Sex), pch = 19, size = 1.2, position = position_dodge(0)) +
    labs(x = '', y = domPropsPretty[i])  + 
    scale_color_manual(values = obj$colors_sex) +
    scale_fill_manual(values = obj$colors_sex) +
    scale_y_continuous(limits = lims) +
    theme(legend.position = 'none', axis.title.x=element_blank())
  
  p2 <- h %>% dplyr::filter(Day %in% 1:4) %>% 
    group_by(GroupNumber, Condition, Sex) %>% 
    dplyr::summarise(prop = mean(prop, na.rm = T)) %>% 
    ggplot(aes(x = Sex, y = prop, fill = Sex)) +
    geom_boxplot(width = 0.5, alpha = 0.8, outlier.shape = NA, aes(col = Sex)) +
    geom_point(pch = 21, size = 2.5, position = position_dodge(0.5)) +
    labs(x = '', y = '')  + 
    scale_y_continuous(limits = lims) +
    scale_fill_manual(values = obj$colors_sex) +
    scale_color_manual(values = obj$colors_sex) +
    theme(legend.position = 'none', 
          axis.text.y = element_blank(),
          axis.title=element_blank()
    )  +
    ggtitle('Baseline')
  
  p3 <- h %>% dplyr::filter(Day %in% 5) %>% 
    group_by(GroupNumber, Condition, Sex) %>% 
    dplyr::summarise(prop = mean(prop, na.rm = T)) %>% 
    ggplot(aes(x = Sex, y = prop, fill = Sex)) +
    geom_boxplot(width = 0.5, alpha = 0.8, outlier.shape = NA, aes(col = Sex)) +
    geom_point(pch = 21, size = 2.5, position = position_dodge(0.5)) +
    labs(x = '', y = '')  + 
    scale_y_continuous(limits = lims) +
    scale_fill_manual(values = obj$colors_sex) +
    scale_color_manual(values = obj$colors_sex) +
    theme(legend.position = 'none', 
          axis.text.y = element_blank(),
          axis.title=element_blank()
    )  +
    ggtitle('Acute stress')
  cowplot::plot_grid(plotlist = list(p1, p2, p3), ncol = 3, align = 'h',
                     label_fontface = 'plain',
                     rel_widths = c(3.4, 1, 1))
})

lastPlot <- lapply(length(domProps), function(i) {
  prop <- domProps[i]
  h$prop <- h[[prop]]
  lims <- c(min(h$prop, na.rm = T), max(h$prop, na.rm = T))
  p1 <- ggplot(h, aes(x = as.factor(Day), y = prop, fill = Sex)) +
    geom_boxplot(width = 0.7, alpha = 0.8, outlier.shape = NA, aes(color = Sex)) +
    geom_dotplot(binaxis='y', stackdir='center', dotsize = 1, position = position_dodge(0.7)) +
    labs(x = 'Day', y = domPropsPretty[i])  + 
    scale_fill_manual(values = obj$colors_sex) +
    scale_color_manual(values = obj$colors_sex) +
    scale_y_continuous(limits = lims) +
    theme(legend.position = 'none')
  
  p2 <- h %>% dplyr::filter(Day %in% 1:4) %>% 
    group_by(GroupNumber, Condition, Sex) %>% 
    dplyr::summarise(prop = mean(prop, na.rm = T)) %>% 
    ggplot(aes(x = Sex, y = prop, fill = Sex)) +
    geom_boxplot(width = 0.5, alpha = 0.8, outlier.shape = NA, aes(col = Sex)) +
    # geom_point(pch = 21, size = 2.5, position = position_dodge(0.5)) +
    geom_dotplot(binaxis='y', stackdir='center', dotsize = 1.3) +
    labs(x = '', y = '')  + 
    scale_y_continuous(limits = lims) +
    scale_fill_manual(values = obj$colors_sex) +
    scale_color_manual(values = obj$colors_sex) +
    theme(legend.position = 'none', axis.text.y = element_blank(), 
          axis.title.y=element_blank()) +
    ggtitle('Baseline')
  
  p3 <- h %>% dplyr::filter(Day %in% 5) %>% 
    group_by(GroupNumber, Condition, Sex) %>% 
    dplyr::summarise(prop = mean(prop, na.rm = T)) %>% 
    ggplot(aes(x = Sex, y = prop, fill = Sex)) +
    geom_boxplot(width = 0.5, alpha = 0.8, outlier.shape = NA, aes(col = Sex)) +
    geom_dotplot(binaxis='y', stackdir='center', dotsize = 1.3) +
    labs(x = '', y = '')  + 
    scale_y_continuous(limits = lims) +
    scale_fill_manual(values = obj$colors_sex) +
    scale_color_manual(values = obj$colors_sex) +
    theme(legend.position = 'none', axis.text.y = element_blank(), 
          axis.title.y=element_blank()) + 
    ggtitle('Acute stress')
  
  cowplot::plot_grid(plotlist = list(p1, p2, p3), ncol = 3, align = 'h', 
                     rel_widths = c(3.4, 1, 1))
})
plots[[length(domProps)]] <- lastPlot[[1]]
cowplot::plot_grid(plotlist = plots, ncol = 1, align = 'v', axis = 'tblr', 
                   labels = c('a', 'b', 'c', 'd'), label_fontface = 'plain',
                   rel_heights = c(rep(1, length(domProps)-1), 1.1))

if(exportFigures) {
  ggsave('../manuscript/figures/Supp_HierarchyPanel.pdf', width = 8.5, height = 11, dpi = 1200)
  ggsave('../manuscript/figures/Supp_HierarchyPanel.png', width = 8.5, height = 11, dpi = 1200)
}

Stats:

domProps <- c('Steepness', 'Despotism', 'DirectionalConsistency', 'Landau')
domPropsPretty <- c('Steepness', 'Despotism', 'Directional Consistency', "Landau's modified h'")
h$Stage <- plyr::mapvalues(h$Day, 1:5, c(rep('Baseline', 4), 'Acute stress'))
h$Stage <- factor(h$Stage, levels = c('Baseline', 'Acute stress'))

h <- h %>% group_by(uGroup, GroupNumber, Sex, Stage) %>% 
  dplyr::summarise_at(.vars = domProps, mean, na.rm = T)

h %>% gather(prop, value, -uGroup, -GroupNumber, -Sex, -Stage) %>% 
  group_by(prop, Sex, Stage) %>% nest() %>% 
  dplyr::mutate(
    SW = map(data, function(d) shapiro.test(d$value)), 
    tSW = map(SW, broom::tidy)
  ) %>% 
  select(-data, -SW) %>% 
  unnest(cols = tSW)

Normality violated only for Landau’s h in female groups during acute stress

Check heteroskedasticity

h %>% gather(prop, value, -uGroup, -GroupNumber, -Sex, -Stage) %>% 
  group_by(prop) %>% nest() %>% 
  dplyr::mutate(
    t = map(data, function(d) {
      car::leveneTest(value ~ Sex * Stage, d)
    })
  )%>% 
  select(-data) %>% 
  unnest(cols = t)

Variances are heteroskedastic for despotism

Test

Steepness

summary(aov(Steepness ~ Sex * Stage + Error(uGroup), data = h))

## 
## Error: uGroup
##           Df Sum Sq Mean Sq F value  Pr(>F)   
## Sex        1 0.2831 0.28311    12.6 0.00201 **
## Residuals 20 0.4494 0.02247                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Within
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## Stage      1 0.00377 0.003772   0.574  0.457
## Sex:Stage  1 0.00000 0.000002   0.000  0.986
## Residuals 20 0.13141 0.006571

Despotism

summary(aov(Despotism ~ Sex * Stage + Error(uGroup), data = h))

## 
## Error: uGroup
##           Df Sum Sq Mean Sq F value  Pr(>F)   
## Sex        1 0.4423  0.4423   10.84 0.00364 **
## Residuals 20 0.8159  0.0408                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Within
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## Stage      1 0.00464 0.004641   0.754  0.396
## Sex:Stage  1 0.01281 0.012810   2.080  0.165
## Residuals 20 0.12317 0.006158

DirectionalConsistency

summary(aov(DirectionalConsistency ~ Sex * Stage + Error(uGroup), data = h))

## 
## Error: uGroup
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## Sex        1 0.6687  0.6687   21.16 0.000173 ***
## Residuals 20 0.6321  0.0316                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Within
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## Stage      1 0.00725 0.007251   1.070  0.313
## Sex:Stage  1 0.00952 0.009523   1.405  0.250
## Residuals 20 0.13557 0.006779

Interesting effect of stage

Males:

summary(aov(DirectionalConsistency ~ Stage + Error(uGroup), data = h[h$Sex == 'Male', ]))

## 
## Error: uGroup
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  9 0.5112  0.0568               
## 
## Error: Within
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## Stage      1 0.01677 0.016766   2.385  0.157
## Residuals  9 0.06328 0.007031

Females:

summary(aov(DirectionalConsistency ~ Stage + Error(uGroup), data = h[h$Sex == 'Female', ]))

## 
## Error: uGroup
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 11  0.121   0.011               
## 
## Error: Within
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## Stage      1 0.00001 0.000008   0.001  0.972
## Residuals 11 0.07229 0.006572

ggplot(h, aes(x = Stage, y = DirectionalConsistency, fill = Sex)) +
  geom_line(aes(group = uGroup, col = Sex)) +
  geom_point(size = 2, pch = 21) +
  labs(x = '')  + 
  scale_fill_manual(values = obj$colors_sex) +
  scale_color_manual(values = obj$colors_sex) +
  facet_grid(~Sex) + 
  theme(legend.position = 'none', axis.text.y = element_blank(), 
        axis.title.y=element_blank())

Landau

summary(aov(Landau ~ Sex * Stage + Error(uGroup), data = h))

## 
## Error: uGroup
##           Df Sum Sq Mean Sq F value Pr(>F)  
## Sex        1 0.2082 0.20824   3.822 0.0647 .
## Residuals 20 1.0896 0.05448                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Within
##           Df Sum Sq  Mean Sq F value Pr(>F)
## Stage      1 0.0015 0.001546   0.072  0.791
## Sex:Stage  1 0.0245 0.024526   1.145  0.297
## Residuals 20 0.4283 0.021416

Baseline non-parametric:

kruskal.test(Landau ~ Sex, data = h[h$Stage == 'Baseline',])

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Landau by Sex
## Kruskal-Wallis chi-squared = 8.2458, df = 1, p-value = 0.004085

Acute stress non-parametric:

kruskal.test(Landau ~ Sex, data = h[h$Stage == 'Acute stress',])

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Landau by Sex
## Kruskal-Wallis chi-squared = 0.19081, df = 1, p-value = 0.6622

PCs 2-3

idVars <- c(colnames(pheno), 'SexCondition')
PCs <- c('PC2', 'PC3')

pp <- obj$analyses$PCA_OutcomeVars_BatchAdj %>%
  select(idVars, PCs) %>%
  pivot_longer(names_to = "PC", values_to = "value", cols = !idVars) %>% 
  group_by(PC) %>% 
  nest() %>%
  mutate(
    m = map2(data, PC, function(d, p) {
      summary(aov(value ~ Condition * Sex, d))
    })
  ) %>% dplyr::select(-data) 

plots <- obj$analyses$PCA_OutcomeVars_BatchAdj %>%
  select(idVars, PCs) %>%
  pivot_longer(names_to = "PC", values_to = "value", cols = !idVars) %>%
  group_by(PC) %>% 
  nest() %>% 
  mutate(
    plot = map2(data, PC, function(d, p) {
      ggplot(d, aes(x = Condition, y = value, fill = SexCondition)) +
        geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, aes(color = SexCondition),
                     position = position_dodge(0.6)) +
        geom_point(pch = 21, position = position_dodge(0.6), size = 2.5) +
        labs(x = '', y = p) + facet_grid(~Sex) +
        scale_fill_manual(values = obj$four_color) +
        scale_color_manual(values = obj$four_color) +
        theme(legend.position = 'none')
    })
  )

Effect of Condition on PC2

pp$m[[1]]

##               Df Sum Sq Mean Sq F value        Pr(>F)    
## Condition      1   0.25    0.25   0.077         0.782    
## Sex            1 156.62  156.62  48.269 0.00000000081 ***
## Condition:Sex  1   4.02    4.02   1.238         0.269    
## Residuals     82 266.06    3.24                          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Effect of Sex on PC3

pp$m[[2]]

##               Df Sum Sq Mean Sq F value Pr(>F)  
## Condition      1   23.2  23.231   5.920 0.0171 *
## Sex            1    7.4   7.371   1.879 0.1742  
## Condition:Sex  1    0.4   0.423   0.108 0.7435  
## Residuals     82  321.8   3.924                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cowplot::plot_grid(plotlist = plots$plot, ncol = 2)

if(exportFigures) {
  ggsave('../manuscript/figures/Supp_PCs.pdf', width = 8.5, height = 4, dpi = 1200)
  ggsave('../manuscript/figures/Supp_PCs.png', width = 8.5, height = 4, dpi = 1200)
}

DS correlations with behaviors

dat <- obj$tests_combined_adj_long
dat$value[dat$has_NA] <- NA
idVars <- c('Color', 'Sex', 'Condition', 'Batch', 'RealMouseNumber', 
            'GroupNumber', 'MouseID', 'GroupID', 'uGroup', 'GroupType')
props <- unique(dat$prop)
props <- props[order(props)]
baseline <- h_cum %>% rename(Dominance = DS)

d <- left_join(dat, baseline)
idVars <- c(idVars, 'Dominance', 'Rank')
d$prop <- plyr::mapvalues(d$prop, props, props_pretty)
d$Sex <- factor(d$Sex, levels = c('Female', 'Male'))
d$Condition <- factor(d$Condition, levels = c('Control', 'CMS'))
d$SexCondition <- paste0(d$Sex, '_', dat$Condition)
d$SexCondition <- factor(d$SexCondition, 
                         levels = c('Female_Control', 'Male_Control', 'Female_CMS', 'Male_CMS'))
dd <- d %>%
  mutate(
    Rank = plyr::mapvalues(Rank, c('alpha', 'beta', 'gamma', 'delta'), 1:4),
    Rank = as.numeric(as.character(Rank))
    ) %>%
  group_by(prop, Sex, Condition) %>%
  filter(!is.na(value)) %>%
  nest() %>% 
  mutate(
    n = map_dbl(data, ~length(.x$value)),
    mod = map(data, ~cor.test(.x$value, .x$Dominance, method = 'pearson')),
    rho = map_dbl(mod, ~.x$estimate), 
    p.value = map_dbl(mod, ~.x$p.value)
  ) %>% 
  group_by(Sex, Condition) %>% 
  mutate(padj = p.adjust(p.value, "BH"))

dd$stars <- ifelse(dd$p.value > 0.05, " ", ifelse(dd$p.value > 0.01, "*", "**"))
dd$lab <- paste0(round(dd$rho, 2), ' (', round(dd$p.value, 3), ')')
dd$lab_ns <- paste0(round(dd$rho, 2))
ggplot(dd, aes(x = Condition, y = prop, fill = rho)) +
  geom_tile(color = 'white', size = .5) +
  # geom_text(aes(label = stars), color = 'white', size = 4, vjust = 0.8) +
  geom_text(data = subset(dd, p.value < 0.05), aes(label = lab), 
            color = 'white', size = 2.9, vjust = 0.5) +
  geom_text(data = subset(dd, p.value > 0.05), aes(label = lab_ns), 
            color = 'white', size = 2.5, vjust = 0.5) +
  scale_fill_gradient2(low = '#4575b4', high = '#d73027') +
  labs(x = NULL, y = NULL, fill = 'Correlation') +
  facet_grid(~Sex) +
  theme(
    legend.position = 'right', legend.direction = "vertical", 
        panel.border = element_blank())

if(exportFigures) {
  ggsave('../manuscript/figures/Supp_DS_beh.pdf', width = 6.5, height = 10, dpi = 1200)
  ggsave('../manuscript/figures/Supp_DS_beh.png', width = 6.5, height = 10, dpi = 1200)
}

Figures

Stoyo Karamihalev & Elena Brivio

Overview

master

bodyweights

organs

tests

tests_combined

Outlier removal

Figure 1 - Baseline Dominance

Baseline Cumulative DS

Correlations with SB behaviors

Stats

Sum of daily changes

Change/no change analysis

DS baseline v stress

Stats

Chases after stress

Stats

Figure panel

Figure 2 - CMS outcomes

Physiological Outcomes

Bodyweight

Coat state

Adrenals

Pretty labels

Figure panel

Figure 3 - PCA

PCA

Variance explained

PCs vs sex & condition

PC1

Correlations with behaviors

Loadings

Selected features

Figure panel

Figure 4

Supplement

DS vs SB behaviors

Hierarchy Panel

PCs 2-3

DS correlations with behaviors

`master`

`bodyweights`

`organs`

`tests`

`tests_combined`