revision

straussed · Mar 7, 2023 · 6eb9b6c · 6eb9b6c
1 parent 49e8ff7
commit 6eb9b6c
Show file tree

Hide file tree

Showing 12 changed files with 149 additions and 37 deletions.
diff --git a/.DS_Store b/.DS_Store
diff --git a/0_pull_data.R → 00_pull_data.R b/0_pull_data.R → 00_pull_data.R
@@ -9,10 +9,10 @@ library(hyenadata) ### Version 1.2.92
 library(dplyr)
 library(here)
 
-data(tblFemaleRanks)
+load('0_create_ranks/female_ranks.RData')
 data(tblLifeHistory.wide)
 data(tblHyenas)
-
+tblFemaleRanks <- left_join(tblFemaleRanks, tblHyenas[,c('id', 'sex')])
 
 ### Assemble annual reproductive success data
 reproduction <- tblLifeHistory.wide %>%
@@ -29,13 +29,15 @@ reproduction <- tblLifeHistory.wide %>%
   summarize(ars = length(id))
 
 ars <- left_join(tblFemaleRanks, reproduction, by = c('id' = 'mom', 'year'))
+ars <- ars %>% filter(sex == 'f') %>% select(clan, year, id, rank, ars)
 ars[is.na(ars$ars),]$ars <- 0
 
 
 ### Assemble rank data
 tblLifeHistory.wide <- left_join(tblLifeHistory.wide, tblHyenas[,c('id', 'sex')])
 
 ranks <- tblFemaleRanks[c('clan', 'year', 'id', 'rank')]
+
 ranks <- ranks %>%
   rename(period = year) %>%
   group_by(clan, period) %>%

diff --git a/0_create_ranks/1.hyena_intx_data.RData b/0_create_ranks/1.hyena_intx_data.RData
diff --git a/0_create_ranks/2.get_ranks_dynarank.R b/0_create_ranks/2.get_ranks_dynarank.R
@@ -0,0 +1,108 @@
+################################################################################
+#                       Get ranks using DynaRank                               #
+#                                                                              #
+#                                                                              #
+#                           By Eli Strauss                                     #
+#                                                                              #
+#                            March 2023                                        #
+################################################################################
+
+
+rm(list = ls())
+options(stringsAsFactors = FALSE)
+setwd('C:\\Users\\strau\\Documents\\code\\inheritance\\0_create_ranks')
+
+library(DynaRankR)
+library(dplyr)
+
+load('1.hyena_intx_data.RData')
+
+set.seed(1989)
+
+
+################################################################################
+#
+#    Talek
+#
+talek_ranks <- informed_matreorder(contestants = contestants.talek,
+                                   convention = 'mri',
+                                   n = 50,
+                                   shuffles = 30,
+                                   require.corroboration = TRUE,
+                                   initial.ranks = initial.ranks.talek,
+                                   interactions = interactions.talek)
+
+plot_ranks(talek_ranks)
+
+
+
+################################################################################
+#
+#    North
+#
+
+north_ranks <- informed_matreorder(contestants = contestants.serena.n,
+                        convention = 'mri',
+                        n = 100,
+                        shuffles = 30,
+                        require.corroboration = TRUE,
+                        initial.ranks = initial.ranks.serena.n,
+                        interactions = interactions.serena.n)
+
+
+plot_ranks(north_ranks)
+
+################################################################################
+#
+#    South
+#
+
+south_ranks <- informed_matreorder(contestants = contestants.serena.s,
+                        convention = 'mri',
+                        n = 100,
+                        shuffles = 30,
+                        require.corroboration = TRUE,
+                        initial.ranks = initial.ranks.serena.s,
+                        interactions = interactions.serena.s)
+
+plot_ranks(south_ranks)
+
+################################################################################
+#
+#    Happy Zebra
+#
+
+hz_ranks <- informed_matreorder(contestants = contestants.happy.zebra,
+                     convention = 'mri',
+                     n = 100,
+                     shuffles = 30,
+                     require.corroboration = TRUE,
+                     initial.ranks = initial.ranks.happy.zebra,
+                     interactions = interactions.happy.zebra)
+
+plot_ranks(hz_ranks)
+
+
+
+
+## Combine data
+talek_ranks$clan <- 'talek'
+north_ranks$clan <- 'serena.n'
+south_ranks$clan <- 'serena.s'
+hz_ranks$clan <- 'happy.zebra'
+
+#filter to 2017
+ranks <- rbind(filter(talek_ranks, period <= 2017),
+               filter(north_ranks, period <= 2017),
+               filter(south_ranks, period <= 2017),
+               filter(hz_ranks,period <= 2017))
+
+tblFemaleRanks <- ranks %>%
+  rename(year = period) %>% 
+  select('clan', 'year', 'id', 'rank')
+
+tblFemaleRanks$year  <- as.numeric(tblFemaleRanks$year)
+
+save(tblFemaleRanks, file = 'female_ranks.RData')
+
+
diff --git a/0_create_ranks/female_ranks.RData b/0_create_ranks/female_ranks.RData
diff --git a/0_hyena_data.RData b/0_hyena_data.RData
diff --git a/2_ars.R b/2_ars.R
@@ -12,7 +12,7 @@ library(here)
 library(rethinking)
 library(ggplot2)
 
