# Introduction

This vignette consists of three brief investigations of mating scenes using package mateable. A mating scene is a bout of mating where the coordinates of participating individuals are defined in one, two, or three dimensions: space, time, and compatibility. From such information we can quantify mating potential, the capacity for sexual reproduction based on the location, reproductive timing, and compatibility of prospective mates. Mating potential can be quantified for a pair of individuals based on the distance between them, the timing of their reproductive activity, and their compatibility. Similarly, mating potential can be defined for an individual within the context of the mating scene or for an entire scene.

library(mateable)
packageDescription("mateable")$Version ## [1] "0.3.1" For any function in mateable you can learn more by typing a question mark before the function name, i.e. ?makeScene. To learn all the functions, type ?mateable and click index at the bottom of the page. # A real dataset In this section we look at flowering during 2012 in a remnant prairie population of Echinacea angustifolia, the narrow-leaved purple coneflower. The data frame eelr2012 is included in the package. str(eelr2012) ## 'data.frame': 44 obs. of 6 variables: ##$ tagNo   : int  1509 7614 8647 9021 9022 9520 9582 10658 12243 13042 ...
##  $heads : num 6 2 1 1 1 2 1 1 3 1 ... ##$ firstDay: Date, format: "2012-06-25" "2012-06-23" ...
##  $lastDay : Date, format: "2012-07-08" "2012-07-03" ... ##$ Ecoord  : num  32.3 174.6 132.3 170.8 160.4 ...
##  $Ncoord : num 36 17.5 68.1 18 19.5 ... This data set includes spatial and temporal information on all 44 plants that flowered in 2012 from the East Elk Lake Road population. Each plant has a tag with a unique number for identification purposes and we listed the count of heads that the plant produced in 2012. These long-lived plants do not produce heads in every year, we removed all plants with zero heads in 2012. The columns firstDay and lastDay indicate the first and last days that each plant produced pollen. Spatial coordinates are in meters. The first thing to do is convert the data frame to a mating scene object. We must indicate which columns contain spatial and temporal coordinates. Note that some columns are integer and numeric. Columns firstDay and lastDay are saved in a Date format. makeScene can read some character formats. If you want to convert character or POSIX to the Date format, read about function as.Date or install package lubridate. eelr <- makeScene(eelr2012, startCol = "firstDay", endCol = "lastDay", xCol = "Ecoord", yCol = "Ncoord", idCol = "tagNo") We can do a lot with this matingScene object. Let’s start by focusing on the spatial dimension of the mating scene. ## The spatial dimension of the mating scene We can look at the spatial “s” dimension with a map. Note that we can use standard R graphical parameters to change the looks of the graph. plotScene(eelr, "s") plotScene(eelr, "s", pch = 1) We can also get the distances between all pairs using function pairDist, and then visualize their distribution. ePair <- pairDist(eelr) # matrix of distances hist(ePair, 40) # visualize histogram We can also calculate the distance from each individual to its nearest neighbors (aka potential mates). Here we visualize the distance of every plant to its 3rd nearest neighbor. eKnn <- kNearNeighbors(eelr, 6) # 1,2,...,6 nearest potential mate str(eKnn) ## num [1:44, 1:6] 29.682 1.738 5.552 0.967 3.23 ... ## - attr(*, "dimnames")=List of 2 ## ..$ : chr [1:44] "1509" "7614" "8647" "9021" ...
##   ..$: chr [1:6] "k1" "k2" "k3" "k4" ... hist(eKnn[, 3], 40, main = "Distance to 3rd nearest neighbor (m)") Function matingSummary calculates many characteristics of the “whole scene.” We can save summaries as a new object eSum and look at characteristics of the entire mating scene, either by indexing or by name. Here are some spatial summaries. eSum <- matingSummary(eelr) eSum[c("minX", "minY", "k1")] # by index ##$minX
## [1] 24.37
##
## $minY ## [1] 11.24 ## ##$k1
## [1] 29.68228
# by name
eSum$minX ## [1] 24.37 eSum$minY
## [1] 11.24
eSum$k1 # mean distance to 1st nearest potential mate ## [1] 29.68228 Distance is a measure of isolation from mates. To characterize mating potential, which is inversely related to distance, we want proximity to mates. For this, we have function proximity. eProx <- proximity(eelr, "maxPropSqrd") eProx$pop
## [1] 0.4970269
hist(eProx$ind$proximity, 30)

