Nothing
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}
install.packages("expp")
: require(rgeos); require(sp); require(spdep); require(deldir) par(mar = c(0,0,1,0) )
r require(expp)
: 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)
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"), ]
breedingDat = lapply(b, SpatialPointsBreeding, coords= ~x+y, id='id', breeding= ~male + female, proj4string = CRS(proj4string(p))) eppDat = lapply(e, eppMatrix, pairs = ~ male + female)
SpatialPointsBreeding
objectpolygonsDat = mapply(DirichletPolygons, x = breedingDat, boundary = split(p, p$year_))
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) }
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)
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
dat$rank = dat$rank - min(dat$rank) table(dat$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)
table(dat$epp, dat$year_) #extra-pair frequency by year.
knitr::kable(table(dat$epp, dat$year_))
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
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.