-plot.dir <- '~/Dropbox/Documents/Research/Full_projects/2023 Inheritancy_mobility/plots/'
+plot.dir <- '~/../Dropbox/Documents/Research/Full_projects/2023 Inheritancy_mobility/plots/'
 load('0_hyena_data.RData')
 
 ars <- ars %>%
@@ -25,15 +25,15 @@ ars <- ars %>%
 ars$clan_num <- as.numeric(as.factor(ars$clan))
 
 ### Number of reproductive events experienced for guiding markov chain simulation
-svg(filename = paste0(plot.dir, '/repro_events_per_year.svg'), width = 6, height = 4)
+pdf(file = paste0(plot.dir, 'repro_events_per_year.pdf'), width = 6, height = 4)
 ars %>%
   group_by(id) %>%
   mutate(repro_scaled = repro_events/clan_size) %>%
   summarize(total_repro_scaled = sum(repro_scaled),
             repro_longevity = length(repro_events))%>%
   ggplot(aes(x = repro_longevity, y = total_repro_scaled))+
   geom_point() + 
-  geom_smooth(color = 'black', method = 'loess')+
+  geom_smooth(color = 'black', method = 'lm')+
   theme_classic()+
   ylab('Recruitment events experienced\n(scaled by clan size)')+
   xlab('Number of years observed reproducing')+
@@ -100,7 +100,7 @@ p.preds <- rbind(p.preds,
                             category = 'mean',
                             rank = 1:30))
 
-svg(paste0(plot.dir, 'prob_reproduce.svg'), width = 6, height = 4)
+pdf(paste0(plot.dir, 'prob_reproduce.pdf'), width = 4, height = 4)
 ggplot(data = filter(p.preds, category == 'post'), aes(x = rank, y = probs, group = sample))+
   geom_line(size = 0.1, alpha = 0.2)+
   geom_line(data = filter(p.preds, category == 'mean'), size = 2)+

diff --git a/2_repro_function.RData b/2_repro_function.RData
diff --git a/3_markov_model_rank_hyenas.R b/3_markov_model_rank_hyenas.R
@@ -18,7 +18,7 @@ library(expm)
 source('1_define_functions.R')
 set.seed(1989)
 
-plot.dir <- '~/Dropbox/Documents/Research/Full_projects/2023 Inheritancy_mobility/plots/'
+plot.dir <- '~/../Dropbox/Documents/Research/Full_projects/2023 Inheritancy_mobility/plots/'
 
 load('0_hyena_data.RData')
 
@@ -52,6 +52,8 @@ ranks.full.life <- ranks %>%
          last.year < last.year.clan)
 length(unique(ranks.full.life$id))
 