Now that we have saved a proximity object, eProx, we can visualize it using function plotPotential. There are many ways to visualize it. The defaults return a three panel figure: pairwise values of proximity, including a histogram of all pairs, a network diagram of a random subset of 9 individuals, and a heat map of all interactions between those same nine. In the network diagram, line width corresponds to the pairwise proximity of the individuals being connected, and the size of label indicates the individual’s mean proximity with all other individuals in the network.

plotPotential(eProx)

If we want to focus on a particular subset of plants, then we can define them. Here, we make only a network diagram of these focal plants. Note that by default the size of the individual labels on the network diagram scales with the individual proximity, but this can be turned off.

focalPlants <- c(17217, 17202, 14582, 15114, 7614, 1509, 17002)
plotPotential(eProx, plotType = "net", sub.ids = focalPlants)

Notice plant 1509 is more isolated from potential mates than the other individuals.

## The temporal dimension of the mating scene

Now we turn from the spatial dimension of the mating scene to the temporal dimension. We can visualize the coordinates of mating activity in time with a mating schedule. The function plotScene works just as well in time as in space, we just need to specify which dimension–use “t” for time.

plotScene(eelr, "t")

plotScene(eelr, "t", sub = 'all', drawQuartiles = FALSE, text.cex = 0.5)

These figures shows the duration of reproductive period during 2012 for each of the 44 individuals in the population with a horizontal line starting on the date the individual first produced pollen and ending on the date when pollen was last produced. Here the individuals are sorted from bottom to top by their first day. The dots indicate the total number of individuals participating in mating on each day.

Just as we can calculate and visualize the spatial distances between all pairs using function pairDist, we can do the same in the temporal dimension with function overlap. Function overlap makes a matrix of overlapping days for all pairs and we can use a histogram to see the distribution.

eOver <- overlap(eelr, compareToSelf = T) # matrix of days overlapping
hist(eOver, 40, main = "Histogram of days overlapping between pair") 

We can look at the number of individuals participating in mating every day.

eRecep <- receptivityByDay(eelr) # T/F receptive on each day
str(eRecep) # matrix
##  logi [1:44, 1:36] FALSE FALSE FALSE FALSE FALSE FALSE ...
##  - attr(*, "dimnames")=List of 2
##   ..$: chr [1:44] "1509" "7614" "8647" "9021" ... ## ..$ : chr [1:36] "1" "2" "3" "4" ...
##  - attr(*, "origin")= Date[1:1], format: "2012-06-13"
dailyReceptivitySummary <- receptivityByDay(eelr, summary = TRUE)
dailyReceptivitySummary # a named integer vector
## 2012-06-13 2012-06-14 2012-06-15 2012-06-16 2012-06-17 2012-06-18
##          1          1          1          1          1          1
## 2012-06-19 2012-06-20 2012-06-21 2012-06-22 2012-06-23 2012-06-24
##          3          3          4          7          9         14
## 2012-06-25 2012-06-26 2012-06-27 2012-06-28 2012-06-29 2012-06-30
##         17         19         25         32         35         39
## 2012-07-01 2012-07-02 2012-07-03 2012-07-04 2012-07-05 2012-07-06
##         38         38         38         35         33         29
## 2012-07-07 2012-07-08 2012-07-09 2012-07-10 2012-07-11 2012-07-12
##         27         19         15         14         10          8
## 2012-07-13 2012-07-14 2012-07-15 2012-07-16 2012-07-17 2012-07-18
##          8          5          3          2          2          2
plot(as.Date(names(dailyReceptivitySummary)), dailyReceptivitySummary,
xlab = 'date', ylab = 'count of receptive individuals')

Folks have calculated synchrony many ways for individuals, pairs, and populations. Our function synchrony does it all; well, it does a lot. We can specify many different synchrony measures, with enticing names, such as these (read the help for details): “augspurger”, “kempenaers”, “sync_prop”, “overlap”, “sync_nn”, “simple1”, “simple2”, and “simple3.” Also, we can calculate measures for different subjects: individual, pairs, and the whole scene (population).

Here, for example, we calculate synchrony for all subjects using the overlap method and save it as eSync. Then we show the mean and median population values in red and blue, respectively.

