--- title: "fishSim-vignette" author: "Shane M Baylis" date: "`r Sys.Date()`" # output: rmarkdown::html_vignette output: rmarkdown::pdf_document: number_sections: true vignette: > %\VignetteIndexEntry{fishSim-vignette} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) options(mc.cores = 2) ``` fishSim is a package for demographic simulation and exploration of kin relationships. It includes functions to build populations and move, mate, kill, age, and switch the sex of individuals within those populations or defined subpopulations. It _doesn't_ include any one-line functions that define and run an entire simulation. Instead, users will have to build their own for-loop over years, with each year containing mating, mortality, ageing, sampling, and/or movement events as needed. With a little planning, this makes fishSim _very_ adaptable. Importantly, fishSim is individual-based and retains parentage information for all individuals, allowing full reconstruction of kin relationships to arbitrary depth. In addition to population simulation features, fishSim includes a set of convenience functions allowing the user to check the scenario population growth rate, find simulation scenarios with null population growth (`PoNG()`), selectively archive subsets of the population for increased computational speed. For a given population, fishSim also includes functions to look up pairs of animals, compare their sets of ancestors, and classify them into kin categories based on their nearest shared ancestor(s), or show patterns of shared ancestors for members of a pair. fishSim is developed in response to the demands of close-kin mark-recapture projects. 'typical' use-cases therefore involve an understanding of both demography and kin relationships. This vignette will cover two use scenarios - one relatively simple one, with: - age-specific (but not sex-specific) fecundity, - constant adult survival, - a single population with no subdivisions (i.e., one 'stock'), - no sex-switching, and - first-year survival set by `PoNG()`, giving a flat population size through time, - one-shot sampling of animals from the population ...and a one designed to show off the bells and whistles, with: - age-specific, sex-specific fecundity with within-season mating exhaustion for females, - age-specific, sex-specific survival, - three subpopulations ('stocks') with markovian movement between subpopulations, - male-to-female sex-switching, - a growing population, - lethal sampling, occurring through time during the simulation. In the 'bells-and-whistles' scenario, we will also cover archiving tools. For both scenarios, we will take a brief look at relationship patterns between pairs of animals. \newpage # The relatively-simple scenario ## Setup In order to simulate a population, we need an object to hold that population's data, and the population will need some founding members. Our data object will have one line per individual, so let's call it `indiv`. ```{r} library(fishSim) indiv <- makeFounders(stocks = c(1)) head(indiv) ``` Now, what's here? `indiv` has one line per individual, and each individual has an ID code (`indiv$Me`), a sex (`indiv$Sex`), a father's ID code (`indiv$Dad` - but which is 'founder' here, because founding animals can't really be said to have parents), a mother's ID code (`indiv$Mum` - ditto), a birth year `indiv$BirthY`, a death year (`indiv$DeathY` - NA for any animals that are still alive), a stock membership (`indiv$Stock`), an age in years (`indiv$AgeLast`), and a flag indicating whether the individual has been sampled (`indiv$SampY`). The age structure, sex ratio, and stock memberships of founders can all be specified in the `makeFounders()` call. We've only changed one of the default values - instead of the default 3 stocks, we're going to have a single panmictic founder population. The other defaults give us 1000 founders, with an even sex ratio, with an age-structure that implies 70% annual survival, but with a hard age-limit of 20 years. ## Mating Having set up our population, we can now get its members to breed, get older, switch sex and die, using `mate()` (or `altMate()` or `bondedMate()`), `birthdays()`, `sexSwitch()`, and `mort()` functions. If we had more than one stock, we could also get individuals to move between stocks with markovian movement probabilities, using `move()`. Let's look at those. ```{r} nrow(indiv) indiv <- mate(indiv, year = 1) nrow(indiv) ## 200 newborns tail(indiv) ## newborns are added to the end ``` Here are some newborns, generated using the default `mate()` settings. Unlike founders, these newborns all have a father's ID and a mother's ID, are all aged 0, and all have a birth year of equal to the `year` argument. There are 200 newborns, because `mate()` generates new members as a proportion of the extant population size (i.e., we specify fecundity for the population, not the individuals), and that proportion (argument `fecundity`) is 0.2 by default. Each individual was born in a clutch, and each clutch contained a Poisson-distributed number of clutch-members, the default mean ('batchSize') of 0.5. `mate()` simply kept producing newborns by mating together random male/female pairs until `nrow(indiv) * 0.2` new offspring were produced. Alternatively, we can mate our population using `bondedMate()`, which stores a mate ID for mature animals and generates offspring within each pair unless one member dies or moves to another stock, or the pair divorces. Re-mating occurs among surviving animals, but this system can give repeated breeding between pair members for many breeding seasons. Finally, we can mate our population using `altMate()`, which specifies individual maturities by age, and a probability distribution of number of offspring per mature female. Let's generate some new founders and try that. ```{r} indiv <- makeFounders(stocks = c(1)) nrow(indiv) ## 1000 founders indiv <- altMate(indiv, firstBreed = 2, batchSize = 0.9, year = 1) tail(indiv) ``` In this call to `altMate()`, females became sexually mature at 2 years of age (`firstBreed = 2`), and each sexually-mature female produced a Poisson-distributed number of offspring with mean set by `batchSize`. It is possible to specify age-specific or sex-specific fecundities, set paternity within-clutch to be single or multiple, to exhaust fathers within each breeding season (females, in this system, always breed to exhaustion each season), and to set a sex ratio for offspring by additional arguments to `altMate()`, and these options can also be set in `mate()` and `bondedMate()`. ## Mortality Now, let's kill some of the population. Mortality probabilities can be flat, age-specific or stock-specific (and with a little extra effort, sex- or age:stock:sex-specific too), or we can randomly kill animals such that the population is reduced to a certain size, and we can set an age past which no animal will survive. For now, let's set a flat 20% mortality rate (i.e., the probability of death is equal for all animals). ```{r} nrow(indiv) indiv <- mort(indiv, year = 1, type = "flat", mortRate = 0.2) nrow(indiv) head(indiv, n = 15) ``` Mortality does not remove dead animals from 'indiv' - we will need to refer to them later - it just updates their death year to the value given in `year`. In general, though, animals with a non-NA value for death year will not move, mate, switch sex, or have birthdays. ## Birthdays On the topic of birthdays: that's the final thing that must be done before we can put this all together and run a full (albeit basic) demographic simulation: we need to be able to increase the age of our animals, and that is what the `birthdays()` function is for. It's very simple: it takes all living members of the population, and increments their age by 1. ```{r} tail(indiv) indiv <- birthdays(indiv) tail(indiv) ``` ## Sampling from the population It is often useful to keep track of which animals in the population have been captured, and if capture is lethal, to mark captured animals as dead through a separate process from the normal 'mort' functions. The function `capture()` takes a sample of individuals, marks them as 'captured' by updating their 'captured' value, and, if capture is lethal, marks them as dead by updating their death year, as in `mort()`. It is possible to make `capture()` sex-specific, age-specific, or age:sex specific, so that only one sex may be captured, only one age-class may be captured, or only one sex at one age-class may be captured, respectively. ```{r} ## non-lethal sampling of 5 females indiv <- capture(indiv, n = 5, year = 1, fatal = FALSE, sex = "F") ## lethal sampling of 5 animals, from either sex indiv <- capture(indiv, n = 5, year = 1, fatal = TRUE) ## non-lethal sampling of 5 one-year-old females indiv <- capture(indiv, n = 5, year = 1, fatal = FALSE, sex = "F", age = 1) ``` ## What's my population doing? Can I make it stay the same size? A couple more functions deserve a mention here: `check_growthrate()` and `PoNG()`. `check_growthrate()` tells you how quickly your population is growing. `PoNG()` tells you what you need first-year survival to be, in order for your population to maintain a constant size (within limits - the estimation is Leslie Matrix-based, and there are ways to fool Leslie Matrices that are out of scope for this vignette but covered briefly in the `PoNG()` and `check_growthrate()` documentation). Let's try both of those, using the mating and mortality settings we have already used. ```{r} check_growthrate(mateType = "flat", mortType = "flat", batchSize = 0.9, firstBreed = 2, mortRate = 0.2) ## with the current settings, we expect our population to grow by about 1.8% annually. PoNG(mateType = "flat", mortType = "flat", batchSize = 0.9, firstBreed = 2, mortRate = 0.2) ## if first-year mortality was about 0.445, rather than 0.2, our population ## would have null growth. You can also read off a range of possible growth rates ## from the plot. ``` ## Turning processes into a simple simulation To turn those processes into a full demographic simulation, all that is needed is to repeat the processes in a loop. Let's do that, but set age-specific mortality so that the long-run average population growth rate is zero using the first-year mortality rate we got from `PoNG()`, so that our population is unlikely to explode or become extinct. ```{r} indiv <- makeFounders(stocks = c(1)) ageMort <- c(0.445, rep(0.2, 100)) ## age-specific mortality is 0.445 for first-years, ## 0.2 for all older age-classes. for (y in 1:60) { indiv <- altMate(indiv, firstBreed = 2, batchSize = 0.9, year = y) ## y for year indiv <- mort(indiv, year = y, type = "age", ageMort = ageMort) ## age-specific mort indiv <- birthdays(indiv) } tail(indiv) ## a population with 60 years of births, deaths, and birthdays nrow(indiv[is.na(indiv$DeathY),]) ## the currently-alive population size. Note that population ## growth only *averages* zero, and variability occurs! ``` ## Looking up relationships between pairs of animals One of the key advantages of fishSim is its ability to report on different relationships between pairs of animals in a sample. Are these two each other's siblings? Half-siblings? Is one the parent of the other? The grandparent? Some obscure half-cousin, once removed? Internally, these relationships are stored in terms of 'shared ancestors at the n'th generation of each member', starting with the 'self' as generation 1. So a pair might have a shared `pairs$TwoThree` ancestor - that ancestor is the parent of one member, and the grandparent of the other. Or the pair might have a shared `pairs$OneTwo` ancestor - in which case, one member's _self_ is the other member's parent. If a pair shares a parent (i.e., `pairs$TwoTwo == 1`), then they also share at least two grandparents (`pairs$ThreeThree >= 2`), four great-grandparents (`pairs$FourFour >= 4`), and so on. This pattern holds generally: a half thiatic pair (`pairs$TwoThree == 1`) will also share at least two `pairs$ThreeFour >= 2` ancestors and at least four `pairs$FourFive >= 4` ancestors. It is also possible, for example, for a pair to share two grandparents, but no parents - that is the case for full cousin pairs. ```{r} ## mark 200 individuals alive at the end of the simulation as 'captured'. indiv <- capture(indiv, n = 200, year = 60, fatal = FALSE) ## look up each animal's ancestors and look for shared ancestors between each ## pair of sampled animals: pairs <- findRelativesPar(indiv = indiv, sampled = TRUE, nCores = 2) POPs <- pairs[pairs$OneTwo == 1,] ## Parent-Offspring pairs GGPs <- pairs[pairs$OneThree == 1,] ## Grandparent-Grandoffspring pairs HSPs <- pairs[pairs$TwoTwo == 1,] ## Half-sibling pairs FSPs <- pairs[pairs$TwoTwo == 2,] ## Full Sibling pairs (self-comparisons ## are automatically excluded) FCPs <- pairs[pairs$ThreeThree == 2 & pairs$TwoTwo != 1,] ## Full Cousin pairs ## look at the number of shared ancestors at each ancestral ## generation, for one of the half-sibling pairs. lookAtPair(HSPs[1,]) relatives <- namedRelatives(pairs) ## shows the number of pairs of each relationship type relatives ``` \newpage # The 'bells and whistles' scenario In this second scenario, we will simulate a metapopulation with: - age-specific, sex-specific fecundity with within-season mating exhaustion for females, - stock-specific, age-specific survival, - three subpopulations ('stocks') with markovian movement between subpopulations, - male-to-female sex-switching, - a growing population. - archiving for speed (?) We set up the population basically as before, with a couple of slight tweaks. First, we should make sure the founding population has multiple stocks (we'll give it three). Second, we will set up an archive matrix. The archive matrix is intended to hold the records of dead animals in large, long-running simulations: dead animals do not take part in any further mating, movement, sex-switching, aging, or mortality events, but if they stay in `indiv`, `indiv` can become huge and unweildy, slowing down all of those processes. There is of course a trade-off, in that writing dead animals to the archive takes system time, so the optimum may be to only archive dead animals once every few 'years'. ## Setting up the data objects ```{r} ## set up founders with three stocks: two that each contain 30% of the population, ## and one that contains the remaining 40%. indiv <- makeFounders(pop = 1000, stocks = c(0.3, 0.3, 0.4)) ## set up archive - just a matrix with zero rows and eight columns archive <- make_archive() ``` ## Parameterising movement, survival, and maturity structures Because we will have inter-stock movement in this sim, we will also need to set up a matrix giving the probability that an individual will move into another stock, given its current stock membership. Because survival will be age- and stock-dependent, we will need a matrix of survival rates with as many columns as stocks, and as many rows as (at least) the age of the oldest plausible population-member. Because we will have age-specific, sex-specific fecundity, we will need a male age-specific maturity curve and a female age-specific maturity curve. Let's set those up now. ```{r} ## Markovian movement matrix stocks <- c(0.3, 0.3, 0.4) admix.m <- matrix(NA, nrow = length(stocks), ncol = length(stocks)) for(i in 1:nrow(admix.m)) { admix.m[i,] <- stocks*stocks[i] } ## admix.m shows movement proportional to starting population sizes. admix.m ## let's tweak those numbers so that animals tend to stay where ## they are, and not move around so much. admix.m <- matrix(c(0.23, 0.03, 0.04, 0.03, 0.23, 0.04, 0.04, 0.04, 0.32), nrow = length(stocks), ncol = length(stocks), byrow = FALSE) admix.m ## Age- and stock-dependent survival ageStockMort <- matrix(c(0.47, 0.37, 0.27, rep(0.23, 97), 0.45, 0.35, 0.25, 0.20, rep(0.22, 96), 0.45, 0.3, 0.3, 0.19, rep(0.2, 96)), ncol = length(stocks), nrow = 100) head(ageStockMort) ## Sex-specific maturity curves maleCurve <- c(0,0,0,0.1,0.5,0.8,0.85,0.9,0.95,rep(1, 91)) femaleCurve <- c(0,0,0.5,0.9,0.95,rep(1,95)) ## maleCurve and femaleCurve should both be long enough that no individuals ## will outlive the curve. head(maleCurve) head(femaleCurve) ``` ## Checking population growth rates with the new parameters Now, what's going to happen with these subpopulations, if we run with those parameters? ```{r} check_growthrate(mateType = "ageSex", mortType = "ageStock", batchSize = 1.6, femaleCurve = femaleCurve, ageStockMort = ageStockMort) ## Not bad. Two of the three populations are increasing. The third will probably be kept viable ## by immigration from the other two. Note that I fiddled 'batchSize' (which is the mean number ## of offspring per mature female per breeding attempt) a bit. ``` ## Running the simulation Now we just run the population simulator, as we did for the simple population, but with a couple of extra things to specify. Let's give it 100 years. ```{r} for (k in 1:100) { ## very rarely, switch some males to females indiv <- sexSwitch(indiv = indiv, direction = "MF", prob = 1e-04) ## move animals according to the markovian matrix we set up before indiv <- move(indiv = indiv, moveMat = admix.m) ## mate animals using the age-specific, sex-specific curves we set up before indiv <- altMate(indiv = indiv, batchSize = 1.6, type = "ageSex", maleCurve = maleCurve, femaleCurve = femaleCurve, year = k) ## kill animals on the basis of their age and stock, as set up before indiv <- mort(indiv = indiv, year = k, type = "ageStock", ageStockMort = ageStockMort) if(k %% 10 == 0) { archive <- archive_dead(indiv = indiv, archive = archive) indiv <- remove_dead(indiv = indiv) cat( sprintf( '\rIteration %i complete', k)) ## occasionally sample from the population and do some clean-up. And report on progress. } if(k %in% c(94:100)) { indiv <- capture(indiv, n = 40, fatal = TRUE, year = k) } indiv <- birthdays(indiv = indiv) } archive <- rbind(archive, indiv) indiv <- archive ## merge 'indiv' and 'archive', since they were only separated for speed. ``` ## Looking up relationships between pairs of animals Now that the simulation is finished, examining kin relationships is the nearly same as the simple scenario. The only difference comes from the fact that we we integrated multiple years of sample-capture into the simulation, so we don't need to call `capture()` as a separate operation. ```{r} ## look up each animal's ancestors and look for shared ancestors between each ## pair of sampled animals: pairs <- findRelativesPar(indiv = indiv, sampled = TRUE, nCores = 2) POPs <- pairs[pairs$OneTwo == 1,] ## Parent-Offspring pairs GGPs <- pairs[pairs$OneThree == 1,] ## Grandparent-Grandoffspring pairs HSPs <- pairs[pairs$TwoTwo == 1,] ## Half-sibling pairs FSPs <- pairs[pairs$TwoTwo == 2,] ## Full Sibling pairs (self-comparisons ## are automatically excluded) FCPs <- pairs[pairs$ThreeThree == 2 & pairs$TwoTwo != 1,] ## Full Cousin pairs ## look at the number of shared ancestors at each ancestral ## generation, for one of the half-sibling pairs. lookAtPair(HSPs[1,]) relatives <- namedRelatives(pairs) ## shows the number of pairs of each relationship type relatives ``` \newpage # Addenda ## A note on mark-dependent survival or movement It is possible to simulate mark-dependent survival or movement by working on subsets of the 'indiv' matrix within the main simulation loop, and applying different mortality or movement operations to each subset. Note that if `mate()`, `altMate()`, `birthdays()`, etc., are run on a subset of the population, the operation will only affect members of that subset. In the case of subsets for marked and unmarked individuals, this would imply that marked animals only mate with other marked animals, and _vice versa_. It is important to ensure that a combined `indiv` object is available for any operations that affect the whole population, or to repeat each opertation for each subset. ```{r, eval = FALSE} indiv <- makeFounders(stocks = c(1)) ageMort <- c(0.305, rep(0.2, 100)) ## age-specific mortality is 0.305 for first-years, ## 0.2 for all older age-classes. for (y in 1:60) { indiv <- altMate(indiv, firstBreed = 2, batchSize = 0.9, year = y) ## 'mate' happens first, affecting the whole population indiv <- capture(indiv, n = 5, year = y, fatal = FALSE) indiv_unmarked <- indiv[is.na(indiv$SampY),] ## all the unmarked animals indiv_marked <- indiv[!is.na(indiv$SampY),] ## all the marked animals indiv_unmarked <- mort(indiv_unmarked, year = y, type = "age", ageMort = ageMort) ## normal mort for unmarked animals indiv_marked <- mort(indiv_marked, year = y, type = "age", ageMort = ageMort*2) ## double mort for marked animals indiv <- rbind(indiv_unmarked, indiv_marked) ## stick them back together and proceed as normal rm(indiv_marked, indiv_unmarked) ## clean up indiv <- birthdays(indiv) } ``` ## Arbitrary variables to affect the simulation Capture status isn't the only variable that can influence survival or movement. It is possible to add arbitrary additional variables to an `indiv` object in order to keep track of anything you might need to know about an individual, and it is possible to make simulation steps depend on those variables. You might want a 'hasBred' column, for instance, keeping track of which animals have produced offspring. Or, alongside a custom 'myMate' function (see 'custom functions', below), you might want to keep track of each animal's pair-bonded mate, so that animals can form pair bonds that persist across years. Here's an example of a species with two morphs, where the spiny morph has a slightly different movement-pattern to the smooth morph: ```{r arbitraryExtras} stocks <- c(0.3, 0.3, 0.4) indiv <- makeFounders( stocks = stocks) indiv[,10] <- sample(c("spiny", "smooth"), nrow(indiv), replace = TRUE) ## an extra column to hold stuff. Extras ## _must_ be in column numbers 10 or greater colnames(indiv)[10] <- "Morph" ## a movement matrix for spiny individuals admix.spiny <- matrix(NA, nrow = length(stocks), ncol = length(stocks)) for(i in 1:nrow(admix.spiny)) { admix.spiny[i,] <- stocks*stocks[i] } ## admix.spiny shows movement proportional to starting population sizes. admix.spiny ## a movement matrix for smooth individuals - let's make them more sedentary, ## with 90% chance of staying where they are and 5% chance of moving to each ## other stock admix.smooth <- matrix(c(0.05), nrow = length(stocks), ncol = length(stocks), byrow = FALSE) diag(admix.smooth) <- 0.9 admix.smooth <- admix.smooth / sum(admix.smooth) ## standardise so the matrix sums to 1 admix.smooth for (y in 1:20) { indiv <- altMate(indiv, firstBreed = 2, batchSize = 0.9, year = y) ## 'mate' happens first, affecting the whole population ## newborns need a 'smooth' or 'spiny' status, because altMate won't give one indiv$Morph[ is.na( indiv$Morph)] <- sample(c("spiny", "smooth"), nrow(indiv[ is.na( indiv$Morph),]), replace = TRUE) indiv <- mort(indiv, year = y, type = "flat", mortRate = 0.2) ## then 'mort', also affecting the whole population indiv_smooth <- indiv[indiv$Morph == "smooth",] ## all the smooth animals indiv_spiny <- indiv[indiv$Morph == "spiny",] ## all the spiny animals indiv_smooth <- move(indiv_smooth, moveMat = admix.smooth) ## one movement pattern for smooth animals indiv_spiny <- move(indiv_spiny, moveMat = admix.spiny) ## another movement pattern for spiny animals indiv <- rbind(indiv_smooth, indiv_spiny) ## stick them back together and proceed as normal rm(indiv_smooth, indiv_spiny) ## clean up indiv <- birthdays(indiv) } ``` ## Custom functions in fishSim The objects used in `fishSim` are deliberately kept simple, to encourage hacking wherever hacking is useful. To that end, most functions in `fishSim` both input and output a `data.frame` that we have called `indiv` throughout this vignette. If your species of interest has a life-history not readily captured by the existing `fishSim` functions, you may wish to write a custom function for just that step in the simulation. For instance, you may have a species in which each female mates with a handful of males in the year preceding each breeding season, and then uses stored sperm from these mates to fertilise a clutch of eggs just before giving birth. Let's call this the 'sperm storage' scenario. The sperm storage scenario can't be coded using the built-in `mate`, `altMate`, or `bondedMate` functions - in those functions, each mating event either has single paternity (i.e., all offspring produced by each female have the same father) or completely random paternity (i.e., the father of each offspring is randomly selected from among all mature males in the same stock), or, using `mate`, within-clutch single paternity but potentially multiple clutches (each with a different father). In the situations handled by `mate` and `altMate`, the expected number of offspring for a female is linearly related to her number of mates, which doesn't fit with our 'sperm storage' scenario. In the situation handled by `bondedMate`, one male and one female are socially-bonded for at least one breeding season, and a female's whole clutch are fathered by her mate. To get around the linear relationship between fecundity and number of mates, we will need to define a custom `mate`-like function, that inputs an `indiv`-like `data.frame`, generates new offspring the way the scenario demands, and outputs an `indiv`-like `data.frame` to be used by the other functions in `fishSim`. Such a function could look like this: ```{r} # A custom myMate function # # parameter 'indiv' is a standard indiv frame, as from makeFounders() # # parameter 'batchSize' sets the number of offspring per breeding female # such that the number of offspring ~Poisson(lambda = batchSize) # # parameter 'breedProb' is the probability that a mature female will # breed this in this mating event # # parameter 'maturityAge' is the knife-edge age at maturity, for both # males and females # # parameter 'year' is the year. It is used to set the $BirthY for new # recruits # # myMate uses a zero-truncated Poisson distribution to determine the # number of mates for each female that breeds. Parameter 'meanMates' # sets the T parameter for that distribution myMate <- function(indiv, batchSize, breedProb, maturityAge, year, meanMates) { library(ids) ## select mature females and males at current time-step matureFemales <- indiv[indiv$Sex == "F" & indiv$AgeLast >= maturityAge & is.na(indiv$DeathY) ,] matureMales <- indiv[indiv$Sex == "M" & indiv$AgeLast >= maturityAge & is.na(indiv$DeathY),] ## select subset of mature females that will breed (given breedProb) iBreed_F <- matureFemales[sample(1:nrow(matureFemales), size = nrow(matureFemales)*breedProb),] iBreed_M <- matureMales ## if no mature females or males available, don't add new individuals to the population if(nrow(iBreed_F) == 0 | nrow(iBreed_M) == 0) { return(indiv) } clutchSize <- rpois(nrow(iBreed_F), batchSize) # how many offsprings per breeding females? nMates <- rTruncPoisson(nrow(iBreed_F), T = meanMates) # how many mates per breeding females? mates <- lapply(nMates, function(i) sample(iBreed_M$Me, i, replace = FALSE)) # draw mates for each female based on nMates # and then assign fathers to each new recruit based on mates identified previously # (note each item in the list corresponds to a unique female) dads <- lapply(1:length(mates), function(i) sample(mates[[i]], clutchSize[[i]], replace = TRUE)) # now set-up table of recruits assembling all variables recruits <- data.frame(Me=uuid(sum(clutchSize), drop_hyphens = TRUE), # unique ID Sex=sample(c("M","F"), sum(clutchSize), prob = c(0.5,0.5), replace = TRUE), Mum=rep(iBreed_F$Me, clutchSize), # add mum vector that match litter sizes assigned earlier Dad=unlist(dads), # unlist dads (length should match new recruits number) BirthY=year, # rest of meta-data is simple for recruits. birth year DeathY=NA, # not dead yet Stock=1, # single stock for now AgeLast=0, # recruits start at age 0 SampY=NA) # not sampled yet outs <- rbind(indiv, recruits) return(outs) } ``` With `myMate` defined, we can go on to write a simulation using the other `fishSim` tools, as in this stripped-down example: ```{r} indiv <- makeFounders(stocks = c(1)) ageMort <- c(0.305, rep(0.2, 100)) ## age-specific mortality as in the ## first, 'simple' example for (y in 1:10) { indiv <- myMate(indiv = indiv, batchSize = 3, breedProb = 1, maturityAge = 4, year = y, meanMates = 3) indiv <- mort(indiv, year = y, type = "age", ageMort = ageMort) indiv <- birthdays(indiv) if(y %in% c(6:10)) { indiv <- capture(indiv, n = 40, fatal = FALSE, year = y) } } pairs <- findRelativesPar(indiv = indiv, sampled = TRUE, nCores = 2) namedRelatives(pairs) ``` ### Pair-bonded mating