+quantile(table(ranks.full.life$id))
+
 dynamics <- ranks.full.life %>% 
   group_by(id) %>%
   summarize(total.passive = abs(sum(delta.passive.decile, na.rm = T)),
@@ -81,9 +83,10 @@ dynamics.plot <- ggplot(dynamics)+
   theme(axis.title.y = element_blank(), legend.position = 'none', 
         axis.text.y = element_text(hjust = 0.5, angle = 90))+
   xlab('Magnitude of net displacement from starting rank')+
-  scale_fill_manual(values = c('black', 'gold2'))+
+  scale_fill_manual(values = c('gray30', 'gray50'))+
   scale_y_discrete(expand = c(0.1,0))+
-  labs(tag = 'A')
+  labs(tag = 'A')+
+  ggtitle('Observed hierarchy dynamics')
 #dev.off()
 
 
@@ -103,21 +106,21 @@ for(tid in unique(ranks$id)){
   changes.passive <- ranks[ranks$id == tid,]$delta.passive.decile[-1]
   changes.active <- ranks[ranks$id == tid,]$delta.active.decile[-1]
 
-  rank.traj.passive <- start.rank + cumsum(changes.passive)
-  rank.traj.active <- start.rank + cumsum(changes.active)
-  rank.traj <- round(start.rank + cumsum(changes.passive) + cumsum(changes.active), 6)
+  rank.traj.passive <- c(start.rank, start.rank + cumsum(changes.passive))
+  rank.traj.active <- c(start.rank, start.rank + cumsum(changes.active))
+  rank.traj <- c(start.rank, round(start.rank + cumsum(changes.passive) + cumsum(changes.active), 6))
 
-  rank.state.list[[tid]] <-  cut(rank.traj, 
+  rank.state.list[[tid]] <-  cut(rank.traj,
                                  breaks = c(-0.001, 0.1, 0.2, 0.3, 0.4,
-                                            0.5, 0.6, 0.7, 0.8, 0.9, 1.001), 
+                                            0.5, 0.6, 0.7, 0.8, 0.9, 1.001),
                                  labels = states)
-  rank.state.passive.list[[tid]] <- cut(rank.traj.passive, 
+  rank.state.passive.list[[tid]] <- cut(rank.traj.passive,
                                         breaks = c(-0.001, 0.1, 0.2, 0.3, 0.4,
-                                                   0.5, 0.6, 0.7, 0.8, 0.9, 1.001), 
+                                                   0.5, 0.6, 0.7, 0.8, 0.9, 1.001),
                                         labels = states)
-  rank.state.active.list[[tid]] <- cut(rank.traj.active, 
+  rank.state.active.list[[tid]] <- cut(rank.traj.active,
                                        breaks = c(-0.001, 0.1, 0.2, 0.3, 0.4,
-                                                  0.5, 0.6, 0.7, 0.8, 0.9, 1.001), 
+                                                  0.5, 0.6, 0.7, 0.8, 0.9, 1.001),
                                        labels = states)
 
   cumulative.changes.list[[tid]] <- data.frame(id = tid, lifetime = length(rank.traj), 
@@ -132,12 +135,12 @@ cumulative.changes <- do.call(rbind, cumulative.changes.list)
 #### Fit markov models
 
 ## Total changes
-transitions.total <- markovchainFit(data = rank.state.list, method = 'map', confint = T)
+transitions.total <- markovchainFit(data = rank.state.list, method = 'map')
 transitions.total.estimate <- transitions.total$estimate@transitionMatrix
 limit.hyena.total <- data.frame(probs = find_limit(transitions.total.estimate),
                                 states)
 
-transitions.passive <- markovchainFit(data = rank.state.passive.list, method = 'map', confint = T)
+transitions.passive <- markovchainFit(data = rank.state.passive.list, method = 'map')
 transitions.passive.estimate <- transitions.passive$estimate@transitionMatrix
 limit.hyena.passive <- data.frame(probs = find_limit(transitions.passive.estimate),
                                 states)
@@ -199,15 +202,16 @@ hyena.total <- ggplot(data = total.predicted.ranks, aes(x = time, y = rank.mean,
   ggtitle('Total dynamics')
 
 hyena.total.limit <- ggplot(data = limit.hyena.total, 
-                             aes(y = states[10:1], x = probs, fill = states[10:1])) +
+                             aes(y = states[10:1], x = probs, fill = states[10:1], color = states[10:1])) +
   geom_bar(stat= 'identity') +
   scale_fill_manual(values = colorRampPalette(c('gold', 'black'))(10))+
+  scale_color_manual(values = colorRampPalette(c('gold', 'black'))(10))+
   theme_void()+
   theme(legend.position = 'none',
         axis.line.x = element_line(), axis.text.x = element_text(),
         axis.title.x = element_text())+
   scale_x_continuous(breaks = c(0, 0.5),
-                     limits = c(0, 0.5), name = 'Steady state')
+                     limits = c(0, 0.6), name = 'Steady state')
 #total
 #dev.off()
 
@@ -229,15 +233,16 @@ hyena.passive <- ggplot(data = passive.predicted.ranks, aes(x = time, y = rank.m
 
 
 hyena.passive.limit <- ggplot(data = limit.hyena.passive, 
-                            aes(y = states[10:1], x = probs, fill = states[10:1])) +
+                            aes(y = states[10:1], x = probs, fill = states[10:1], color = states[10:1])) +
   geom_bar(stat= 'identity') +
   scale_fill_manual(values = colorRampPalette(c('gold', 'black'))(10))+
+  scale_color_manual(values = colorRampPalette(c('gold', 'black'))(10))+
   theme_void()+
   theme(legend.position = 'none',
         axis.line.x = element_line(), axis.text.x = element_text(),
         axis.title.x = element_text())+
   scale_x_continuous(breaks = c(0, 0.5),
-                     limits = c(0, 0.51), name = 'Steady state')
+                     limits = c(0, 0.6), name = 'Steady state')
 #passive
 #dev.off()
 
@@ -258,19 +263,21 @@ hyena.active <- ggplot(data = active.predicted.ranks, aes(x = time, y = rank.mea
   ggtitle('Active dynamics')
 
 hyena.active.limit <- ggplot(data = limit.hyena.active, 
-                              aes(y = states[10:1], x = probs, fill = states[10:1])) +
+                              aes(y = states[10:1], x = probs, fill = states[10:1], color = states[10:1])) +
   geom_bar(stat= 'identity') +
   scale_fill_manual(values = colorRampPalette(c('gold', 'black'))(10))+
+  scale_color_manual(values = colorRampPalette(c('gold', 'black'))(10))+
   theme_void()+
   theme(legend.position = 'none',
         axis.line.x = element_line(), axis.text.x = element_text(),
         axis.title.x = element_text())+
   scale_x_continuous(breaks = c(0, 0.5),
-                     limits = c(0, 0.5), name = 'Steady state')
+                     limits = c(0, 0.6), name = 'Steady state')
+
 #active
 #dev.off()
 
-svg(paste0(plot.dir, '/hyena_dynamics.svg'), height = 4.8, width = 9)
+pdf(paste0(plot.dir, 'hyena_dynamics.pdf'), height = 4.8, width = 9)
 dynamics.plot + hyena.total + hyena.total.limit + 
   hyena.passive + hyena.passive.limit + 
   hyena.active + hyena.active.limit + 

diff --git a/5_markov_model_rank_simulated.R b/5_markov_model_rank_simulated.R
@@ -17,7 +17,7 @@ library(gridExtra)
 library(expm)
 
 rm(list = ls())
-plot.dir <- '~/Dropbox/Documents/Research/Full_projects/2023 Inheritancy_mobility/plots/'
+plot.dir <- '~/../Dropbox/Documents/Research/Full_projects/2023 Inheritancy_mobility/plots/'
 source('1_define_functions.R')
 load('4_rank_data_simulated.Rdata')
 load('2_repro_function.RData')

diff --git a/5_sim_markov_chains.Rdata b/5_sim_markov_chains.Rdata
diff --git a/6_assess_inequality.R → 6_sim_plot_queens_inequality.R b/6_assess_inequality.R → 6_sim_plot_queens_inequality.R
@@ -7,6 +7,7 @@ library(patchwork )
 library(gridExtra)
 library(ineq)
 library(tidyr)
+library(markovchain)
 
 rm(list = ls())
 source('1_define_functions.R')
@@ -15,7 +16,7 @@ load('2_repro_function.RData')
 load('5_sim_markov_chains.Rdata')
 set.seed(1989)
 
-plot.dir <- '~/Dropbox/Documents/Research/Full_projects/2023 Inheritancy_mobility/plots/'
+plot.dir <- '~/../Dropbox/Documents/Research/Full_projects/2023 Inheritancy_mobility/plots/'
 
 
 #### Queen effects
@@ -211,20 +212,14 @@ mocor.queen <- ggplot(data = filter(mocorqueen.predicted.ranks, start %in% queen
 
 mocor.queen
 
-
-# pdf('plot_grid_queens.pdf', width = 4, height = 4)
-# mri.youngest.queen + mri.youngest.norank.queen +
-#   none.queen + none.norank.queen
-# dev.off()
-
 mri.youngest + 
   geom_line(data = filter(mri.youngestqueen.predicted.ranks, start == queen.states[1]), col = 'coral', lty = 2, size = 1)
 
 none + 
   geom_line(data = filter(nonequeen.predicted.ranks, start == queen.states[1]), col = 'coral', lty = 2, size = 1)
 
 
-pdf(paste(plot.dir, 'simulations.pdf'), width = 8, height = 4)
+pdf(paste0(plot.dir, 'simulations.pdf'), width = 8, height = 4)
 mri.youngest + geom_line(data = filter(mri.youngestqueen.predicted.ranks, start == queen.states[1]), col = 'coral', lty = 2, size = 1) + 
 mri.youngest.limit +
   mri.youngest.norank + geom_line(data = filter(mri.youngest.norankqueen.predicted.ranks, start == queen.states[1]), col = 'coral', lty = 2, size = 1) +  mri.youngest.norank.limit +
@@ -301,7 +296,7 @@ reproduction <- ggplot(repro.plot, aes(x = condition, y = gini.u))+
 #dev.off()
 
 
-svg(paste0(plot.dir, 'reproductive_inequality.svg'), width = 9, height = 5)       
+pdf(paste0(plot.dir, 'reproductive_inequality.pdf'), width = 9, height = 5)       
 reproduction + 
   inset_element(none.queen, 0.08, 0.6, 0.30, 0.9, align_to = 'full') +
   inset_element(mri.youngest.queen, 0.30, 0.6, 0.52, 0.9, align_to = 'full')+