eSync <- synchrony(eelr, "overlap")
hist(eSync$ind[, 2], 30) abline(v = eSync$pop, col ="red", lwd = 2)
abline(v = synchrony(eelr, "overlap", averageType = "median")$pop, col = "blue", lwd = 2) Individual overlap indicates the proportion of potentital mates in the scene that were flowering averaged over all days that the focal individual participated in mating. Just like we visualized spatial mating potential, we can visualize mating potential in the temporal dimension. Here we emphasize the same focal plants as we did above. plotPotential(eSync, sub.ids = focalPlants) Notice that 1509 is not isolated in time from potential mates. ## Time and space together Function plotScene is useful because it enables us to visualize all dimensions of the mating scene (in this data set we only have two dimensions). In the below figure we also highlight the focal plants we selected above so we can see where and when they participate in mating. plotScene(eelr, c('s','t'), sub = focalPlants, N = 4, text.cex = 0.5) Function plot3DScene is also useful because it enables us to visualize two dimensions of the mating scene in one panel. plot3DScene(eelr) plot3DScene(eelr, pt.cex = .5, sub = focalPlants, N = 4) Function matingSummary calculates many characteristics of the entire mating scene in all three dimensions. Though in this example data set we have only dimensions. eSum <- matingSummary(eelr) eSum[c("meanSD", "sdSD", "meanDur", "sdDur", "peak")] # index eSum[c("minX", "minY", "k1")] eSum$minX # by name
eSum$minY eSum$k1
## $meanSD ## [1] "2012-06-27" ## ##$sdSD
## [1] 4.75488
##
## $meanDur ## [1] 12.25 ## ##$sdDur
## [1] 3.133651
##
## $peak ## [1] "2012-07-01" ## ##$minX
## [1] 24.37
##
## $minY ## [1] 11.24 ## ##$k1
## [1] 29.68228
##
## [1] 24.37
## [1] 11.24
## [1] 29.68228

# Compatibility

In animals, females are compatible mates with males but they are not compatible with other females. Similarly, males are only compatible with females. We say that mating potential between a pair of males is zero, between a pair of females is zero, and between a mixed-sex pair it is 1. Most animal species, though not all, have a breeding system like this. It’s more complicated in plants. Plants have many breeding systems, including self-compatibility and dioecy–dioecy is just like the simple case in animals. About half of all plants species have some kind of self-incompatibility system which makes it possible for some pairs to be mating incompatible, i.e. have mating potential of zero. Package mateable allows us to model the animal breeding system and the breeding system in Echinacea (sporophytic self-incompatibility). In the following examples, we show examples of each. We intend to add capability for handling more breeding systems in future releases of mateable.

# A simulated dataset

Along with using data from real populations, we can also simulate scenes. Here we simulate a scene using values from the eelr summary as inputs for the simulation parameters.

# make scene based off eelr summary information
simScene <- simulateScene(size = nrow(eelr), meanSD = eSum$meanSD, sdSD = eSum$sdSD, meanDur = eSum$meanDur, sdDur = eSum$sdDur, xRange = c(eSum$minX, eSum$maxX),
yRange = c(eSum$minY, eSum$maxY))

A simulated scene can be treated the same way as a scene made from real data, meaning all functions work in the same way. In addition, simulated scenes will always have all three dimensions (space, time, and compatibility), so you can examine any aspect of mating potential.

sProx <- proximity(simScene, "maxPropSqrd")
sSync <- synchrony(simScene, "augspurger")
sComp <- compatibility(simScene, "si_echinacea")

The coordinates in the spatial and temporal dimensions are generated from uniform and normal distributions, respectively. The coordinates in the compatibility dimension are, for each individual, two alleles selected at random from the number of total alleles in the scene (default is 10). In the third panel of the second figure allele shows the alleles for each individual; they are labeled 1, 2, 3, … , 10.

plotScene(simScene)
plot3DScene(simScene, sub = 'random')

plotPotential(sSync)
plotPotential(sProx)

Mating potentials in space and time are continuous because space and time are continuous. In contrast, mating compatibility between two individuals, as we have modeled it here, is either possible or not. Therefore mating potential for a pair is either 1 or 0.

plotPotential(sComp, density = FALSE)

plotPotential(sComp, plotType = c('net','heat'), density = FALSE)

Plot3DPotential gives us an idea of relationships between different dimensions of potential.

plot3DPotential(list(sProx, sSync), subject = "ind")

plot3DPotential(list(sProx, sSync), subject = "ind", sample = 'all')

plot3DPotential(list(sProx, sSync), subject = "ind", sample = 'random')

plot3DPotential(list(sProx, sSync, sComp))

plot3DPotential(list(sProx, sSync, sComp), sub.ids = c(21,4))

Because these scenes are simulated, we do not expect to see correlations in potential between dimensions.

## Compatibility in animals

