Extra-pair paternity in Blue Tits (_Cyanistes caeruleus_): a case study from Westerholz, Bavaria, Germany

In expp: Spatial Analysis of Extra-Pair Paternity

---

title: "Extra-pair paternity in Blue Tits (Cyanistes caeruleus): a case study from Westerholz, Bavaria, Germany" author: "Mihai Valcu" date: "r Sys.Date()" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Blue Tits Case study} %\VignetteEngine{knitr::rmarkdown}

---

Supplement to Schlicht, Valcu and Kempenaers Schlicht, Lotte, Mihai Valcu, and Bart Kempenaers. "Spatial patterns of extra-pair paternity: beyond paternity gains and losses." Journal of Animal Ecology 84.2 (2015): 518-531.

---

1. Getting started

• Download and install R.
• Open R, and install expp with install.packages("expp") :

require(rgeos); require(sp); require(spdep); require(deldir)
par(mar = c(0,0,1,0) )

• To load expp type r require(expp):

2. Load datasets

For info on the data-sets type:

help(bluetit_breeding)
help(bluetit_epp)
help(bluetit_boundary)

data(bluetit_breeding)
head(bluetit_breeding[bluetit_breeding$year_ == 2011, ])

data(bluetit_epp)
head(bluetit_epp[bluetit_epp$year_ == 2011, ])

data(bluetit_epp)
knitr::kable(head(bluetit_epp[bluetit_epp == 2011, ]))

data(bluetit_boundary)
summary(bluetit_boundary)

3. Prepare data

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

b = split(bluetit_breeding, bluetit_breeding$year_)
e = split(bluetit_epp, bluetit_epp$year_) 

# sample sizes by year

# number of breeding pairs
sapply(b, nrow)

# number of extra-pair events
sapply(e, nrow)

# For the sake of conciseness only two years are used in the folowing analyses
b = b[c("2009", "2010")]
e = e[c("2009", "2010")]
p = bluetit_boundary[bluetit_boundary$year_ %in% c("2009", "2010"), ]

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, 
  proj4string = CRS(proj4string(p)))

eppDat = lapply(e, eppMatrix, pairs = ~ male + female)

3.3. Compute Dirichlet polygons based on the SpatialPointsBreeding object

polygonsDat = mapply(DirichletPolygons, x = breedingDat, boundary = split(p, p$year_)) 

---

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

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

op = par(mar = c(0,0,2,0))

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

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

year = '2010'
box = 110
eppOut10 = eppOut[[year]]
plot(eppOut10 , zoom = box, maxlag = 2,cex = .7,  border = 'white', col = 'grey70', zoom.col = "bisque")

par(op)

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

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

par(op)

5. Fitting a glmm

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

dat = lapply(eppOut, as.data.frame) # a list of data.frame(s)
dat = do.call(rbind, dat)
dat$year_ = dat$year__MALE; dat$year__FEMALE = NULL

5.2. Data transformations prior to modelling.

Rescale rank; rank 1 becames rank 0

dat$rank = dat$rank - min(dat$rank)
table(dat$rank)

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.

knitr::kable(table(dat$epp, dat$year_))

Run the glmm model (this may take a while depending on your system!).

require(lme4)
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 (Laplace Approximation) ['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 df.resid 
## 599.4    658.8   -291.7    583.4    12406 

## Scaled residuals: 
## Min     1Q Median     3Q    Max 
## -0.568 -0.048 -0.017 -0.006 96.102 

## Random effects:
## Groups Name        Variance  Std.Dev. 
## male   (Intercept) 1.245e+00 1.116e+00
## female (Intercept) 8.376e-02 2.894e-01
## year_  (Intercept) 3.121e-15 5.586e-08
## Number of obs: 12414, groups:  male, 121; female, 118; year_, 2

## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                -3.325738   0.331005 -10.047  < 2e-16 ***
## rank                       -1.166547   0.132700  -8.791  < 2e-16 ***
## male_age_MALEjuv           -1.380823   0.418108  -3.303 0.000958 ***
## relative_asynchrony_MALE   -0.476106   0.402514  -1.183 0.236876    
## relative_asynchrony_FEMALE -0.004861   0.376569  -0.013 0.989700    
## ---

## Correlation of Fixed Effects:
## (Intr) rank   m__MAL r__MAL
## rank        -0.272                     
## ml_g_MALEjv -0.198  0.025              
## rltv_s_MALE  0.075  0.022  0.004       
## rlt__FEMALE  0.006 -0.004 -0.003 -0.393

expp documentation built on June 20, 2021, 5:06 p.m.