################################## ### Electronic Supplementary Material: Code ### "The effect of dispersal on rates of cumulative cultural evolution" ### Author: Mark Dyble ################################## library("stringr"); library("QRM"); library("ggplot2"); library("kinship2") # Load the required packages (install if required) is.even <- function(x) x %% 2 == 0 ############################### #### Key model parameters ############################## ### Demographic G <- 20 # How many groups? N <- 24 # How many agents per group? (Groups are two generations of N/2) Popn_gen_group <- N/2 # How many individuals per group? (Must be multiple of 4). NB: Actual group size will be double this as we keep two generations alive at any one time P <- G*Popn_gen_group ## Overall size of a generation ### Culutral alpha <- 3 ## Erosian parameter beta <- 1.5 ## Distribution parameter Starting_z <- 10 DispersalSystem <- "Bilocal" ## Can be either "Bilocal"or "Patrilocal" ### Simulations Gens <- 100 # How many generations to run the model for? runs <- 5 ## How many runs? ############################### #### Start of model ############################## for(run in 1:runs){ mean_z_pat <- mean_z_mat <- mean_z <- rep(NA,times=Gens) # For storing results for(generation in 1:Gens){ ########################### ## Initialise new generation ########################### Popn <- as.data.frame(matrix(0,P,8)); colnames(Popn) <- c("ID","Mother","Father","Sex","Group","Zi","Mat_Z","Pat_Z") # Dataframe to store agent 'data' Popn$ID <- 1:P # Assign ID numbers Popn$Group <- rep(1:G,each=Popn_gen_group) # Assign natal groups Popn$Sex <- rep(c(1,2),times=(0.5*nrow(Popn))) ## Assign gender (male = 1) if(generation==1){ Popn[,2:3] <- NA # No parents in generation 1 Popn[,6:8] <- Starting_z # Give default starting skill level in generation 1 } #################### ## Parentage and inheritance of skills #################### if(generation>1){ res_mo <- res_fa <- res_bi <- res_maleZ <- res_femaleZ <- 1:P # Data store Parents_Name <- Parents_gen$ID; Parents_Groups <- Parents_gen$Group;Popn_Groups <- Popn$Group; ParentsZi <- Parents_gen$Zi;ParentsPatz <- Parents_gen$Pat_Z;ParentsMatz <- Parents_gen$Mat_Z # Pull out data used to save time in loop with lookups for(i in 1:P){ # Give each a random z from one male and one female in their group with error in inheritance ### Assign parents to new individuals same_g <- which(Parents_Groups==Popn_Groups[i]) # Which agents from the previous generation are in this new agents group? mother <- sample(same_g[is.even(same_g)],1) # Choose a random female from the above to be the mother father <- sample(same_g[!is.even(same_g)],1) # Choose a random male from the above to be the father res_mo[i] <- Parents_Name[mother]; res_fa[i] <- Parents_Name[father] # Store the data ### New agents vertically interit from parents res_bi[i] <- rGumbel(1, mu = ((max(c(ParentsZi[mother],ParentsZi[father])))-alpha), sigma = beta) ## Bilineal inheritance sampled from the Gumbel distribution. res_maleZ[i] <- rGumbel(1, mu = (ParentsPatz[father]-alpha), sigma= beta) ## Patrilineal inheritance sampled from the Gumbel distribution. res_femaleZ[i] <- rGumbel(1, mu = (ParentsMatz[mother]-alpha), sigma= beta) ## Matrilineal inheritance sampled from the Gumbel distribution. } Popn$Zi <- res_bi; Popn$Mat_Z <- res_femaleZ; Popn$Pat_Z <- res_maleZ; Popn$Mother <- res_mo; Popn$Father <- res_fa # Store the results in the main dataframe (Popn) Popn$Zi[which(Popn$Zi<=0)] <- 0; Popn$Mat_Z[which(Popn$Mat_Z<=0)] <- 0; Popn$Pat_Z[which(Popn$Pat_Z<=0)] <- 0 ## Make sure no skills level fall below 0 } #################### #### Dispersal #################### if(DispersalSystem=="Patrilocal"){ Movers_women <- which(Popn$Sex==2) # All the females move Popn$Group[Movers_women] <- sample(Popn$Group[Movers_women]) # Shuffle the movers between camps randomly } if(DispersalSystem=="Bilocal"){ Movers <- rep(((Popn_gen_group/2)+1):Popn_gen_group,times=G) + rep((0:(G-1))*Popn_gen_group,each=(Popn_gen_group/2)) ## Half the group disperse. Equal number of Males and Females MoversMale <- Movers[which(Popn$Sex[Movers]==1)] MoversFemale <- Movers[which(Popn$Sex[Movers]==2)] Popn$Group[MoversMale] <- sample(Popn$Group[MoversMale]) # Shuffle the dispersing males between camps randomly Popn$Group[MoversFemale] <- sample(Popn$Group[MoversFemale]) # Shuffle the dispersing females between camps randomly } ####################################### ## Inherit skills horizontally from group members ####################################### Popn$Pat_Z[which(is.even(1:P))] <- NA ; Popn$Mat_Z[which(!is.even(1:P))] <- NA; Group <- Popn$Group # Pull out data to save time Res <- Skills <- Popn$Zi; ResP <- PatSkills <- Popn$Pat_Z; ResM <- MatSkills <- Popn$Mat_Z # For saving data for(i in 1:nrow(Popn)){ # Iterate through all individuals own_pat_skill <- PatSkills[i] ; own_mat_skill <- MatSkills[i] ; own_skill <-Skills[i] ## What are their skills levels? names_group <- which(Group==Group[i]) # Which inds are in their group? skills_group <- Skills[names_group[which(Skills[names_group]>own_skill)]] # Skill levels of those in their group skills_male <- PatSkills[names_group[!is.even(names_group)][which(PatSkills[names_group[!is.even(names_group)]]>own_pat_skill)]] # Skill levels of those in their group skills_female <- MatSkills[names_group[is.even(names_group)][which(MatSkills[names_group[is.even(names_group)]]>own_mat_skill)]] # Skill levels of those in their group if(length(skills_group)>1) Res[i] <- sample(skills_group,1,prob = (skills_group-own_skill)) else if(length(skills_group)==1) Res[i] <- skills_group ## Horizontal transmission from more skilled agents if(!is.even(i)){if(length(skills_male)>1) ResP[i] <- sample(skills_male,1,prob = (skills_male-own_pat_skill)) else if(length(skills_male)==1) ResP[i] <- skills_male} ## Horizontal transmission from more skilled agents if(is.even(i)){if(length(skills_female)>1) ResM[i] <- sample(skills_female,1,prob = (skills_female-own_mat_skill)) else if(length(skills_female)==1) ResM[i] <- skills_female} ## Horizontal transmission from more skilled agents } Popn$Zi <- Res; Popn$Mat_Z <- ResM; Popn$Pat_Z <- ResP ## Update z values mean_z[generation] <- mean(Popn$Zi) # Record the mean z values for the population mean_z_mat[generation] <- mean(Popn$Mat_Z[which(Popn$Sex==2)]) mean_z_pat[generation] <- mean(Popn$Pat_Z[which(Popn$Sex==1)]) #################### ## Housekeeping #################### ## The below adds the latest generation to a list of all individuals (All_Inds). To avoid repeating ID numbers, move all exsiting IDs up by P if(generation==1) All_Inds <- Popn if(generation>1) All_Inds <- rbind(Popn,All_Inds) All_Inds$ID <- All_Inds$ID+P; All_Inds$Mother <- All_Inds$Mother+P; All_Inds$Father <- All_Inds$Father+P Parents_gen <- All_Inds[1:P,] } ###################### ## End of iteration ###################### ### Store results in 'results' dataframe: if(run==1){ results <- as.data.frame(matrix(0,Gens*3,3)); colnames(results) <- c("Inheritance","Gen","Z_score") results$Inheritance <- rep(c("Bilineal","Matrilineal","Patrilineal"),each=Gens) results$Gen <- as.numeric(rep(1:Gens,times=3)) results$Z_score <- as.numeric(as.character(c(mean_z,mean_z_mat,mean_z_pat))) } if(run>1){ temp <- as.data.frame(matrix(0,Gens*3,3)); colnames(temp) <- c("Inheritance","Gen","Z_score") temp$Inheritance <- rep(c("Bilineal","Matrilineal","Patrilineal"),each=Gens) temp$Gen <- as.numeric(rep(1:Gens,times=3)) temp$Z_score <- as.numeric(as.character(c(mean_z,mean_z_mat,mean_z_pat))) results <- rbind(temp,results) } } ################################## # End of model ################################## ################################## ### Visulaise results results <- results[which(results$Gen%in%((0:(Gens/5))*5)),] results_agg <- aggregate(results,by=list(results$Inheritance,results$Gen),FUN=mean) Gen_labels <- unique(results_agg$Gen) plot(Gen_labels,results_agg$Z_score[which(results_agg$Group.1=="Bilineal")],ylab="Skill level (mean z)",bty='l',ylim = c(0,max(results_agg$Z_score)),xlab = "Generation",pch=0,type="o",lwd=1.2) points(Gen_labels,results_agg$Z_score[which(results_agg$Group.1=="Matrilineal")],pch=1,type="o",lwd=1.2,col="red") points(Gen_labels,results_agg$Z_score[which(results_agg$Group.1=="Patrilineal")],pch=2,type="o",lwd=1.2,col="blue") ##################################