#====================================================== # Code for generating community composition # Authors of the original model: Mark Dyble and Gul Deniz Salali # contact: mark.dyble@cantab.net # Last updated by author 26/06/2017 # Note - this is a modified version which is much faster than the one used for the original analysis and also produces the pie charts in R # You are welcome to use this code freely but please cite the following paper if you do so: # M. Dyble et al (2015) "Sex equality can explain the unique social structure of hunter-gatherer bands". Science 348.6236: 796-798. # See the Supplimentary data for additional description of the model #===================================================== library(igraph) # igraph package is required ## Library kin_names <- c("father","mother","son","daughter","brother","sister") bilocal <- c(.125,.125,.125,.125,.25,.25) # Probability of drawing different kin patrilocal <- c(.25,0,.25,0,.5,0)# Probability of drawing different kin dispersal <- "bilocal" # can be patrilocal or bilocal popsize= 20 # Must be even. Popn size of 20 used in Fig1 in the paper runs <- 100 # How many simulations? normal_master= matrix(0,10,runs); normal_master=as.data.frame(normal_master); colnames(normal_master)=c(1:runs) all_results <- matrix(0,9,runs); rownames(all_results) <- c("S","PK","DK","SP","SPK","PKS","SDK","OA","NR") for(run in 1:runs){ #=========================== # MAKING OUR FICTIONAL COMMUNITY #=========================== d=matrix(0,2,5) # create an initial dataframe with 2 rows for the first two people in camp colnames(d)=c("id","fatherid","motherid", "spid","sex") d=as.data.frame(d) ### Populate the initial camp with two people d[,1]= c(1,2) # ID names of original two camp inhabitants d[1,2:3]=c(101,102) # Parents of the first person d[2,2:3]=c(103,104) # Parents of the second person d[1:2,4]=c(2,1) # Names of the spouses d[1:2,5]=c(1,2) # Sex (male = 1, fem = 2) n=2; new= 3; new.par=105 # Set up before the main loop # this loop is the core of the model - it builds us a fictional camp while(n <= popsize-2){ # Draw an existing agent to choose a kin memeber to join the camp if(dispersal=="bilocal") agent= sample(d[,1],1) # choose an existing agent at random if(dispersal=="patrilocal") agent= sample(d[(which(d$sex==1)),1],1) ## Draw a kin member of this agent to join the camp agents_father <- d$fatherid[which(d$id==agent)] # see below - this is part of avoiding duplicating people. i.e. if your brother already 'brought' your father, you can't bring him again if(dispersal=="bilocal") if(!agents_father%in%d$id) kin <- sample(kin_names,1,prob = bilocal) else kin <- sample(kin_names[3:6],1,prob = bilocal[3:6]) if(dispersal=="patrilocal") if(!agents_father%in%d$id) kin <- sample(kin_names,1,prob = patrilocal) else kin <- sample(kin_names[3:6],1,prob = patrilocal[3:6]) if(kin=="brother"){ row <- row2 <- rep(0,ncol(d)) # new agent and their spouse to bind to dataframe row[1:5]= c(new,d$fatherid[which(d$id==agent)],d$motherid[which(d$id==agent)],new+1,1) row2[1:5] <- c(new+1,new.par,new.par+1,new,2) new= new + 2; new.par=new.par+2 d=rbind(d,row,row2) # add the agent to the dataframe } else if(kin=="sister"){ row <- row2 <- rep(0,ncol(d)) # new agent and their spouse to bind to dataframe row[1:5]= c(new,d$fatherid[which(d$id==agent)],d$motherid[which(d$id==agent)],new+1,2) row2[1:5] <- c(new+1,new.par,new.par+1,new,1) new= new + 2; new.par=new.par+2 d=rbind(d,row,row2) # add the agent to the dataframe } else if(kin=="father"){ row <- row2 <- rep(0,ncol(d)) # new agent and their spouse to bind to dataframe row[1:5] <- c(d$fatherid[which(d$id==agent)],new.par,(new.par+1),d$motherid[which(d$id==agent)],1) # new parent info row2[1:5] <- c(d$motherid[which(d$id==agent)],new.par+2,new.par+3,d$fatherid[which(d$id==agent)],2) new.par <- new.par + 4 # add 1 everytime you use an id no d=rbind(d,row,row2) # add the two new agents to the dataframe } else if(kin=="mother"){ row <- row2 <- rep(0,ncol(d)) # new agent and their spouse to bind to dataframe row[1:5] <- c(d$motherid[which(d$id==agent)],new.par,(new.par+1),d$fatherid[which(d$id==agent)],2) # new parent info row2[1:5] <- c(d$fatherid[which(d$id==agent)],new.par+2,new.par+3,d$motherid[which(d$id==agent)],1) new.par <- new.par + 4 # add 1 everytime you use an id no d=rbind(d,row,row2) # add the two new agents to the dataframe } else if(kin=="son"){ row <- row2 <- rep(0,ncol(d)) # new agent and their spouse to bind to dataframe row[c(1,4,5)] <- c(new,new+1,1) if(d$sex[which(d$id==agent)]==1) { # if agent is male row[2:3]= c(agent,d$spid[which(d$id==agent)]) } else { # if agent is female row[2:3]= c(d$spid[which(d$id==agent)],agent) } row2[1:5]= c(new+1,new.par,new.par+1,new,2) # details of the new sppouse new.par= new.par + 2; new=new+2 d=rbind(d,row,row2) # add the agent to the dataframe } else if(kin=="daughter"){ row <- row2 <- rep(0,ncol(d)) # new agent and their spouse to bind to dataframe row[c(1,4,5)] <- c(new,new+1,2) if(d$sex[which(d$id==agent)]==1) { # if agent is male row[2:3]= c(agent,d$spid[which(d$id==agent)]) } else { # if agent is female row[2:3]= c(d$spid[which(d$id==agent)],agent) } row2[1:5]= c(new+1,new.par,new.par+1,new,1) # details of the new sppouse new.par= new.par + 2;new=new+2 d=rbind(d,row,row2) # add the agent to the dataframe } n=n+2 } ################################### # Code relatedness between our agents ################################### pk_mat <- matrix(0,popsize,popsize) # This will be a matrix with a '1' between all sibling and parent-child dyads for(i in 1:popsize){ for(j in 1:popsize){ if(d$fatherid[i]==d$fatherid[j]) pk_mat[i,j] <- 1 # if they share the same father, they are sibs if(d$fatherid[i]==d$id[j]) pk_mat[i,j] <- 1 # father of 1? if(d$motherid[i]==d$id[j]) pk_mat[i,j] <- 1 # mother of 1? if(d$id[i]==d$motherid[j]) pk_mat[i,j] <- 1 # i is mother of j? if(d$id[i]==d$fatherid[j]) pk_mat[i,j] <- 1 # i is father of j? } } cat_mat <- mar_mat <- matrix(0,popsize,popsize);rownames(cat_mat) <-rownames(cat_mat) <-rownames(mar_mat) <-colnames(mar_mat) <- rownames(pk_mat) <-colnames(pk_mat) <- d$id for(i in 1:popsize){ # This mar_mat will have a '1' between all spouses ego <- rownames(mar_mat)[i] alter <- d$spid[which(d$id==ego)] mar_mat[i,as.character(alter)] <- 1 } both_mat <- mar_mat + pk_mat # combine the two # construct networks - this is what we need igraph for pk_net <- graph.adjacency(pk_mat,weighted=TRUE, mode="undirected") mar_net <- graph.adjacency(mar_mat,weighted=TRUE, mode="undirected") both_net <- graph.adjacency(both_mat,weighted=TRUE, mode="undirected") ## shortest matrices too pk.short.net<- shortest.paths(pk_net, v=V(pk_net), to=V(pk_net),mode = c("all", "out", "in"),weights = NA, algorithm = c("automatic", "unweighted","dijkstra", "bellman-ford","johnson")) mar.short.net<- shortest.paths(mar_net, v=V(mar_net), to=V(mar_net),mode = c("all", "out", "in"),weights = NA, algorithm = c("automatic", "unweighted","dijkstra", "bellman-ford","johnson")) both.short.net<- shortest.paths(both_net, v=V(both_net), to=V(both_net),mode = c("all", "out", "in"),weights = NA, algorithm = c("automatic", "unweighted","dijkstra", "bellman-ford","johnson")) # to get the spouse info, we need to swap lines. So 1,3,5 become 2,4,6 and vice versa. spouse_kin_mat <- matrix(0,popsize,popsize) evens <- (1:(popsize/2))*2; odds <- ((1:(popsize/2))*2)-1 spouse_kin_mat[odds,] <- pk.short.net[evens,] spouse_kin_mat[evens,] <- pk.short.net[odds,] # The order we do these in is important cat_mat[which(spouse_kin_mat==(both.short.net-1))] <- "SDK" # Spouse's distant kin, will inc other that are replaced below cat_mat[which(both.short.net==2)] <- "PKS" # to start with this will inc SPK and DK ties too but they are replaced below cat_mat[which(both.short.net==2 & pk.short.net>2 & spouse_kin_mat == 1)] <- "SPK" cat_mat[which(pk.short.net<6 & pk.short.net==both.short.net)] <- "DK" cat_mat[which(pk_mat==1)] <- "PK" ## Primary Kin for(i in 1:popsize){ cat_mat[i,i] <- "self"} ## Code self cat_mat[which(mar_mat==1)] <- "SP" ## Spouse cat_mat[which(both.short.net>5)] <- "NR" # Unrelated people cat_mat[which(cat_mat==0)] <- "OA" # The remainder are other affines #### Count Up all_results[1,run] <- length(which(cat_mat=="self")) all_results[2,run] <- length(which(cat_mat=="PK")) all_results[3,run] <- length(which(cat_mat=="DK")) all_results[4,run] <- length(which(cat_mat=="SP")) all_results[5,run] <- length(which(cat_mat=="SPK")) all_results[6,run] <- length(which(cat_mat=="PKS")) all_results[7,run] <- length(which(cat_mat=="SDK")) all_results[8,run] <- length(which(cat_mat=="OA")) all_results[9,run] <- length(which(cat_mat=="NR")) } ## Visualise mean_res <- rowMeans(all_results) pie(mean_res, main="Categorical relatedness",clockwise=T,col=c("grey30","grey30","grey30","grey90","grey90","grey90","grey90","white","white"))