For those who work with animals, have no fear, you can create mating scene objects too! When using simulated data, use the option sAlleles = 2 and choose either the number 1 or 2 to represent female or male (it doesn’t matter which, they’re randomly generated. M is 2 by default). If using real data, create a column with 1 or 2 for female or male and specify that column as s1Col in makeScene. When using compatibility, choose the method “dioecious” and compatibility will be 1 if they are not the same sex (i.e. not both 1 in column s1).


simulatedCoral <- simulateScene(10, sAlleles = 2)
plotScene(simulatedCoral)

plot3DScene(simulatedCoral)

plotPotential(compatibility(simulatedCoral, "dioecious"))

We can look at the matingScene object called “simulatedCoral”. Column s2 is ignored.

simulatedCoral
##    id start end duration        x         y s1 s2
## 1   1     3  14       12 61.35159 57.347986  2  1
## 2   2    12  26       15 73.33633 64.840381  1  2
## 3   3     8  19       12 85.51832 75.432383  1  2
## 4   4     1  11       11 73.64774 38.774773  1  2
## 5   5     9  23       15  2.18483 71.542488  1  2
## 6   6     5  20       16 14.25105  6.018506  1  2
## 7   7     1  12       12 92.90814 74.096211  2  1
## 8   8     5  18       14 10.05030 35.491184  1  2
## 9   9    11  17        7 39.18539 76.007530  2  1
## 10 10    22  29        8 51.62646 88.184330  2  1
str(simulatedCoral)
## 'data.frame':    10 obs. of  8 variables:
##  $id : int 1 2 3 4 5 6 7 8 9 10 ##$ start   : num  3 12 8 1 9 5 1 5 11 22
##  $end : num 14 26 19 11 23 20 12 18 17 29 ##$ duration: num  12 15 12 11 15 16 12 14 7 8
##  $x : num 61.35 73.34 85.52 73.65 2.18 ... ##$ y       : num  57.3 64.8 75.4 38.8 71.5 ...
##  $s1 : Factor w/ 2 levels "1","2": 2 1 1 1 1 1 2 1 2 2 ##$ s2      : Factor w/ 2 levels "1","2": 1 2 2 2 2 2 1 2 1 1
##  - attr(*, "t")= logi TRUE
##  - attr(*, "s")= logi TRUE
##  - attr(*, "mt")= logi TRUE
##  - attr(*, "originalNames")= chr  "id" "start" "end" "x" ...
##  - attr(*, "origin")= Date, format: "2012-07-04"
summary(simulatedCoral)
##        id            start           end           duration
##  Min.   : 1.00   Min.   : 1.0   Min.   :11.00   Min.   : 7.00
##  1st Qu.: 3.25   1st Qu.: 3.5   1st Qu.:14.75   1st Qu.:11.25
##  Median : 5.50   Median : 6.5   Median :18.50   Median :12.00
##  Mean   : 5.50   Mean   : 7.7   Mean   :18.90   Mean   :12.20
##  3rd Qu.: 7.75   3rd Qu.:10.5   3rd Qu.:22.25   3rd Qu.:14.75
##  Max.   :10.00   Max.   :22.0   Max.   :29.00   Max.   :16.00
##        x                y          s1    s2
##  Min.   : 2.185   Min.   : 6.019   1:6   1:4
##  1st Qu.:20.485   1st Qu.:43.418   2:4   2:6
##  Median :56.489   Median :68.191
##  Mean   :50.406   Mean   :58.774
##  3rd Qu.:73.570   3rd Qu.:75.098
##  Max.   :92.908   Max.   :88.184

# Multi-year scenes

Here we investigate the dynamics of a population’s mating scene over multiple seasons. Multi-year scenes must be formatted as lists, with each list element representing one mating scene (or one year). If we had a multi-year dataset, function makeScene would make it into a multi-year scene automatically using the argument multiYear = TRUE. Alternatively, we can simulate several mating scenes and combine them in a list by hand, as in the example below.

simScene1 <- simulateScene(size = nrow(eelr), meanSD = eSum$meanSD, sdSD = eSum$sdSD, meanDur = eSum$meanDur, sdDur = eSum$sdDur, xRange = c(eSum$minX, eSum$maxX),
yRange = c(eSum$minY, eSum$maxY))
simScene2<- simulateScene(size = 1.5*nrow(eelr), meanSD = eSum$meanSD + 365, sdSD = eSum$sdSD, meanDur = eSum$meanDur, sdDur = eSum$sdDur, xRange = c(eSum$minX, eSum$maxX),
yRange = c(eSum$minY, eSum$maxY))
simScene3 <- simulateScene(size = 0.8*nrow(eelr), meanSD = eSum$meanSD + 730, sdSD = eSum$sdSD, meanDur = eSum$meanDur, sdDur = eSum$sdDur, xRange = c(eSum$minX, eSum$maxX),
yRange = c(eSum$minY, eSum$maxY))

