## Supplementary analyses for the project on multiple nest box great tit populations---- ## to investigate viability fluctuating (density dependent) selection on exploration behaviour ## Analyse whether locally and immigrant recruited individuals have similar exploration behavior ## to explore the possibility that recruitment bias might be larger for fast explorers. ## Populations: # Starnberg (Bavaria, Germany): Data from Niels Dingemanse and Alexia Mouchet LMU Munich; years 2010 to 2014, 12 plots # Lauwersmeer (Groningen, the Netherlands): Data from Marion Nicolaus and Joost Tinbergen; years 2006 to 2009, 12 plots # Westerheide (Wageningen NIOO, the Netherlands): Data from Kees van Oers; years 1999 to 2016, 1 plot # Boshoek, Antwerp (Belgium): Data from Erik Matthysen and Frank Adriaensen; years 2006 to 2017, 9 plots; Co-authors: Erik Matthysen and Mathijs van Overveld # Wytham (Oxford, UK): Data from John Quinn; years 2005 to 2016, 1 plot- Other co-authors: Ela Colle and Allison Roth ## Author/Date: Alexia Mouchet (LMU Munich), April 2021 --------------------------------------------------------------------------------------------------------------------------------------------------------------- # Libraries---- library (tidyverse) library(plyr) library(dplyr) library(lme4) library(INLA) ## Starnberg population---- dat_sta = read.table("Dataset_S1_Starnberg.csv", sep = ",", header = T) dat_sta_local = read.table("Dataset_S1.2_LocalRecruits_Starnberg.csv", sep = ",", header = T) # Merge behavioral data with immigration status data dat_sta_rec = merge (dat_sta, dat_sta_local[, c("BirdID", "Immigrant", "BirthYear")], by = "BirdID", all.x=TRUE) # Factorize age field and select only yearlings (age= 5) dat_sta_rec$Age = factor(dat_sta_rec$Age) dat_sta_age5 = dat_sta_rec[dat_sta_rec$Age=="5",] # Add value for immigrants (local dataset contained data only for local recruits) dat_sta_age5 = dat_sta_age5 %>% replace_na(list(Immigrant = "1")) # Add Population ID to dataset dat_sta_age5$PopulationID = rep(1, nrow(dat_sta_age5)) # Calculate exploration score corrected for time of year (July 1st as reference) # Exploration July date dat_sta_age5$expdate = as.Date(dat_sta_age5$ExpDate,"%d/%m/%Y") dat_sta_age5$JulyYear = ifelse (as.numeric(format(dat_sta_age5$expdate, "%m"))<7, as.numeric(format(dat_sta_age5$expdate, "%Y"))-1, as.numeric(format(dat_sta_age5$expdate, "%Y"))) dat_sta_age5$July1st = as.Date(paste(dat_sta_age5$JulyYear, "07", "01", sep = "-"), origin = "1899-12-30","%Y-%m-%d") dat_sta_age5$ExpJulyDate = dat_sta_age5$expdate - dat_sta_age5$July1st # Correct exploration score for time of year (July 1st as reference) exp.year.change = lm(dat_sta_age5$ExpLabScore ~ dat_sta_age5$ExpJulyDate, data = dat_sta_age5) dat_sta_age5$ExpLabScore.Corrected = dat_sta_age5$ExpLabScore - coef(exp.year.change)[2] #calculate sqrt of corrected Exploration score sqrtExpLabScore.Cor = sqrt(dat_sta_age5$ExpLabScore.Corrected) dat_sta_age5$sqrtExpLabScore.Cor = round(sqrtExpLabScore.Cor, digit = 2) # standardise corrected exploration score per Plot-Year PlotYearsqrtExpCorMean = aggregate(dat_sta_age5$sqrtExpLabScore.Cor, list(dat_sta_age5$PlotYear), mean) colnames(PlotYearsqrtExpCorMean) = c("PlotYear", "PlotYearsqrtExpCorMean") summary(PlotYearsqrtExpCorMean) PlotYearsqrtExpCorSD = aggregate(dat_sta_age5$sqrtExpLabScore.Cor, list(dat_sta_age5$PlotYear), sd) colnames(PlotYearsqrtExpCorSD) = c("PlotYear", "PlotYearsqrtExpCorSD") summary(PlotYearsqrtExpCorSD) dat_sta_age5 = merge(dat_sta_age5, PlotYearsqrtExpCorMean, by = "PlotYear") dat_sta_age5 = merge(dat_sta_age5, PlotYearsqrtExpCorSD, by = "PlotYear") dat_sta_age5$PYstdsqrtExpScoreCor = (dat_sta_age5$sqrtExpLabScore.Cor - dat_sta_age5$PlotYearsqrtExpCorMean)/dat_sta_age5$PlotYearsqrtExpCorSD # standardize corrected exploration score over the whole dataset dat_sta_age5$DSsqrtExpCorMean = mean(dat_sta_age5$sqrtExpLabScore.Cor) dat_sta_age5$DSsqrtExpCorSD = sd(dat_sta_age5$sqrtExpLabScore.Cor) dat_sta_age5$DSstdsqrtExpScoreCor = (dat_sta_age5$sqrtExpLabScore.Cor-dat_sta_age5$DSsqrtExpCorMean)/dat_sta_age5$DSsqrtExpCorSD ## Antwerp population---- dat_ant = read.table("Dataset_S2_Antwerp.csv", sep = ",", header = T) dat_ant_rec = read.table("Dataset_S2.2_LocalRecruits_Antwerp.csv", sep = ",", header = T) # Merge behavioral data with immigration status data dat_ant_rec = merge(dat_ant, dat_ant_rec[, c("BirdID", "Immigrant")], by = "BirdID") # Factorize age field and select only yearlings (age= 1) dat_ant_rec = as.factor(dat_ant_rec$Age) dat_ant_age5 = dat_ant_rec[dat_ant_rec$Age=="1",] # Add Population ID to dataset dat_ant_age5$PopulationID = rep(2, nrow(dat_ant_age5)) # Calculate exploration score corrected for time of year (July 1st as reference) # Exploration July date dat_ant_age5$expdate = as.Date(dat_ant_age5$ExpDate,"%d/%m/%Y") dat_ant_age5$JulyYear = ifelse (as.numeric(format(dat_ant_age5$expdate, "%m"))<7, as.numeric(format(dat_ant_age5$expdate, "%Y"))-1, as.numeric(format(dat_ant_age5$expdate, "%Y"))) dat_ant_age5$July1st = as.Date(paste(dat_ant_age5$JulyYear, "07", "01", sep = "-"), origin = "1899-12-30","%Y-%m-%d") dat_ant_age5$ExpJulyDate = dat_ant_age5$expdate - dat_ant_age5$July1st # Correct exploration score for time of year (July 1st as reference) exp.year.change = lm(dat_ant_age5$ExpLabScore ~ dat_ant_age5$ExpJulyDate, data = dat_ant_age5) dat_ant_age5$ExpLabScore.Corrected = dat_ant_age5$ExpLabScore - coef(exp.year.change)[2] #calculate sqrt of corrected Exploration score sqrtExpLabScore.Cor = sqrt(dat_ant_age5$ExpLabScore.Corrected) dat_ant_age5$sqrtExpLabScore.Cor = round(sqrtExpLabScore.Cor, digit = 2) # standardise corrected exploration score over the whole dataset dat_ant_age5$DSsqrtExpCorMean = mean(dat_ant_age5$sqrtExpLabScore.Cor) dat_ant_age5$DSsqrtExpCorSD = sd(dat_ant_age5$sqrtExpLabScore.Cor) dat_ant_age5$DSstdsqrtExpScoreCor = (dat_ant_age5$sqrtExpLabScore.Cor-dat_ant_age5$DSsqrtExpCorMean)/dat_ant_age5$DSsqrtExpCorSD ## Lauwersmeer population---- dat_lauw = read.table("Dataset_S3_Lauwersmeer.csv", sep = ",", header = T) dat_lauw_local = read.table("Dataset_S3.2_LocalRecruits_Lauwersmeer.csv", sep = ",", header = T) # Merge behavioral data with immigration status data dat_lauw_rec = merge (dat_lauw, dat_lauw_local[, c("BirdID", "Immigrant", "BirthYear")], by = "BirdID", all.x = TRUE) # Factorize age field and select only yearlings (age= 5) dat_lauw_rec$Age = factor(dat_lauw_rec$Age) dat_lauw_age5 = dat_lauw_rec[dat_lauw_rec$Age=="5",] # Add Population ID to dataset dat_lauw_age5$PopulationID = rep(3, nrow(dat_lauw_age5)) # Rename Plots for they do not have same values than in Seewiesen population NewPlotNames = c("22", "23", "24", "25", "26", "27", "28", "29", "30", "31", "32", "33") dat_lauw_age5$Plot = factor(as.numeric(dat_lauw_age5$Plot), labels = NewPlotNames) # Calculate exploration score corrected for time of year (July 1st as reference) # Exploration July date dat_lauw_age5$expdate = as.Date(dat_lauw_age5$ExpDate,"%d/%m/%Y") dat_lauw_age5$JulyYear = ifelse (as.numeric(format(dat_lauw_age5$expdate, "%m"))<7, as.numeric(format(dat_lauw_age5$expdate, "%Y"))-1, as.numeric(format(dat_lauw_age5$expdate, "%Y"))) dat_lauw_age5$July1st = as.Date(paste(dat_lauw_age5$JulyYear, "07", "01", sep = "-"), origin = "1899-12-30","%Y-%m-%d") dat_lauw_age5$ExpJulyDate = dat_lauw_age5$expdate - dat_lauw_age5$July1st # Correct exploration score for time of year (July 1st as reference) exp.year.change = lm(dat_lauw_age5$ExpLabScore ~ dat_lauw_age5$ExpJulyDate, data = dat_lauw_age5) dat_lauw_age5$ExpLabScore.Corrected = dat_lauw_age5$ExpLabScore - coef(exp.year.change)[2] #calculate sqrt of corrected Exploration score sqrtExpLabScore.Cor = sqrt(dat_lauw_age5$ExpLabScore.Corrected) dat_lauw_age5$sqrtExpLabScore.Cor = round(sqrtExpLabScore.Cor, digit = 2) # standardize corrected exploration score over the whole dataset dat_lauw_age5$DSsqrtExpCorMean = mean(dat_lauw_age5$sqrtExpLabScore.Cor) dat_lauw_age5$DSsqrtExpCorSD = sd(dat_lauw_age5$sqrtExpLabScore.Cor) dat_lauw_age5$DSstdsqrtExpScoreCor = (dat_lauw_age5$sqrtExpLabScore.Cor-dat_lauw_age5$DSsqrtExpCorMean)/dat_lauw_age5$DSsqrtExpCorSD ## Westerheide population---- dat_west = read.table("Dataset_S4_Westerheide.csv", sep = ",", header = T) dat_west_local = read.table("Dataset_S4.2_Westerheide.csv", sep = ",", header = T) # Merge behavioral data with immigration status data dat_west_rec = merge (dat_west, dat_west_local[, c("BirdID", "Immigrant", "BirthYear")], by = "BirdID", all.x=TRUE) # Select yearlings dat_west_age5 = subset (dat_west_rec, Age== "5") # Set Immigrant category as =1 (local recruits dataset contained only data for local recruits) dat_west_age5 = dat_west_age5 %>% replace_na(list(Immigrant = "1")) # Create Population ID column dat_west_age5$PopulationID = rep(4, nrow(dat_west_age5)) # Calculate exploration score corrected for time of year (July 1st as reference) # Exploration July date dat_west_age5$expdate = as.Date(dat_west_age5$ExpDate,"%d/%m/%Y") dat_west_age5$JulyYear = ifelse (as.numeric(format(dat_west_age5$expdate, "%m"))<7, as.numeric(format(dat_west_age5$expdate, "%Y"))-1, as.numeric(format(dat_west_age5$expdate, "%Y"))) dat_west_age5$July1st = as.Date(paste(dat_west_age5$JulyYear, "07", "01", sep = "-"), origin = "1899-12-30","%Y-%m-%d") dat_west_age5$ExpJulyDate = dat_west_age5$expdate - dat_west_age5$July1st # Correct exploration score for time of year (July 1st as reference) exp.year.change = lm(dat_west_age5$ExpLabScore ~ dat_west_age5$ExpJulyDate, data = dat_west_age5) dat_west_age5$ExpLabScore.Corrected = dat_west_age5$ExpLabScore - coef(exp.year.change)[2] #calculate sqrt of corrected Exploration score sqrtExpLabScore.Cor = sqrt(dat_west_age5$ExpLabScore.Corrected) dat_west_age5$sqrtExpLabScore.Cor = round(sqrtExpLabScore.Cor, digit = 2) # standardise corrected exploration score over the whole dataset dat_west_age5$DSsqrtExpCorMean = mean(dat_west_age5$sqrtExpLabScore.Cor) dat_west_age5$DSsqrtExpCorSD = sd(dat_west_age5$sqrtExpLabScore.Cor) dat_west_age5$DSstdsqrtExpScoreCor = (dat_west_age5$sqrtExpLabScore.Cor-dat_west_age5$DSsqrtExpCorMean)/dat_west_age5$DSsqrtExpCorSD # Wytham population---- dat_Wyt = read.table("Dataset_S5_Wytham.csv", sep = ",", header = T) dat_wyt_chicks = read.table("Dataset_S5.2_LocalRecruits_Wytham.csv", sep = ",", header = T) # Merge behavioral data with immigration status data dat_wyt_rec = merge(dat_wyt, dat_wyt_chicks[, c("BirdID", "BirthYear")], by = "BirdID", all.x =TRUE) # Select yearlings dat_wyt_age5 = subset (dat_wyt_rec, Age== "1") # Remove records where exploration score = 0 as should be 1 at minimum dat_wyt_age5 = subset (dat_wyt_rec, ExpLabScore!= 0) # Categorize birds as local recruits or immigrants dat_wyt_age5$Immigrant = ifelse(dat_wyt_age5$ChickYear!= "NA", 0) dat_wyt_age5 = dat_wyt_age5 %>% replace_na(list(Immigrant = "1")) dat_wyt_age5 = dat_wyt_age5 %>% replace_na(list(AnnualRecruits = 0)) # Create Population ID column dat_wyt_age5$PopulationID = rep(5, nrow(dat_wyt_age5)) # Calculate exploration score corrected for time of year (July 1st as reference) # Exploration July date dat_wyt_age5$expdate = as.Date(dat_wyt_age5$ExpDate,"%d/%m/%Y") dat_wyt_age5$JulyYear = ifelse (as.numeric(format(dat_wyt_age5$expdate, "%m"))<7, as.numeric(format(dat_wyt_age5$expdate, "%Y"))-1, as.numeric(format(dat_wyt_age5$expdate, "%Y"))) dat_wyt_age5$July1st = as.Date(paste(dat_wyt_age5$JulyYear, "07", "01", sep = "-"), origin = "1899-12-30","%Y-%m-%d") dat_wyt_age5$ExpJulyDate = dat_wyt_age5$expdate - dat_wyt_age5$July1st # Correct exploration score for time of year (July 1st as reference) exp.year.change = lm(dat_wyt_age5$ExpLabScore ~ dat_wyt_age5$ExpJulyDate, data = dat_wyt_age5) dat_wyt_age5$ExpLabScore.Corrected = dat_wyt_age5$ExpLabScore - coef(exp.year.change)[2] #calculate sqrt of corrected Exploration score sqrtExpLabScore.Cor = sqrt(dat_wyt_age5$ExpLabScore.Corrected) dat_wyt_age5$sqrtExpLabScore.Cor = round(sqrtExpLabScore.Cor, digit = 2) # standardise corrected exploration score over the whole dataset dat_wyt_age5$DSsqrtExpCorMean = mean(dat_wyt_age5$sqrtExpLabScore.Cor) dat_wyt_age5$DSsqrtExpCorSD = sd(dat_wyt_age5$sqrtExpLabScore.Cor) dat_wyt_age5$DSstdsqrtExpScoreCor = (dat_wyt_age5$sqrtExpLabScore.Cor-dat_wyt_age5$DSsqrtExpCorMean)/dat_wyt_age5$DSsqrtExpCorSD # Select columns required for analyses within each dataset Sta_Rec = dplyr::select(dat_sta_age5, BirdID, Plot, PlotYear, PopulationID, BroodYear, Sex, Immigrant, ExpLabScore, DSstdsqrtExpScoreCor, ExpDate, AnnualRecruits) Lauw_Rec = dplyr::select(dat_lauw_age5, BirdID, Plot, PlotYear, PopulationID, BroodYear, Sex, Immigrant, ExpLabScore, DSstdsqrtExpScoreCor, ExpDate, AnnualRecruits) Ant_Rec = dplyr::select(dat_ant_age5, BirdID, Plot, PlotYear, PopulationID, BroodYear, Sex, Immigrant, ExpLabScore, DSstdsqrtExpScoreCor, ExpDate, AnnualRecruits) West_Rec = dplyr::select(dat_west_age5, BirdID, Plot, PlotYear, PopulationID, BroodYear, Sex, Immigrant, ExpLabScore, DSstdsqrtExpScoreCor, ExpDate, AnnualRecruits) Wyt_Rec = dplyr::select(dat_wyt_age5, BirdID, Plot, PlotYear, PopulationID, BroodYear, Sex, Immigrant, ExpLabScore, DSstdsqrtExpScoreCor, ExpDate, AnnualRecruits) # Combine data from all populations into one dataset dat_allpop = rbind (Sta_Rec, Lauw_Rec, Ant_Rec, West_Rec, Wyt_Rec) # Create variable Population x Year dat_allpop$PopYear = paste(dat_allpop$PopulationID,dat_allpop$BroodYear, sep = "") # Factorize data" colFactor = c("BirdID", "Sex", "BroodYear", "Plot", "PlotYear", "PopYear", "PopulationID", "Immigrant") dat_allpop[,colFactor] = lapply(dat_allpop[,colFactor], factor) # Model---- sdExp = sd(dat_allpop$DSstdsqrtExpScoreCor) pcprior <- list(prec = list(prior="pc.prec", param = c((sdExp/0.31), 0.01))) formula_immi = DSstdsqrtExpScoreCor ~ Immigrant + f(py, model = "iid", hyper = pcprior) + f(p, model = "iid", hyper = pcprior) + f(yr, model = "iid", hyper = pcprior) + f(pop, model = "iid", hyper = pcprior) + f(popyr, model = "iid", hyper = pcprior) modINLA_immi = inla (formula_immi, data = dat_allpop, family = "gaussian", control.compute=list(config = TRUE), control.predictor=list(compute=T)) summary(modINLA_immi) # Population variance intercept tau.int.Pop= modINLA_immi$internal.marginals.hyperpar$"Log precision for pop" s.marginal.pop <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Pop) ranef.int.Pop = as.data.frame (inla.zmarginal(s.marginal.pop)) # Plot variance intercept tau.int.Plot= modINLA_immi$internal.marginals.hyperpar$"Log precision for p" s.marginal.p <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Plot) ranef.int.Plot = as.data.frame (inla.zmarginal(s.marginal.p)) # Year variance intercept tau.int.Yr = modINLA_immi$internal.marginals.hyperpar$"Log precision for yr" s.marginal.yr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Yr) ranef.int.Yr = as.data.frame (inla.zmarginal(s.marginal.yr)) # Population-year variance intercept tau.int.PopYr= modINLA_immi$internal.marginals.hyperpar$"Log precision for popyr" s.marginal.popyr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PopYr) ranef.int.PopYr = as.data.frame (inla.zmarginal(s.marginal.popyr)) # Plot-year variance intercept tau.int.PY= modINLA_immi$internal.marginals.hyperpar$"Log precision for py" s.marginal.py <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PY) ranef.int.PY = as.data.frame (inla.zmarginal(s.marginal.py)) # Residual tau.gausobs= modINLA_immi$internal.marginals.hyperpar$"Log precision for the Gaussian observations" s.marginal.gausobs <- inla.tmarginal(function(x) (1/exp(x)), tau.gausobs) ranef.GausObs = as.data.frame (inla.zmarginal(s.marginal.gausobs))