# Extra-pair paternity in Blue Tits (Cyanistes caeruleus): a case study from Westerholz, Bavaria

## Supplement to Schlicht, Valcu and Kempenaers “Spatial patterns of extra-pair paternity: beyond paternity gains and losses” (in prep.)

### 1. Getting started

• Open R, and install expp by copying the following line into your R console:
install.packages("expp")

require(expp)


For info on the data-sets type:

help(westerholzBreeding)
help(westerholzEPP)

data(westerholzBreeding)

data(westerholzEPP)

year_ id x y layingDate male_age female_age female male
2010 214 4417785 5334719 111 adult juv f63 m35

year_ female male
2009 f29 m16
2009 f53 m16
2009 f36 m16
2009 f47 m41
2009 f51 m79
2009 f38 m13

### 3. Prepare data

#### 3.1 Split by year (each year needs to be processed separately)

b = split(westerholzBreeding, westerholzBreeding$year_) e = split(westerholzEPP, westerholzEPP$year_)

# sample sizes by year
lapply(b, nrow)

## $2009 ## [1] 56 ## ##$2010
## [1] 96

lapply(e, nrow)

## $2009 ## [1] 29 ## ##$2010
## [1] 32


#### 3.2 Run a couple of helper functions on both breeding data and extra-pair paternity data

breedingDat = lapply(b, SpatialPointsBreeding, coords= ~x+y, id='id', breeding= ~male + female)
eppDat = lapply(e, eppMatrix, pairs = ~ male + female)



#### 3.3. Compute Dirichlet polygons based on the SpatialPointsBreeding object

polygonsDat = lapply(breedingDat, DirichletPolygons)


### 4. All the objects are now ready to be processed by the epp function.

maxlag = 10
O = mapply(FUN = epp, breedingDat, polygonsDat, eppDat, maxlag)

## Warning: maxlag of 10 is exceeding the number of possible lags, reverting
## to 8 lags.

# op = par(mfrow = c(1,2))

for (year in c("2009", "2010")) {
plot(O[[year]], cex = 0.7, lwd = 0.5, border = "navy")
title(main = year)
}

## Warning: EP male 'm41' not found.



# par(op)


#### Select one nest-box of a given year and zoom in.

year = "2010"
box = 110
O10 = O[[year]]
plot(O10, zoom = box, maxlag = 2, cex = 0.7, border = "white", col = "grey70",
zoom.col = "bisque")


op = par(mfrow = c(1, 2))

barplot(O[[1]], relativeValues = TRUE, main = 2009)
legend(x = "topright", legend = c("Observed", "Potential"), lty = c(1, 2), bty = "n")
barplot(O[[2]], relativeValues = TRUE, main = 2010)



par(op)


### 5. Fitting a glmm

#### 5.1 Convert O (a list of 2 epp objects) into a data.frame.

dat = lapply(O, as.data.frame)  # a list of data.frame-s
dat = do.call(rbind, dat)
dat$year_ = dat$year__MALE
table(dat$rank)  ## ## 0 1 2 3 4 5 6 7 8 9 ## 776 1404 1862 2012 1890 1542 1180 774 434 230  ##### Center and re-scale breeding asynchrony (i.e. the difference in laying data between male and female) within each rank. center = function(x) { return(x - mean(x, na.rm = TRUE)) } scale2 = function(x) { return(x/(2 * sd(x, na.rm = TRUE))) } # Compute asynchrony dat$asynchrony = abs(dat$layingDate_MALE - dat$layingDate_FEMALE)

# a Compute relative within-rank asynchrony
MALE_splitBy = paste(dat$year_, dat$id_MALE, dat$male, dat$rank, sep = "_")
dat$relative_asynchrony_MALE = unsplit(lapply(split(dat$asynchrony, MALE_splitBy),
center), MALE_splitBy)
dat$relative_asynchrony_MALE = scale2(dat$relative_asynchrony_MALE)

FEMALE_splitBy = paste(dat$year_, dat$id_FEMALE, dat$female, dat$rank, sep = "_")
dat$relative_asynchrony_FEMALE = unsplit(lapply(split(dat$asynchrony, FEMALE_splitBy),
center), FEMALE_splitBy)
dat$relative_asynchrony_FEMALE = scale2(dat$relative_asynchrony_FEMALE)


#### 5.3 Run glmm

##### Check if sample size is sufficient for the number of variables we aim to include into the model.
table(dat$epp, dat$year_)  #extra-pair frequency by year.

id 2009 2010
0 3053 8991
1 27 33
##### Run the glmm model (this may take a while depending on your system!).
require(lme4)
dat$age2 = ifelse(dat$male_age_MALE == "juv", 1, 2)
fm = glmer(epp ~ rank + male_age_MALE + relative_asynchrony_MALE + relative_asynchrony_FEMALE +
(1 | male) + (1 | female) + (1 | year_), data = dat, family = binomial)
summary(fm)

## Generalized linear mixed model fit by maximum likelihood ['glmerMod']
##  Family: binomial ( logit )
## Formula: epp ~ rank + male_age_MALE + relative_asynchrony_MALE + relative_asynchrony_FEMALE +      (1 | male) + (1 | female) + (1 | year_)
##    Data: dat
##
##      AIC      BIC   logLik deviance
##    588.3    647.5   -286.2    572.3
##
## Random effects:
##  Groups Name        Variance Std.Dev.
##  male   (Intercept) 1.12e+00 1.05790
##  female (Intercept) 9.87e-02 0.31423
##  year_  (Intercept) 1.54e-05 0.00393
## Number of obs: 12104, groups: male, 119; female, 117; year_, 2
##
## Fixed effects:
##                            Estimate Std. Error z value Pr(>|z|)
## (Intercept)                  -3.275      0.245  -13.39  < 2e-16 ***
## rank                         -1.237      0.156   -7.95  1.9e-15 ***
## male_age_MALEjuv             -1.365      0.449   -3.04   0.0023 **
## relative_asynchrony_MALE     -0.386      0.408   -0.95   0.3442
## relative_asynchrony_FEMALE    0.076      0.376    0.20   0.8397
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
##             (Intr) rank   m__MAL r__MAL
## rank        -0.449
## ml_g_MALEjv -0.382  0.010
## rltv_s_MALE  0.091  0.017  0.000
## rlt__FEMALE -0.008 -0.001 -0.024 -0.345