multiYearScene <- list('2012' = simScene1,'2013' = simScene2, '2014' = simScene3)

All of the different analysis and plotting methods can be applied directly to multi-year scenes. We can make multi-panel plots of the matingScene over years. The plot limits are consistent across years, making it easier to compare differences.

plotScene(multiYearScene,sub = c(1,6,12,18,13,24,55,45,60), text.cex = 0.8)

We can also combine mating dimensions for multi-year plots using plot3DScene.

plot3DScene(multiYearScene, pt.cex = 1.2, sub = c(1,6,12,18,13,24,55,45,60))

The functions synchrony, proximity, and compatibility also work on multi-year scenes, returning a list of potentials objects for each year.

syncMulti <- synchrony(multiYearScene, method = 'augs')
proxMulti <- proximity(multiYearScene, method = 'maxPropSqrd')
compatMulti <- compatibility(multiYearScene, method = 'si_echinacea')

str(syncMulti) # a list of lists
## List of 3
##  $2012:List of 3 ## ..$ pop : num 0.685
##   ..$ind :'data.frame': 44 obs. of 2 variables: ## .. ..$ id       : int [1:44] 1 2 3 4 5 6 7 8 9 10 ...
##   .. ..$synchrony: num [1:44] 0.7818 0.8118 0.0426 0.6237 0.7359 ... ## ..$ pair: num [1:44, 1:44] 1 1 0 0.636 0.857 ...
##   .. ..- attr(*, "idOrder")= int [1:44] 1 2 3 4 5 6 7 8 9 10 ...
##   ..- attr(*, "t")= logi TRUE
##   ..- attr(*, "s")= logi FALSE
##   ..- attr(*, "c")= logi FALSE
##  $2013:List of 3 ## ..$ pop : num 0.553
##   ..$ind :'data.frame': 66 obs. of 2 variables: ## .. ..$ id       : int [1:66] 1 2 3 4 5 6 7 8 9 10 ...
##   .. ..$synchrony: num [1:66] 0.599 0.69 0.327 0.612 0.57 ... ## ..$ pair: num [1:66, 1:66] 1 1 0.5 1 0.846 ...
##   .. ..- attr(*, "idOrder")= int [1:66] 1 2 3 4 5 6 7 8 9 10 ...
##   ..- attr(*, "t")= logi TRUE
##   ..- attr(*, "s")= logi FALSE
##   ..- attr(*, "c")= logi FALSE
##  $2014:List of 3 ## ..$ pop : num 0.517
##   ..$ind :'data.frame': 35 obs. of 2 variables: ## .. ..$ id       : int [1:35] 1 2 3 4 5 6 7 8 9 10 ...
##   .. ..$synchrony: num [1:35] 0.394 0.606 0.603 0.473 0.658 ... ## ..$ pair: num [1:35, 1:35] 1 0.154 0.6 0 0.364 ...
##   .. ..- attr(*, "idOrder")= int [1:35] 1 2 3 4 5 6 7 8 9 10 ...
##   ..- attr(*, "t")= logi TRUE
##   ..- attr(*, "s")= logi FALSE
##   ..- attr(*, "c")= logi FALSE

The functions plotPotential and plot3DPotential will then work on multi-year potential objects. For example, we can visualize synchrony over multiple years. Note that these functions will try to select the same sample of individuals across years, so if there is a year when few individuals in that sample are participating in the mating scene, there will be fewer individuals displayed in that year’s heatmap and network diagram.

plotPotential(syncMulti, sub.ids = c(1,6,12,18,13,24,55,44,60))

And we can visualize proximity over multiple years.

plotPotential(proxMulti,sub.ids = c(1,6,12,18,13,24,55,45,60))

And, like before, we can visualize combinations of synchrony, proximity, and compatibility over multiple years.

plot3DPotential(list(syncMulti, proxMulti, compatMulti), subject = 'ind',
pt.cex = 1, sub.ids = c(1,6,12,18,13,24,55,45,60))

## Conclusion

We find mateable useful and hope you do too. It can get better, so we are improving it. A development version of mateable is available via gitHub. We welcome user suggestions for improvements. Please submit bugs and feature requests to the mateable development page on github or contact Stuart directly.