## Analyses of multiple nest box great tit populations---- ## to investigate heterogeneous selection on exploration behaviour at multiple spatiotemporal scales. ## Populations: # Starnberg, Bavaria (LMU Munich, Germany): Data from Niels Dingemanse and Alexia Mouchet; 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 Mathijs van Overveld; years 2006 to 2017, 9 plots; Data manager: Frank Adriaensen # Wytham (Oxford, UK): Data from John Quinn, Ella Cole, Allison Roth; years 2005 to 2016, 1 plot ## Author/Date: Alexia Mouchet (LMU Munich), April 2021 --------------------------------------------------------------------------------------------------------------------------------------------------- ### Libraries---- library (tidyverse) library(plyr) library(dplyr) library(INLA) library(ggplot2) library(ggthemes) ### Dataset and data preparation for all populations ## Starnberg population---- ## Data files---- dat_sta = read.table("Dataset_S1_Starnberg.csv", sep = ",", header = T) densdat_sta = read.table("Dataset_S1.1_Density_Starnberg.csv", sep = ";", header = T) ## Data preparation---- # Calculate density per hectare densdat_sta$DensHa = round((densdat_sta$Density)/densdat_sta$AreaHa, digits = 2) # Add Population ID to dataset dat_sta$PopulationID = rep(1, nrow(dat_sta)) # Merge behavioral and fitness data with density data alldat_sta = merge(dat_sta, densdat_sta, by = c("PlotYear", "Plot", "BroodYear")) # factorisation of variables alldat_sta[,1:7] = lapply(alldat_sta[,1:7], factor) # Remove data from plot-year 182010 because only one record alldat_sta = alldat_sta[alldat_sta$PlotYear != "182010",] alldat_sta$PlotYear = factor(alldat_sta$PlotYear) alldat_sta$PopulationID=factor(alldat_sta$PopulationID) # Give unique value for BroodYear variable so it is unique accross populations (because fitted as random factor) alldat_sta$PopYear = as.factor(paste(alldat_sta$PopulationID, alldat_sta$BroodYear, sep = "")) # Integrative fitness metrics (annual fitness) # with annual recruits alldat_sta$AnnualFitness = alldat_sta$Survived + (alldat_sta$AnnualRecruits/2) # with annual fledglings alldat_sta$AnnualFitnessFled = alldat_sta$Survived + (alldat_sta$AnnualFledged/2) ## Calculate exploration score corrected for time of year (July 1st as reference) # Exploration July date alldat_sta$expdate = as.Date(alldat_sta$ExpDate,"%d/%m/%Y") alldat_sta$JulyYear = ifelse (as.numeric(format(alldat_sta$expdate, "%m"))<7, as.numeric(format(alldat_sta$expdate, "%Y"))-1, as.numeric(format(alldat_sta$expdate, "%Y"))) alldat_sta$July1st = as.Date(paste(alldat_sta$JulyYear, "07", "01", sep = "-"), origin = "1899-12-30","%Y-%m-%d") alldat_sta$ExpJulyDate = alldat_sta$expdate - alldat_sta$July1st ## Correct exploration score for time of year (July 1st as reference) exp.year.change = lm(alldat_sta$ExpLabScore ~ alldat_sta$ExpJulyDate, data = alldat_sta) alldat_sta$ExpLabScore.Corrected = alldat_sta$ExpLabScore - coef(exp.year.change)[2] #calculate sqrt of corrected Exploration score sqrtExpLabScore.Cor = sqrt(alldat_sta$ExpLabScore.Corrected) alldat_sta$sqrtExpLabScore.Cor = round(sqrtExpLabScore.Cor, digit = 2) # standardise corrected exploration score per Plot-Year PlotYearsqrtExpCorMean = aggregate(alldat_sta$sqrtExpLabScore.Cor, list(alldat_sta$PlotYear), mean) colnames(PlotYearsqrtExpCorMean) = c("PlotYear", "PlotYearsqrtExpCorMean") summary(PlotYearsqrtExpCorMean) PlotYearsqrtExpCorSD = aggregate(alldat_sta$sqrtExpLabScore.Cor, list(alldat_sta$PlotYear), sd) colnames(PlotYearsqrtExpCorSD) = c("PlotYear", "PlotYearsqrtExpCorSD") summary(PlotYearsqrtExpCorSD) # Add new variables to dataset alldat_sta = merge(alldat_sta, PlotYearsqrtExpCorMean, by = "PlotYear") alldat_sta = merge(alldat_sta, PlotYearsqrtExpCorSD, by = "PlotYear") alldat_sta$PYstdsqrtExpScoreCor = (alldat_sta$sqrtExpLabScore.Cor - alldat_sta$PlotYearsqrtExpCorMean)/alldat_sta$PlotYearsqrtExpCorSD # standardise corrected exploration score over the whole dataset alldat_sta$DSsqrtExpCorMean = mean(alldat_sta$sqrtExpLabScore.Cor) alldat_sta$DSsqrtExpCorSD = sd(alldat_sta$sqrtExpLabScore.Cor) alldat_sta$DSstdsqrtExpScoreCor = (alldat_sta$sqrtExpLabScore.Cor-alldat_sta$DSsqrtExpCorMean)/alldat_sta$DSsqrtExpCorSD # center density over whole dataset alldat_sta$densDSMean = mean(alldat_sta$DensHa) alldat_sta$densCenDSMean = alldat_sta$DensHa - alldat_sta$densDSMean ## Relative integrative fitness---- # Mean integrative fitness (calculated with recruits) per Plot-Year PlotYearFitnessMean = aggregate(alldat_sta$AnnualFitness, list(alldat_sta$PlotYear), mean) colnames(PlotYearFitnessMean) = c("PlotYear", "PYFitnessMean") alldat_sta = merge(alldat_sta, PlotYearFitnessMean, by = "PlotYear") # Calculate relative integrative fitness by dividing by mean integrative fitness alldat_sta$wAnnualFitness = alldat_sta$AnnualFitness/alldat_sta$PYFitnessMean ## Relative Recruitment---- # Mean recruitment per Plot-Year PlotYearRecruitMean = aggregate(alldat_sta$AnnualRecruits, list(alldat_sta$PlotYear), mean) colnames(PlotYearRecruitMean) = c("PlotYear", "PYRecruitMean") alldat_sta = merge(alldat_sta, PlotYearRecruitMean, by = "PlotYear") # Calculate relative recruitment by dividing by mean recruitment alldat_sta$wRecruit = alldat_sta$AnnualRecruits/alldat_sta$PYRecruitMean ## Relative Survival---- # Mean survival per Plot-Year PlotYearSurvMean = aggregate(alldat_sta$Survived, list(alldat_sta$PlotYear), mean) colnames(PlotYearSurvMean) = c("PlotYear", "PYSurvMean") alldat_sta = merge(alldat_sta, PlotYearSurvMean, by = "PlotYear") # Calculate relative survival by dividing by mean survival alldat_sta$wSurvival = alldat_sta$Survived/alldat_sta$PYSurvMean ## Antwerp population---- ## Data files---- dat_ant = read.table("Dataset_S2_Antwerp.csv", sep = ",", header = T) densdat_ant = read.table("Dataset_S2.1_Density_Antwerp.csv", sep = ";", header = T) ## Data preparation---- # Data selection and new variables # remove last year of monitoring because fitness data not known for that year densdat_ant = subset(densdat_ant, BroodYear < 2018) # calculate density per hectare densdat_ant$DensHa = round((densdat_ant$Density)/densdat_ant$AreaHa, digits = 2) # Add Population ID dat_ant$PopulationID=rep(2, nrow(dat_ant)) # Merge behavioral and fitness data with density data alldat_ant = merge(dat_ant, densdat_ant, by = c("BroodYear", "PlotYear", "Plot")) # factorisation of variables alldat_ant[,1:7] = lapply(alldat_ant[,1:7], factor) alldat_ant$PopulationID=factor(alldat_ant$PopulationID) # Give BroodYear a new value for it is unique across populations (because fitted as random effect) alldat_ant$PopYear = as.factor(paste(alldat_ant$PopulationID, alldat_ant$BroodYear, sep = "")) ## Integrative fitness metrics (annual fitness) # with annual recruits alldat_ant$AnnualFitness = alldat_ant$Survived + (alldat_ant$AnnualRecruits/2) # with annual fledglings alldat_ant$AnnualFitnessFled = alldat_ant$Survived + (alldat_ant$AnnualFledged/2) ## Calculate exploration score corrected for time of year (July 1st as reference) # Exploration July date alldat_ant$expdate = as.Date(alldat_ant$ExpDate,"%d/%m/%Y") alldat_ant$JulyYear = ifelse (as.numeric(format(alldat_ant$expdate, "%m"))<7, as.numeric(format(alldat_ant$expdate, "%Y"))-1, as.numeric(format(alldat_ant$expdate, "%Y"))) alldat_ant$July1st = as.Date(paste(alldat_ant$JulyYear, "07", "01", sep = "-"), origin = "1899-12-30","%Y-%m-%d") alldat_ant$ExpJulyDate = alldat_ant$expdate - alldat_ant$July1st ## Correct exploration score for time of year (July 1st as reference) exp.year.change = lm(alldat_ant$ExpLabScore ~ alldat_ant$ExpJulyDate, data = alldat_ant) alldat_ant$ExpLabScore.Corrected = alldat_ant$ExpLabScore - coef(exp.year.change)[2] #calculate sqrt of corrected Exploration score sqrtExpLabScore.Cor = sqrt(alldat_ant$ExpLabScore.Corrected) alldat_ant$sqrtExpLabScore.Cor = round(sqrtExpLabScore.Cor, digit = 2) # standardise corrected exploration score per PlotYear PlotYearsqrtExpCorMean = aggregate(alldat_ant$sqrtExpLabScore.Cor, list(alldat_ant$PlotYear), mean) colnames(PlotYearsqrtExpCorMean) = c("PlotYear", "PlotYearsqrtExpCorMean") summary(PlotYearsqrtExpCorMean) PlotYearsqrtExpCorSD = aggregate(alldat_ant$sqrtExpLabScore.Cor, list(alldat_ant$PlotYear), sd) colnames(PlotYearsqrtExpCorSD) = c("PlotYear", "PlotYearsqrtExpCorSD") summary(PlotYearsqrtExpCorSD) # Add new variables to dataset alldat_ant = merge(alldat_ant, PlotYearsqrtExpCorMean, by = "PlotYear") alldat_ant = merge(alldat_ant, PlotYearsqrtExpCorSD, by = "PlotYear") alldat_ant$PYstdsqrtExpScoreCor = (alldat_ant$sqrtExpLabScore.Cor - alldat_ant$PlotYearsqrtExpCorMean)/alldat_ant$PlotYearsqrtExpCorSD # standardize corrected exploration score over the whole dataset alldat_ant$DSsqrtExpCorMean = mean(alldat_ant$sqrtExpLabScore.Cor) alldat_ant$DSsqrtExpCorSD = sd(alldat_ant$sqrtExpLabScore.Cor) alldat_ant$DSstdsqrtExpScoreCor = (alldat_ant$sqrtExpLabScore.Cor-alldat_ant$DSsqrtExpCorMean)/alldat_ant$DSsqrtExpCorSD # center density over whole dataset alldat_ant$densDSMean = mean(alldat_ant$DensHa) alldat_ant$densCenDSMean = alldat_ant$DensHa-alldat_ant$densDSMean ## Relative integrative fitness---- # Mean integrative fitness per Plot-Year PlotYearFitnessMean = aggregate(alldat_ant$AnnualFitness, list(alldat_ant$PlotYear), mean) colnames(PlotYearFitnessMean) = c("PlotYear", "PYFitnessMean") alldat_ant = merge(alldat_ant, PlotYearFitnessMean, by = "PlotYear") # Calculate relative integrative fitness by dividing by mean integrative fitness alldat_ant$wAnnualFitness = alldat_ant$AnnualFitness/alldat_ant$PYFitnessMean ## Relative recruitment---- # Mean recruitment per Plot-Year PlotYearRecruitMean = aggregate(alldat_ant$AnnualRecruits, list(alldat_ant$PlotYear), mean) colnames(PlotYearRecruitMean) = c("PlotYear", "PYRecruitMean") alldat_ant = merge(alldat_ant, PlotYearRecruitMean, by = "PlotYear") # Calculate relative recruitment by dividing by mean recruitment alldat_ant$wRecruit = alldat_ant$AnnualRecruits/alldat_ant$PYRecruitMean ## Relative survival---- # Mean survival per Plot-Year PlotYearSurvMean = aggregate(alldat_ant$Survived, list(alldat_ant$PlotYear), mean) colnames(PlotYearSurvMean) = c("PlotYear", "PYSurvMean") alldat_ant = merge(alldat_ant, PlotYearSurvMean, by = "PlotYear") alldat_ant$wSurvival = alldat_ant$Survived/alldat_ant$PYSurvMean summary(alldat_ant) ## Lauwersmeer population---- ## Data files---- dat_lauw = read.table("Dataset_S3_Lauwersmeer_MS.csv", sep = ",", header = T) densdat_lauw = read.table("Dataset_S3.1_Density_Lauwersmeer.csv", sep = ";", header = T) ## Data preparation---- # Calculate density per hectare densdat_lauw$DensHa = densdat_lauw$Density/densdat_lauw$Area # Merge behavioural and fitness data with density data alldat_lauw = merge(dat_lauw, densdat_lauw, by = c("PlotYear", "Plot", "BroodYear")) # factorisation of variables alldat_lauw[,1:7 ]= lapply(alldat_lauw[,1:7], factor) # Add Population ID alldat_lauw$PopulationID=rep(3, nrow(alldat_lauw)) # Rename Plots for they do not have same values than in Starnberg population NewPlotNames = c("22", "23", "24", "25", "26", "27", "28", "29", "30", "31", "32", "33") alldat_lauw$Plot = factor(as.numeric(alldat_lauw$Plot), labels = NewPlotNames) # Rename Years for they also have unique values across populations alldat_lauw$PopYear = as.factor(paste(alldat_lauw$PopulationID, alldat_lauw$BroodYear, sep = "")) ## Integrative fitness metrics (annual fitness) # with annual recruits alldat_lauw$AnnualFitness = alldat_lauw$Survived + (alldat_lauw$AnnualRecruits/2) # with annual fledglings alldat_lauw$AnnualFitnessFled = alldat_lauw$Survived + (alldat_lauw$AnnualFledged/2) ## Calculate exploration score corrected for time of year (July 1st as reference) # Exploration July date alldat_lauw$expdate = as.Date(alldat_lauw$ExpDate,"%d/%m/%Y") alldat_lauw$JulyYear = ifelse (as.numeric(format(alldat_lauw$expdate, "%m"))<7, as.numeric(format(alldat_lauw$expdate, "%Y"))-1, as.numeric(format(alldat_lauw$expdate, "%Y"))) alldat_lauw$July1st = as.Date(paste(alldat_lauw$JulyYear, "07", "01", sep = "-"), origin = "1899-12-30","%Y-%m-%d") alldat_lauw$ExpJulyDate = alldat_lauw$expdate - alldat_lauw$July1st ## Correct exploration score for time of year (July 1st as reference) exp.year.change = lm(alldat_lauw$ExpLabScore ~ alldat_lauw$ExpJulyDate, data = alldat_lauw) alldat_lauw$ExpLabScore.Corrected = alldat_lauw$ExpLabScore - coef(exp.year.change)[2] #calculate sqrt of corrected Exploration score sqrtExpLabScore.Cor = sqrt(alldat_lauw$ExpLabScore.Corrected) alldat_lauw$sqrtExpLabScore.Cor = round(sqrtExpLabScore.Cor, digit = 2) # standardise corrected exploration score per PlotYear PlotYearsqrtExpCorMean = aggregate(alldat_lauw$sqrtExpLabScore.Cor, list(alldat_lauw$PlotYear), mean) colnames(PlotYearsqrtExpCorMean) = c("PlotYear","PlotYearsqrtExpCorMean") summary(PlotYearsqrtExpCorMean) PlotYearsqrtExpCorSD = aggregate(alldat_lauw$sqrtExpLabScore.Cor, list(alldat_lauw$PlotYear), sd) colnames(PlotYearsqrtExpCorSD) = c("PlotYear","PlotYearsqrtExpCorSD") summary(PlotYearsqrtExpCorSD) alldat_lauw = merge(alldat_lauw, PlotYearsqrtExpCorMean, by = "PlotYear") alldat_lauw = merge(alldat_lauw, PlotYearsqrtExpCorSD, by = "PlotYear") alldat_lauw$PYstdsqrtExpScoreCor = (alldat_lauw$sqrtExpLabScore.Cor - alldat_lauw$PlotYearsqrtExpCorMean)/alldat_lauw$PlotYearsqrtExpCorSD # standardize corrected exploration score over the whole dataset alldat_lauw$DSsqrtExpCorMean = mean(alldat_lauw$sqrtExpLabScore.Cor) alldat_lauw$DSsqrtExpCorSD = sd(alldat_lauw$sqrtExpLabScore.Cor) alldat_lauw$DSstdsqrtExpScoreCor = (alldat_lauw$sqrtExpLabScore.Cor-alldat_lauw$DSsqrtExpCorMean)/alldat_lauw$DSsqrtExpCorSD # center density over whole dataset alldat_lauw$densDSMean = mean(alldat_lauw$DensHa) alldat_lauw$densCenDSMean = alldat_lauw$DensHa - alldat_lauw$densDSMean ## Relative integrative fitness---- # Mean integrative fitness (calculated with recruits) per Plot-Year PlotYearFitnessMean = aggregate(alldat_lauw$AnnualFitness, list(alldat_lauw$PlotYear), mean) colnames(PlotYearFitnessMean) = c("PlotYear", "PYFitnessMean") alldat_lauw = merge(alldat_lauw, PlotYearFitnessMean, by = "PlotYear") # Calculate relative integrative fitness by dividing by mean integrative fitness alldat_lauw$wAnnualFitness = alldat_lauw$AnnualFitness/alldat_lauw$PYFitnessMean ## Relative recruitment---- # Mean recruitment per Plot-Year PlotYearRecruitMean = aggregate(alldat_lauw$AnnualRecruits, list(alldat_lauw$PlotYear), mean) colnames(PlotYearRecruitMean) = c("PlotYear", "PYRecruitMean") alldat_lauw = merge(alldat_lauw, PlotYearRecruitMean, by = "PlotYear") # Calculate relative recruitment by dividing by mean recruitment alldat_lauw$wRecruit = alldat_lauw$AnnualRecruits/alldat_lauw$PYRecruitMean ## Relative survival---- # Mean survival per Plot-Year PlotYearSurvMean = aggregate(alldat_lauw$Survived, list(alldat_lauw$PlotYear), mean) colnames(PlotYearSurvMean) = c("PlotYear", "PYSurvMean") alldat_lauw = merge(alldat_lauw, PlotYearSurvMean, by = "PlotYear") # Calculate relative survival by dividing by mean survival alldat_lauw$wSurvival = alldat_lauw$Survived/alldat_lauw$PYSurvMean summary (alldat_lauw) ## Westerheide population---- ## Data files---- dat_west = read.table("Dataset_S4_Westerheide.csv", sep=",", header=T) densdat_west = read.table("Dataset_S4.1_Density_Westerheide.csv", sep=";", header=T) ## Data preparation---- # Data selection: remove exploration scores of 0 as should be 1 as minimum dat_west = subset(dat_west, dat_west$ExpLabScore != 0) # Add Population ID dat_west$PopulationID=rep(4, nrow(dat_west)) # Merge behavioral and fitness data with density data alldat_west = merge(dat_west, densdat_west, by = c("PlotYear", "BroodYear")) # factorisation of variables alldat_west[,1:7] = lapply(alldat_west[,1:7], factor) alldat_west$PopulationID=factor(alldat_west$PopulationID) # Give BroodYear a new value for it is unique accross populations (because fitted as random effect) alldat_west$PopYear = as.factor(paste(alldat_west$PopulationID, alldat_west$BroodYear, sep = "")) # Give unique ID for plot-year as only one plot in this population alldat_west$PlotYear = rep(98, nrow(alldat_west)) alldat_west$PlotYear = as.factor(alldat_west$PlotYear) ## Integrative fitness metrics (annual fitness) # with annual recruits alldat_west$AnnualFitness = alldat_west$Survived + (alldat_west$AnnualRecruits/2) # with annual fledglings alldat_west$AnnualFitnessFled = alldat_west$Survived + (alldat_west$AnnualFledged/2) ## Calculate exploration score corrected for time of year (July 1st as reference) # Exploration July date alldat_west$expdate = as.Date(alldat_west$ExpDate,"%d/%m/%Y") alldat_west$JulyYear = ifelse (as.numeric(format(alldat_west$expdate, "%m"))<7, as.numeric(format(alldat_west$expdate, "%Y"))-1, as.numeric(format(alldat_west$expdate, "%Y"))) alldat_west$July1st = as.Date(paste(alldat_west$JulyYear, "07", "01", sep = "-"), origin = "1899-12-30","%Y-%m-%d") alldat_west$ExpJulyDate = alldat_west$expdate - alldat_west$July1st ## Correct exploration score for time of year (July 1st as reference) exp.year.change = lm(alldat_west$ExpLabScore ~ alldat_west$ExpJulyDate, data = alldat_west) alldat_west$ExpLabScore.Corrected = alldat_west$ExpLabScore - coef(exp.year.change)[2] #calculate sqrt of corrected Exploration score sqrtExpLabScore.Cor = sqrt(alldat_west$ExpLabScore.Corrected) alldat_west$sqrtExpLabScore.Cor = round(sqrtExpLabScore.Cor, digit = 2) # standardise corrected exploration score per alldat_west PlotYearsqrtExpCorMean = aggregate(alldat_west$sqrtExpLabScore.Cor, list(alldat_west$PlotYear), mean) colnames(PlotYearsqrtExpCorMean) = c("PlotYear", "PlotYearsqrtExpCorMean") summary(PlotYearsqrtExpCorMean) PlotYearsqrtExpCorSD = aggregate(alldat_west$sqrtExpLabScore.Cor, list(alldat_west$PlotYear), sd) colnames(PlotYearsqrtExpCorSD) = c("PlotYear", "PlotYearsqrtExpCorSD") summary(PlotYearsqrtExpCorSD) # Add new variables to dataset alldat_west = merge(alldat_west, PlotYearsqrtExpCorMean, by = "PlotYear") alldat_west = merge(alldat_west, PlotYearsqrtExpCorSD, by = "PlotYear") alldat_west$PYstdsqrtExpScoreCor = (alldat_west$sqrtExpLabScore.Cor - alldat_west$PlotYearsqrtExpCorMean)/alldat_west$PlotYearsqrtExpCorSD # standardize corrected exploration score over the whole dataset alldat_west$DSsqrtExpCorMean = mean(alldat_west$sqrtExpLabScore.Cor) alldat_west$DSsqrtExpCorSD = sd(alldat_west$sqrtExpLabScore.Cor) alldat_west$DSstdsqrtExpScoreCor = (alldat_west$sqrtExpLabScore.Cor-alldat_west$DSsqrtExpCorMean)/alldat_west$DSsqrtExpCorSD # center density over whole dataset alldat_west$densDSMean=mean(alldat_west$DensHa) alldat_west$densCenDSMean= alldat_west$DensHa-alldat_west$densDSMean alldat_west$YearMeanDeviation = alldat_west$densCenDSMean ## Relative integrative fitness---- # Mean integrative fitness (calculated with recruits) per Plot-Year PlotYearFitnessMean = aggregate(alldat_west$AnnualFitness, list(alldat_west$PlotYear), mean) colnames(PlotYearFitnessMean) = c("PlotYear", "PYFitnessMean") alldat_west = merge(alldat_west, PlotYearFitnessMean, by = "PlotYear") # Calculate relative integrative fitness by dividing by mean integrative fitness alldat_west$wAnnualFitness = alldat_west$AnnualFitness/alldat_west$PYFitnessMean ## Relative recruitment---- # Mean recruitment per PlotYear PlotYearRecruitMean = aggregate(alldat_west$AnnualRecruits, list(alldat_west$PlotYear), mean) colnames(PlotYearRecruitMean) = c("PlotYear", "PYRecruitMean") alldat_west = merge(alldat_west, PlotYearRecruitMean, by = "PlotYear") # Calculate relative recruitment by dividing by mean recruitment alldat_west$wRecruit = alldat_west$AnnualRecruits/alldat_west$PYRecruitMean ## Relative survival---- # Mean survival per PlotYear PlotYearSurvMean = aggregate(alldat_west$Survived, list(alldat_west$PlotYear), mean) colnames(PlotYearSurvMean) = c("PlotYear", "PYSurvMean") alldat_west = merge(alldat_west, PlotYearSurvMean, by = "PlotYear") # Calculate relative survival by dividing by mean survival alldat_west$wSurvival = alldat_west$Survived/alldat_west$PYSurvMean summary(alldat_west) ## Wytham population---- ## Data files---- dat_wyt = read.table("Dataset_S5_Wytham.csv", sep=",", header=T) densdat_wyt = read.table("Dataset_S5.1_Density_Wytham.csv", sep=";", header=T) ## Data preparation---- # Calculate density per hectare densdat_wyt$DensHa = round((densdat_wyt$Density/densdat_wyt$AreaHa),2) # Add Population ID dat_wyt$PopulationID=rep(5, nrow(dat_wyt)) # Merge behavioral and fitness data with density data alldat_wyt = merge(dat_wyt, densdat_wyt, by = c("PlotYear", "BroodYear")) # factorisation of variables alldat_wyt[,1:7] = lapply(alldat_wyt[,1:7], factor) alldat_wyt$PopulationID = factor(alldat_wyt$PopulationID) # Remove exploration scores = 0 as they should not exist (minimum should be 1) alldat_wyt = subset(alldat_wyt, alldat_wyt$ExpLabScore != 0) # Give BroodYear a new value for it is unique accorss populations (because fitted as random effect) alldat_wyt$PopYear = as.factor(paste(alldat_wyt$PopulationID, alldat_wyt$BroodYear, sep = "")) # Create Plot-Year variable - same ID for all records because only one plot, # and for it does not contribute to Plot-Year variance estimation alldat_wyt$PlotYear = rep(99, nrow(alldat_wyt)) alldat_wyt$PlotYear = as.factor(alldat_wyt$PlotYear) ## Integrative fitness metrics (annual fitness) # with annual recruits alldat_wyt$AnnualFitness = alldat_wyt$Survived + (alldat_wyt$AnnualRecruits/2) # with annual fledglings alldat_wyt$AnnualFitnessFled = alldat_wyt$Survived + (alldat_wyt$AnnualFledged/2) ## Calculate exploration score corrected for time of year (July 1st as reference) # Exploration July date alldat_wyt$expdate = as.Date(alldat_wyt$ExpDate,"%d/%m/%Y") alldat_wyt$JulyYear = ifelse (as.numeric(format(alldat_wyt$expdate, "%m"))<7, as.numeric(format(alldat_wyt$expdate, "%Y"))-1, as.numeric(format(alldat_wyt$expdate, "%Y"))) alldat_wyt$July1st = as.Date(paste(alldat_wyt$JulyYear, "07", "01", sep = "-"), origin = "1899-12-30","%Y-%m-%d") alldat_wyt$ExpJulyDate = alldat_wyt$expdate - alldat_wyt$July1st ## Correct exploration score for time of year (July 1st as reference) exp.year.change = lm(alldat_wyt$ExpLabScore ~ alldat_wyt$ExpJulyDate, data = alldat_wyt) alldat_wyt$ExpLabScore.Corrected = alldat_wyt$ExpLabScore - coef(exp.year.change)[2] #calculate sqrt of corrected Exploration score sqrtExpLabScore.Cor = sqrt(alldat_wyt$ExpLabScore.Corrected) alldat_wyt$sqrtExpLabScore.Cor = round(sqrtExpLabScore.Cor, digit = 2) # standardise corrected exploration score per PlotYear PlotYearsqrtExpCorMean = aggregate(alldat_wyt$sqrtExpLabScore.Cor, list(alldat_wyt$PlotYear), mean) colnames(PlotYearsqrtExpCorMean) = c("PlotYear", "PlotYearsqrtExpCorMean") summary(PlotYearsqrtExpCorMean) PlotYearsqrtExpCorSD = aggregate(alldat_wyt$sqrtExpLabScore.Cor, list(alldat_wyt$PlotYear), sd) colnames(PlotYearsqrtExpCorSD) = c("PlotYear", "PlotYearsqrtExpCorSD") summary(PlotYearsqrtExpCorSD) # Add new variables to dataset alldat_wyt = merge(alldat_wyt, PlotYearsqrtExpCorMean, by = "PlotYear") alldat_wyt = merge(alldat_wyt, PlotYearsqrtExpCorSD, by = "PlotYear") alldat_wyt$PYstdsqrtExpScoreCor = (alldat_wyt$sqrtExpLabScore.Cor - alldat_wyt$PlotYearsqrtExpCorMean)/alldat_wyt$PlotYearsqrtExpCorSD # standardize corrected exploration score over the whole dataset alldat_wyt$DSsqrtExpCorMean = mean(alldat_wyt$sqrtExpLabScore.Cor) alldat_wyt$DSsqrtExpCorSD = sd(alldat_wyt$sqrtExpLabScore.Cor) alldat_wyt$DSstdsqrtExpScoreCor = (alldat_wyt$sqrtExpLabScore.Cor-alldat_wyt$DSsqrtExpCorMean)/alldat_wyt$DSsqrtExpCorSD # center density over whole dataset alldat_wyt$densDSMean=mean(alldat_wyt$DensHa) alldat_wyt$densCenDSMean= alldat_wyt$DensHa-alldat_wyt$densDSMean alldat_wyt$YearMeanDeviation = alldat_wyt$densCenDSMean ## Relative integrative fitness---- # Mean integrative fitness (calculated with recruits) per Plot-Year PlotYearFitnessMean = aggregate(alldat_wyt$AnnualFitness, list(alldat_wyt$PlotYear), mean) colnames(PlotYearFitnessMean) = c("PlotYear", "PYFitnessMean") alldat_wyt = merge(alldat_wyt, PlotYearFitnessMean, by = "PlotYear") # Calculate relative integrative fitness by dividing by mean integrative fitness alldat_wyt$wAnnualFitness = alldat_wyt$AnnualFitness/alldat_wyt$PYFitnessMean ## Relative recruitment---- # Mean recruitment per Plot-Year PlotYearRecruitMean = aggregate(alldat_wyt$AnnualRecruits, list(alldat_wyt$PlotYear), mean) colnames(PlotYearRecruitMean) = c("PlotYear", "PYRecruitMean") alldat_wyt = merge(alldat_wyt, PlotYearRecruitMean, by = "PlotYear") # Calculate relative recruitment by dividing by mean recruitment alldat_wyt$wRecruit = alldat_wyt$AnnualRecruits/alldat_wyt$PYRecruitMean ## Relative Survival---- # Mean survival per Plot-Year PlotYearSurvMean = aggregate(alldat_wyt$Survived, list(alldat_wyt$PlotYear), mean) colnames(PlotYearSurvMean) = c("PlotYear", "PYSurvMean") alldat_wyt = merge(alldat_wyt, PlotYearSurvMean, by = "PlotYear") # Calculate relative survival by dividing by mean survival alldat_wyt$wSurvival = alldat_wyt$Survived/alldat_wyt$PYSurvMean summary(alldat_wyt) ## Populations combined in one DS - Standardisation and centering within population---- Sta = dplyr::select(alldat_sta, PlotYear, BroodYear, PopYear, Plot, PopulationID, PYstdsqrtExpScoreCor, DSstdsqrtExpScoreCor, DensHa, densDSMean, densCenDSMean, AnnualRecruits, AnnualFledged, Survived, AnnualFitness, AnnualFitnessFled, wRecruit, wAnnualFitness, wSurvival) Lauw = dplyr::select(alldat_lauw, PlotYear, BroodYear, PopYear, Plot, PopulationID, PYstdsqrtExpScoreCor, DSstdsqrtExpScoreCor, DensHa, densDSMean, densCenDSMean, AnnualRecruits, AnnualFledged, Survived, AnnualFitness, AnnualFitnessFled, wRecruit, wAnnualFitness, wSurvival) Antwerp = dplyr::select(alldat_ant, PlotYear, BroodYear, PopYear, Plot, PopulationID, PYstdsqrtExpScoreCor, DSstdsqrtExpScoreCor, DensHa, densDSMean, densCenDSMean, AnnualRecruits, AnnualFledged, Survived, AnnualFitness, AnnualFitnessFled, wRecruit, wAnnualFitness, wSurvival) West = dplyr::select(alldat_west, PlotYear, BroodYear, PopYear, Plot, PopulationID, PYstdsqrtExpScoreCor, DSstdsqrtExpScoreCor, DensHa, densDSMean, densCenDSMean, AnnualRecruits, AnnualFledged, Survived, AnnualFitness, AnnualFitnessFled, wRecruit, wAnnualFitness, wSurvival) Wytham = dplyr::select(alldat_wyt, PlotYear, BroodYear, PopYear, Plot, PopulationID, PYstdsqrtExpScoreCor, DSstdsqrtExpScoreCor, DensHa, densDSMean, densCenDSMean, AnnualRecruits, AnnualFledged, Survived, AnnualFitness, AnnualFitnessFled, wRecruit, wAnnualFitness, wSurvival) datPopComb = rbind.fill(Sta, Lauw, West, Antwerp, Wytham) # center density over the dataset with all populations to investigate among-population differences in the effect of density on survival datPopComb$MeanPopDensMeans = mean (datPopComb$densDSMean) datPopComb$PopDensCen = datPopComb$densDSMean - datPopComb$MeanPopDensMeans # To check for overdispersion: to be added as random effect in poisson model; if variance = 0 means no overdispersion datPopComb$obsid = 1:nrow(datPopComb) summary(datPopComb) # Variable definition for iid2d INLA models---- pop = as.numeric(datPopComb$PopulationID) n.pop = max(pop) N.Pop = 2*n.pop pop2 = as.numeric(datPopComb$PopulationID) + n.pop n.plot = nlevels(datPopComb$Plot) N.Plot = 2*n.plot p = as.numeric(datPopComb$Plot) p2 = as.numeric(datPopComb$Plot) + n.plot n.year = nlevels(datPopComb$BroodYear) N.Year = 2*n.year yr = as.numeric(datPopComb$BroodYear) yr2 = as.numeric(datPopComb$BroodYear) + n.year n.popyear = nlevels(datPopComb$PopYear) N.PopYear = 2*n.popyear popyr = as.numeric(datPopComb$PopYear) popyr2 = as.numeric(datPopComb$PopYear) + n.popyear n.PY = nlevels(datPopComb$PlotYear) N.PY = 2*n.PY py = as.numeric(datPopComb$PlotYear) py2 = as.numeric(datPopComb$PlotYear) + n.PY datPopComb$pop = pop datPopComb$pop2 = pop2 datPopComb$p = p datPopComb$p2 = p2 datPopComb$yr = yr datPopComb$yr2 = yr2 datPopComb$popyr = popyr datPopComb$popyr2 = popyr2 datPopComb$py = py datPopComb$py2 = py2 #### Analyses for results presented in Tables (main text and SI) and graph for figure 2---- ### Table 1: Standardized selection gradient (exploration behavior standardized over dataset)---- ## Integrative fitness---- # Model formula_IndFit_DSstd_w = wAnnualFitness ~ DSstdsqrtExpScoreCor + I(DSstdsqrtExpScoreCor^2) + f(pop, model = "iid2d", n=N.Pop) + f(pop2, DSstdsqrtExpScoreCor, copy = "pop") + f(p, model = "iid2d", n=N.Plot) + f(p2, DSstdsqrtExpScoreCor, copy = "p") + f(yr, model = "iid2d", n=N.Year) + f(yr2, DSstdsqrtExpScoreCor, copy = "yr") + f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, DSstdsqrtExpScoreCor, copy = "popyr") + f(py, model = "iid2d", n=N.PY) + f(py2, DSstdsqrtExpScoreCor, copy = "py") modINLA_IndFit_DSstd_w = inla (formula_IndFit_DSstd_w, data = datPopComb, family = "gaussian", control.compute=list(cpo = TRUE, config = TRUE), control.predictor=list(compute=T)) summary(modINLA_IndFit_DSstd_w) # Linear and nonlinear selection gradient fixefTable = modINLA_IndFit_DSstd_w$summary.fixed ## Offspring recruitment---- # Model formula_wRec = wRecruit ~ DSstdsqrtExpScoreCor + I(DSstdsqrtExpScoreCor^2) + f(pop, model = "iid2d", n=N.Pop) + f(pop2, DSstdsqrtExpScoreCor, copy = "pop") + f(p, model = "iid2d", n=N.Plot) + f(p2, DSstdsqrtExpScoreCor, copy = "p") + f(yr, model = "iid2d", n=N.Year) + f(yr2, DSstdsqrtExpScoreCor, copy = "yr") + f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, DSstdsqrtExpScoreCor, copy = "popyr") + f(py, model = "iid2d", n=N.PY) + f(py2, DSstdsqrtExpScoreCor, copy = "py") modINLA_wRec = inla (formula_wRec, data = datPopComb, family = "gaussian", control.compute=list(cpo = TRUE, config = TRUE), control.predictor=list(compute=T)) summary(modINLA_wRec) # Linear and nonlinear selection gradient fixefTable = modINLA_wRec$summary.fixed ## Adult survival---- # Model formula_wSurv = wSurvival ~ DSstdsqrtExpScoreCor + I(DSstdsqrtExpScoreCor^2) + f(py, model = "iid2d", n=N.PY) + f(py2, DSstdsqrtExpScoreCor, copy = "py") + f(p, model = "iid2d", n=N.Plot) + f(p2, DSstdsqrtExpScoreCor, copy = "p") + f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, DSstdsqrtExpScoreCor, copy = "popyr") + f(pop, model = "iid2d", n=N.Pop) + f(pop2, DSstdsqrtExpScoreCor, copy = "pop") + f(yr, model = "iid2d", n=N.Year) + f(yr2, DSstdsqrtExpScoreCor, copy = "yr") modINLA_wSurv = inla (formula_wSurv, data = datPopComb, family = "gaussian", control.compute=list(cpo = TRUE, config = TRUE), control.predictor=list(compute=T)) summary(modINLA_wSurv) # Linear and nonlinear selection gradient fixefTable = modINLA_wSurv$summary.fixed ### Table 2 and Table S1---- ### Unstandardized selection gradients and proportion of variance ## Integrative fitness---- # Model formula_IndFit_DSstd = AnnualFitness ~ DSstdsqrtExpScoreCor + I(DSstdsqrtExpScoreCor^2) + f(pop, model = "iid2d", n=N.Pop) + f(pop2, DSstdsqrtExpScoreCor, copy = "pop") + f(p, model = "iid2d", n=N.Plot) + f(p2, DSstdsqrtExpScoreCor, copy = "p") + f(yr, model = "iid2d", n=N.Year) + f(yr2, DSstdsqrtExpScoreCor, copy = "yr") + f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, DSstdsqrtExpScoreCor, copy = "popyr") + f(py, model = "iid2d", n=N.PY) + f(py2, DSstdsqrtExpScoreCor, copy = "py") + modINLA_IndFit_DSstd = inla (formula_IndFit_DSstd, data = datPopComb, family = "gaussian", control.compute=list(cpo = TRUE, config = TRUE), control.predictor=list(compute=T)) summary(modINLA_IndFit_DSstd) # Table S1:Linear and nonlinear unstandardized selection gradient fixefTable = modINLA_IndFit_DSstd$summary.fixed # Table S1: Variance (intercept and slope) and intercept-slope correlation estimates # For macro-spatial variation (among populations) tau.int.Pop= modINLA_IndFit_DSstd$internal.marginals.hyperpar$"Log precision for pop (component 1)" s.marginal.pop <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Pop) ranef.int.Pop = as.data.frame (inla.zmarginal(s.marginal.pop)) ranefPopmean = ranef.int.Pop$mean tau.slope.Pop= modINLA_IndFit_DSstd$internal.marginals.hyperpar$"Log precision for pop (component 2)" s.marginal.slpop <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Pop) ranef.slope.Pop = as.data.frame (inla.zmarginal(s.marginal.slpop)) ranefslPopmean = ranef.slope.Pop$mean rho.int.Pop= modINLA_IndFit_DSstd$internal.marginals.hyperpar$"Rho_internal1:2 for pop" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Pop) ranef.rho.int.Pop = as.data.frame (inla.zmarginal(s.marginal)) # For micro-spatial variation (among plots) tau.int.Plot= modINLA_IndFit_DSstd$internal.marginals.hyperpar$"Log precision for p (component 1)" s.marginal.p <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Plot) ranef.int.Plot = as.data.frame (inla.zmarginal(s.marginal.p)) ranefPlotmean = ranef.int.Plot$mean tau.slope.Plot= modINLA_IndFit_DSstd$internal.marginals.hyperpar$"Log precision for p (component 2)" s.marginal.slp <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Plot) ranef.slope.Plot = as.data.frame (inla.zmarginal(s.marginal.slp)) ranefslPlotmean = ranef.slope.Plot$mean rho.int.Plot= modINLA_IndFit_DSstd$internal.marginals.hyperpar$"Rho_internal1:2 for p" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Plot) ranef.rho.int.Plot = as.data.frame (inla.zmarginal(s.marginal)) # For temporal variation (among years) tau.int.Yr = modINLA_IndFit_DSstd$internal.marginals.hyperpar$"Log precision for yr (component 1)" s.marginal.yr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Yr) ranef.int.Yr = as.data.frame (inla.zmarginal(s.marginal.yr)) ranefYrmean = ranef.int.Yr$mean tau.slope.Yr= modINLA_IndFit_DSstd$internal.marginals.hyperpar$"Log precision for yr (component 2)" s.marginal.slyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Yr) ranef.slope.Yr = as.data.frame (inla.zmarginal(s.marginal.slyr)) ranefslYrmean = ranef.slope.Yr$mean rho.int.Yr= modINLA_IndFit_DSstd$internal.marginals.hyperpar$"Rho_internal1:2 for yr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Yr) ranef.rho.int.Yr = as.data.frame (inla.zmarginal(s.marginal)) # For macro-scale temporal variation (among population x year) tau.int.PopYr= modINLA_IndFit_DSstd$internal.marginals.hyperpar$"Log precision for popyr (component 1)" s.marginal.popyr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PopYr) ranef.int.PopYr = as.data.frame (inla.zmarginal(s.marginal.popyr)) ranefPopYrmean = ranef.int.PopYr$mean tau.slope.PopYr= modINLA_IndFit_DSstd$internal.marginals.hyperpar$"Log precision for popyr (component 2)" s.marginal.slpopyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PopYr) ranef.slope.PopYr = as.data.frame (inla.zmarginal(s.marginal.slpopyr)) ranefslPopYrmean = ranef.slope.PopYr$mean rho.int.PopYr= modINLA_IndFit_DSstd$internal.marginals.hyperpar$"Rho_internal1:2 for popyr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PopYr) ranef.rho.int.PopYr = as.data.frame (inla.zmarginal(s.marginal)) # For micro-scale temporal variation (among plot x year) tau.int.PY= modINLA_IndFit_DSstd$internal.marginals.hyperpar$"Log precision for py (component 1)" s.marginal.py <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PY) ranef.int.PY = as.data.frame (inla.zmarginal(s.marginal.py)) ranefPYmean = ranef.int.PY$mean tau.slope.PY= modINLA_IndFit_DSstd$internal.marginals.hyperpar$"Log precision for py (component 2)" s.marginal.slpy <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PY) ranef.slope.PY = as.data.frame (inla.zmarginal(s.marginal.slpy)) ranefslPYmean = ranef.slope.PY$mean rho.int.PY= modINLA_IndFit_DSstd$internal.marginals.hyperpar$"Rho_internal1:2 for py" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PY) ranef.rho.int.PY = as.data.frame (inla.zmarginal(s.marginal)) # Residuals tau.resi= modINLA_IndFit_DSstd$internal.marginals.hyperpar$"Log precision for the Gaussian observations" s.marginal.resi <- inla.tmarginal(function(x) (1/exp(x)), tau.resi) ranef.resi = as.data.frame (inla.zmarginal(s.marginal.resi)) # Table 2: Proportion of variance attributable to each hierarchical level (i.e., random factor) # total variance in slopes totvar = s.marginal.slpop + s.marginal.slp + s.marginal.slyr + s.marginal.slpopyr + s.marginal.slpy # Macro-spatial variation: proportion of variance in slope explained by population r.slpop = s.marginal.slpop/totvar r.slpop.mean = mean(r.slpop[,1]) ci.pop = as.data.frame(quantile(r.slpop[,1], prob = c(0.025, 0.975))) r.slpop.df = paste0(round(r.slpop.mean,2)," ", "(",round(ci.pop[1,],2), ",", round(ci.pop[2,],2),")") # Micro-spatial variation: proportion of variance in slope explained by plot r.slp = s.marginal.slp/totvar r.slp.mean = mean(r.slp[,1]) ci.p = as.data.frame(quantile(r.slp[,1], prob = c(0.025, 0.975))) r.slp.df = paste0(round(r.slp.mean,2)," ", "(",round(ci.p[1,],2), ",", round(ci.p[2,],2),")") # Temporal variation: proportion of variance in slope explained by year r.slyr = s.marginal.slyr/totvar r.slyr.mean = mean(r.slyr[,1]) ci.yr = as.data.frame(quantile(r.slyr[,1], prob = c(0.025, 0.975))) r.slyr.df = paste0(round(r.slyr.mean,2)," ", "(",round(ci.yr[1,],2), ",", round(ci.yr[2,],2),")") # Macro-scale temporal variation: proportion of variance in slope explained by population-year r.slpopyr = s.marginal.slpopyr/totvar r.slpopyr.mean = mean(r.slpopyr[,1]) ci.popyr = as.data.frame(quantile(r.slpopyr[,1], prob = c(0.025, 0.975))) r.slpopyr.df = paste0(round(r.slpopyr.mean,2)," ", "(",round(ci.popyr[1,],2), ",", round(ci.popyr[2,],2),")") # Micro-scale temporal variation: proportion of variance in slope explained by plot-year r.slpy = s.marginal.slpy/totvar r.slpy.mean = mean(r.slpy[,1]) ci.py = as.data.frame(quantile(r.slpy[,1], prob = c(0.025, 0.975))) r.slpy.df = paste0(round(r.slpy.mean,2)," ", "(",round(ci.py[1,],2), ",", round(ci.py[2,],2),")") # Gather all proportions together into one data frame r.all.df = as.data.frame(rbind(r.slpop.df, r.slp.df, r.slyr.df, r.slpopyr.df, r.slpy.df)) ## Offspring recruitment---- # Define alpha and beta parameters for poisson model apar <- 0.5 bpar = apar*var(datPopComb$AnnualRecruits) lgprior <- list(prec = list(prior="loggamma", param = c(apar,bpar))) # Model formula_rndquad_wodens_DSstd = AnnualRecruits ~ DSstdsqrtExpScoreCor + I(DSstdsqrtExpScoreCor^2) + f(pop, model = "iid2d", n=N.Pop) + f(pop2, DSstdsqrtExpScoreCor, copy = "pop") + f(p, model = "iid2d", n=N.Plot) + f(p2, DSstdsqrtExpScoreCor, copy = "p") + f(yr, model = "iid2d", n=N.Year) + f(yr2, DSstdsqrtExpScoreCor, copy = "yr") + f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, DSstdsqrtExpScoreCor, copy = "popyr") + f(py, model = "iid2d", n = N.PY) + f(py2, DSstdsqrtExpScoreCor, copy = "py") + f(obsid, model = "iid", hyper = lgprior) modINLArndquad_wodens_DSstd = inla (formula_rndquad_wodens_DSstd, data = datPopComb, family = "poisson", control.compute=list(cpo=TRUE), control.predictor=list(compute=T)) # Check if there are any failures of the model (data point for which estimation poorly performed) sum(modINLArndquad_wodens_DSstd$cpo$failure) # recompute model failures manually modINLArndquad_wodens_DSstd = inla.cpo(modINLArndquad_wodens_DSstd) summary(modINLArndquad_wodens_DSstd) # Table S1 (Recruitment):Linear and nonlinear unstandardized selection gradient fixefTable = modINLArndquad_wodens_DSstd$summary.fixed # Table S1 (Recruitment): Variance (intercept and slope) and intercept-slope correlation estimates # For macro-spatial variation: among populations tau.int.Pop= modINLArndquad_wodens_DSstd$internal.marginals.hyperpar$"Log precision for pop (component 1)" s.marginal.pop <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Pop) ranef.int.Pop = as.data.frame (inla.zmarginal(s.marginal.pop)) ranefPopmean = ranef.int.Pop$mean tau.slope.Pop= modINLArndquad_wodens_DSstd$internal.marginals.hyperpar$"Log precision for pop (component 2)" s.marginal.slpop <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Pop) ranef.slope.Pop = as.data.frame (inla.zmarginal(s.marginal.slpop)) ranefslPopmean = ranef.slope.Pop$mean rho.int.Pop= modINLArndquad_wodens_DSstd$internal.marginals.hyperpar$"Rho_internal1:2 for pop" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Pop) ranef.rho.int.Pop = as.data.frame (inla.zmarginal(s.marginal)) # For micro-spatial variation: among plots tau.int.Plot= modINLArndquad_wodens_DSstd$internal.marginals.hyperpar$"Log precision for p (component 1)" s.marginal.p <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Plot) ranef.int.Plot = as.data.frame (inla.zmarginal(s.marginal.p)) ranefPlotmean = ranef.int.Plot$mean tau.slope.Plot= modINLArndquad_wodens_DSstd$internal.marginals.hyperpar$"Log precision for p (component 2)" s.marginal.slp <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Plot) ranef.slope.Plot = as.data.frame (inla.zmarginal(s.marginal.slp)) ranefslPlotmean = ranef.slope.Plot$mean rho.int.Plot= modINLArndquad_wodens_DSstd$internal.marginals.hyperpar$"Rho_internal1:2 for p" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Plot) ranef.rho.int.Plot = as.data.frame (inla.zmarginal(s.marginal)) # For temporal variation: among years tau.int.Yr= modINLArndquad_wodens_DSstd$internal.marginals.hyperpar$"Log precision for yr (component 1)" s.marginal.yr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Yr) ranef.int.Yr = as.data.frame (inla.zmarginal(s.marginal.yr)) ranefYrmean = ranef.int.Yr$mean tau.slope.Yr= modINLArndquad_wodens_DSstd$internal.marginals.hyperpar$"Log precision for yr (component 2)" s.marginal.slyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Yr) ranef.slope.Yr = as.data.frame (inla.zmarginal(s.marginal.slyr)) ranefslYrmean = ranef.slope.Yr$mean rho.int.Yr= modINLArndquad_wodens_DSstd$internal.marginals.hyperpar$"Rho_internal1:2 for yr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Yr) ranef.rho.int.Yr = as.data.frame (inla.zmarginal(s.marginal)) # For macro-scale temporal variation: among population-years tau.int.PopYr= modINLArndquad_wodens_DSstd$internal.marginals.hyperpar$"Log precision for popyr (component 1)" s.marginal.popyr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PopYr) ranef.int.PopYr = as.data.frame (inla.zmarginal(s.marginal.popyr)) ranefPopYrmean = ranef.int.PopYr$mean tau.slope.PopYr= modINLArndquad_wodens_DSstd$internal.marginals.hyperpar$"Log precision for popyr (component 2)" s.marginal.slpopyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PopYr) ranef.slope.PopYr = as.data.frame (inla.zmarginal(s.marginal.slpopyr)) ranefslPopYrmean = ranef.slope.PopYr$mean rho.int.PopYr= modINLArndquad_wodens_DSstd$internal.marginals.hyperpar$"Rho_internal1:2 for popyr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PopYr) ranef.rho.int.PopYr = as.data.frame (inla.zmarginal(s.marginal)) # For micro-scale temporal variation: among plot-years tau.int.PY= modINLArndquad_wodens_DSstd$internal.marginals.hyperpar$"Log precision for py (component 1)" s.marginal.py <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PY) ranef.int.PY = as.data.frame (inla.zmarginal(s.marginal.py)) ranefPYmean = ranef.int.PY$mean tau.slope.PY= modINLArndquad_wodens_DSstd$internal.marginals.hyperpar$"Log precision for py (component 2)" s.marginal.slpy <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PY) ranef.slope.PY = as.data.frame (inla.zmarginal(s.marginal.slpy)) ranefslPYmean = ranef.slope.PY$mean rho.int.PY= modINLArndquad_wodens_DSstd$internal.marginals.hyperpar$"Rho_internal1:2 for py" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PY) ranef.rho.int.PY = as.data.frame (inla.zmarginal(s.marginal)) # Overdispersion estimate tau.int.obsid= modINLArndquad_wodens_DSstd$internal.marginals.hyperpar$"Log precision for obsid" s.marginal.obsid <- inla.tmarginal(function(x) (1/exp(x)), tau.int.obsid) ranef.int.obsid = as.data.frame (inla.zmarginal(s.marginal.obsid)) # Table S1 (Recruitment) - Permutation tests (permuting labels of each level at which selection can vary) # to test whether slopes estimated at each level are significantly different from each other # Macro-spatial variation (among populations) significance: Population label permutation sformula_dsstd_Recruitwodens_popRnd = AnnualRecruits ~ DSstdsqrtExpScoreCor + I(DSstdsqrtExpScoreCor^2) + f(spop, model = "iid2d", n=sN.Pop) + f(spop2, DSstdsqrtExpScoreCor, copy = "spop") + f(p, model = "iid2d", n=N.Plot) + f(p2, DSstdsqrtExpScoreCor, copy = "p") + f(yr, model = "iid2d", n=N.Year) + f(yr2, DSstdsqrtExpScoreCor, copy = "yr") + f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, DSstdsqrtExpScoreCor, copy = "popyr") + f(py, model = "iid2d", n = N.PY) + f(py2, DSstdsqrtExpScoreCor, copy = "py") + f(obsid, model = "iid", hyper = lgprior) n.sim <- 1000 respop_dsstd_recruit<-matrix(NA, n.sim,2) for(i in 1:n.sim){ datPopComb$sPopulationID<-sample(datPopComb$PopulationID,nrow(datPopComb))## randomize plot labels spop = as.numeric(datPopComb$sPopulationID) sn.pop = max(spop) sN.Pop = 2*sn.pop spop2 = as.numeric(datPopComb$sPopulationID) + sn.pop datPopComb$spop = spop datPopComb$spop2 = spop2 modRecruit_dsstd_popS.Inla = inla (sformula_dsstd_Recruitwodens_popRnd, data = datPopComb, family = "poisson") tau.int.Pop= modRecruit_dsstd_popS.Inla$internal.marginals.hyperpar$"Log precision for spop (component 1)" s.marginal <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Pop) ranef1 = as.data.frame (inla.zmarginal(s.marginal)) respop_dsstd_recruit[i,1]<-ranef1$mean tau.slope.Pop= modRecruit_dsstd_popS.Inla$internal.marginals.hyperpar$"Log precision for spop (component 2)" s.marginal <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Pop) ranef2 = as.data.frame (inla.zmarginal(s.marginal)) respop_dsstd_recruit[i,2]<-ranef2$mean } # Permutation p-value: Amount of variation expected by chance larger than observed variation in selection table(respop_dsstd_recruit[,1]>ranefPopmean)/n.sim table(respop_dsstd_recruit[,2]>ranefslPopmean)/n.sim # Save model output saveRDS(respop_dsstd_recruit[,1], file = "PopPermutationIntercept_DSstd_Recruit.rds") saveRDS(respop_dsstd_recruit[,2], file = "PopPermutationSlope_DSstd_Recruit.rds") # Micro-spatial variation (among plots) significance: Plot label permutation sformula_dsstd_Recruitwodens_plotRnd = AnnualRecruits ~ DSstdsqrtExpScoreCor + I(DSstdsqrtExpScoreCor^2) + f(pop, model = "iid2d", n=N.Pop) + f(pop2, DSstdsqrtExpScoreCor, copy = "pop") + f(sp, model = "iid2d", n=sN.Plot) + f(sp2, DSstdsqrtExpScoreCor, copy = "sp") + f(yr, model = "iid2d", n=N.Year) + f(yr2, DSstdsqrtExpScoreCor, copy = "yr") + f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, DSstdsqrtExpScoreCor, copy = "popyr") + f(py, model = "iid2d", n = N.PY) + f(py2, DSstdsqrtExpScoreCor, copy = "py") + f(obsid, model = "iid", hyper = lgprior) n.sim <- 1000 resplot_dsstd_recruit<-matrix(NA, n.sim,2) for(i in 1:n.sim){ datPopComb$sPlot<-sample(datPopComb$Plot,nrow(datPopComb))## randomize plot labels sn.plot = nlevels(datPopComb$sPlot) sN.Plot = 2*sn.plot sp = as.numeric(datPopComb$sPlot) sp2 = as.numeric(datPopComb$sPlot) + sn.plot datPopComb$sp = sp datPopComb$sp2 = sp2 modRecruit_dsstd_plotS.Inla = inla (sformula_dsstd_Recruitwodens_plotRnd, data = datPopComb, family = "poisson") tau.int.Plot= modRecruit_dsstd_plotS.Inla$internal.marginals.hyperpar$"Log precision for sp (component 1)" s.marginal <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Plot) ranef1 = as.data.frame (inla.zmarginal(s.marginal)) resplot_dsstd_recruit[i,1]<-ranef1$mean tau.slope.Plot= modRecruit_dsstd_plotS.Inla$internal.marginals.hyperpar$"Log precision for sp (component 2)" s.marginal.sl <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Plot) ranef2 = as.data.frame (inla.zmarginal(s.marginal.sl)) resplot_dsstd_recruit[i,2]<-ranef2$mean } # Permutation p-value: Amount of variation expected by chance larger than observed variation in selection table(resplot_dsstd_recruit[,1]>ranefPlotmean)/n.sim table(resplot_dsstd_recruit[,2]>ranefslPlotmean)/n.sim # Save model output saveRDS(resplot_dsstd_recruit[,1], file = "PlotPermutationIntercept_DSstd_Recruit.rds") saveRDS(resplot_dsstd_recruit[,2], file = "PlotPermutationSlope_DSstd_Recruit.rds") # Temporal variation (among years) significance: Year label permutation sformula_dsstd_Recruitwodens_yrRnd = AnnualRecruits ~ PYstdsqrtExpScoreCor + I(PYstdsqrtExpScoreCor^2) + f(pop, model = "iid2d", n=N.Pop) + f(pop2, PYstdsqrtExpScoreCor, copy = "pop") + f(p, model = "iid2d", n=N.Plot) + f(p2, PYstdsqrtExpScoreCor, copy = "p") + f(syr, model = "iid2d", n=sN.Year) + f(syr2, PYstdsqrtExpScoreCor, copy = "syr") + f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, PYstdsqrtExpScoreCor, copy = "popyr") + f(py, model = "iid2d", n = N.PY) + f(py2, PYstdsqrtExpScoreCor, copy = "py") + f(obsid, model = "iid", hyper = lgprior) n.sim <- 1000 resyr_dsstd_recruit<-matrix(NA, n.sim,2) for(i in 1:n.sim){ datPopComb$sBroodYear<-sample(datPopComb$BroodYear,nrow(datPopComb))## randomize plot labels syr = as.numeric(datPopComb$sBroodYear) sn.year = max(syr) sN.Year = 2*sn.year syr2 = as.numeric(datPopComb$sBroodYear) + sn.year datPopComb$syr = syr datPopComb$syr2 = syr2 modRecruit_dsstd_yrS.Inla = inla (sformula_dsstd_Recruitwodens_yrRnd, data = datPopComb2, family = "poisson", control.inla = list(strategy = "laplace", npoints = 21)) tau.int.Yr= modRecruit_dsstd_yrS.Inla$internal.marginals.hyperpar$"Log precision for syr (component 1)" s.marginal <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Yr) ranef1 = as.data.frame (inla.zmarginal(s.marginal)) resyr_dsstd_recruit[i,1]<-ranef1$mean tau.slope.Yr= modRecruit_dsstd_yrS.Inla$internal.marginals.hyperpar$"Log precision for syr (component 2)" s.marginal <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Yr) ranef2 = as.data.frame (inla.zmarginal(s.marginal)) resyr_dsstd_recruit[i,2]<-ranef2$mean } hist(resyr_dsstd_recruit[,1], xlim = c(min(resyr_dsstd_recruit), ranefYrmean), xlab = "Variation in intercept", main = "", las = 1) abline(v=ranefYrmean, lwd=4) hist(resyr_dsstd_recruit[,1], xlim = c(min(resyr_dsstd_recruit[,1]), ranefYrmean), xlab = "Variation in intercept", main = "", las = 1) abline(v=ranefYrmean, lwd=4, col="orange") hist(resyr_dsstd_recruit[,2], xlim = c(min(resyr_dsstd_recruit[,2]), ranefslYrmean), xlab = "Variation in selection", main = "", las = 1) abline(v=ranefslYrmean, lwd=4) hist(resyr_dsstd_recruit[,2],xlim = c(min(resyr_dsstd_recruit[,2]), ranefslYrmean), xlab = "Variation in selection", main = "", las = 1) abline(v=ranefslYrmean, lwd=4, col = "orange") # Permutation p-value: Amount of variation expected by chance larger than observed variation in selection table(resyr_dsstd_recruit[,1]>ranefYrmean)/n.sim table(resyr_dsstd_recruit[,2]>ranefslYrmean)/n.sim # Save model output saveRDS(resyr_dsstd_recruit[,1], file = "YearPermutationIntercept_DSstd_Recruit.rds") saveRDS(resyr_dsstd_recruit[,2], file = "YearPermutationSlope_DSstd_Recruit.rds") # Macro-scale temporal variation (among population-years) significance: Population-Year label permutation. sformula_dsstd_Recruit_popyrRnd = AnnualRecruits ~ DSstdsqrtExpScoreCor + I(DSstdsqrtExpScoreCor^2) + f(pop, model = "iid2d", n=N.Pop) + f(pop2, DSstdsqrtExpScoreCor, copy = "pop") + f(p, model = "iid2d", n=N.Plot) + f(p2, DSstdsqrtExpScoreCor, copy = "p") + f(yr, model = "iid2d", n=N.Year) + f(yr2, DSstdsqrtExpScoreCor, copy = "yr") + f(spopyr, model = "iid2d", n=sN.PopYear) + f(spopyr2, DSstdsqrtExpScoreCor, copy = "spopyr") + f(py, model = "iid2d", n = N.PY) + f(py2, DSstdsqrtExpScoreCor, copy = "py") + f(obsid, model = "iid", hyper = lgprior) n.sim <- 1000 respopyr_dsstd_recruit<-matrix(NA, n.sim,2) for(i in 1:n.sim){ datPopComb$sPopYear<-sample(datPopComb$PopYear,nrow(datPopComb))## randomize year labels sn.popyear = nlevels(datPopComb$sPopYear) sN.PopYear = 2*sn.popyear spopyr = as.numeric(datPopComb$sPopYear) spopyr2 = as.numeric(datPopComb$sPopYear) + sn.popyear datPopComb$spopyr = spopyr datPopComb$spopyr2 = spopyr2 modRecruit_dsstd_popyrS.Inla = inla (sformula_dsstd_Recruit_popyrRnd, data = datPopComb, family = "poisson") tau.int.PopYear= modRecruit_dsstd_popyrS.Inla$internal.marginals.hyperpar$"Log precision for spopyr (component 1)" s.marginal <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PopYear) ranef1 = as.data.frame (inla.zmarginal(s.marginal)) respopyr_dsstd_recruit[i,1]<-ranef1$mean tau.slope.PopYear= modRecruit_dsstd_popyrS.Inla$internal.marginals.hyperpar$"Log precision for spopyr (component 2)" s.marginal <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PopYear) ranef2 = as.data.frame (inla.zmarginal(s.marginal)) respopyr_dsstd_recruit[i,2]<-ranef2$mean } # Permutation p-value: Amount of variation expected by chance larger than observed variation in selection table(respopyr_dsstd_recruit[,1]>ranefYrmean)/n.sim table(respopyr_dsstd_recruit[,2]>ranefslYrmean)/n.sim # Save model output saveRDS(respopyr_dsstd_recruit[,1], file = "PopYearPermutationIntercept_DSstd_Recruit.rds") saveRDS(respopyr_dsstd_recruit[,2], file = "PopYearPermutationSlope_DSstd_Recruit.rds") # Micro-scale temporal variation (among plot-years) significance: Plot-Year label permutation sformula_dsstd_Recruitwodens_pyRnd = AnnualRecruits ~ DSstdsqrtExpScoreCor + I(DSstdsqrtExpScoreCor^2) + f(pop, model = "iid2d", n=N.Pop) + f(pop2, DSstdsqrtExpScoreCor, copy = "pop") + f(p, model = "iid2d", n=N.Plot) + f(p2, DSstdsqrtExpScoreCor, copy = "p") + f(yr, model = "iid2d", n=N.Year) + f(yr2, DSstdsqrtExpScoreCor, copy = "yr") + f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, DSstdsqrtExpScoreCor, copy = "popyr") + f(spy, model = "iid2d", n = sN.PY) + f(spy2, DSstdsqrtExpScoreCor, copy = "spy") + f(obsid, model = "iid", hyper = lgprior) n.sim <- 1000 respy_dsstd_recruit<-matrix(NA, n.sim,2) for(i in 1:n.sim){ datPopComb$sPlotYear<-sample(datPopComb$PlotYear,nrow(datPopComb))## randomize plot labels sn.PY = nlevels(datPopComb$sPlotYear) sN.PY = 2*sn.PY spy = as.numeric(datPopComb$sPlotYear) spy2 = as.numeric(datPopComb$sPlotYear) + sn.PY datPopComb$spy = spy datPopComb$spy2 = spy2 modRecruit_dsstd_pyS.Inla = inla (sformula_dsstd_Recruitwodens_pyRnd, data = datPopComb, family = "poisson") tau.int.PY= modRecruit_dsstd_pyS.Inla$internal.marginals.hyperpar$"Log precision for spy (component 1)" s.marginal <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PY) ranef1 = as.data.frame (inla.zmarginal(s.marginal)) respy_dsstd_recruit[i,1]<-ranef1$mean tau.slope.PY= modRecruit_dsstd_pyS.Inla$internal.marginals.hyperpar$"Log precision for spy (component 2)" s.marginal <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PY) ranef2 = as.data.frame (inla.zmarginal(s.marginal)) respy_dsstd_recruit[i,2]<-ranef2$mean } # Permutation p-value: Amount of variation expected by chance larger than observed variation in selection table(respy_dsstd_recruit[,1]>ranefPYmean)/n.sim table(respy_dsstd_recruit[,2]>ranefslPYmean)/n.sim # Save model output saveRDS(respy_dsstd_recruit[,1], file = "PlotYearPermutationIntercept_DSstd_Recruit.rds") saveRDS(respy_dsstd_recruit[,2], file = "PlotYearPermutationSlope_DSstd_Recruit.rds") # Table 2 (Recruitment): Proportion of variance attributable to each hierarchical level (i.e., random factor) # total variance slopes totvar = s.marginal.slpy + s.marginal.slp + s.marginal.slpopyr + s.marginal.slpop + s.marginal.slyr # Macro-spatial variation: proportion of variance in slope explained by population r.slpop = s.marginal.slpop/totvar r.slpop.mean = mean(r.slpop[,1]) ci.pop = as.data.frame(quantile(r.slpop[,1], prob = c(0.025, 0.975))) r.slpop.df = paste0(round(r.slpop.mean,2)," ", "(",round(ci.pop[1,],2), ",", round(ci.pop[2,],2),")") # Micro-spatial variation: proportion of variance in slope explained by plot r.slp = s.marginal.slp/totvar r.slp.mean = mean(r.slp[,1]) ci.p = as.data.frame(quantile(r.slp[,1], prob = c(0.025, 0.975))) r.slp.df = paste0(round(r.slp.mean,2)," ", "(",round(ci.p[1,],2), ",", round(ci.p[2,],2),")") # Temporal variation: proportion of variance in slope explained by year r.slyr = s.marginal.slyr/totvar r.slyr.mean = mean(r.slyr[,1]) ci.yr = as.data.frame(quantile(r.slyr[,1], prob = c(0.025, 0.975))) r.slyr.df = paste0(round(r.slyr.mean,2)," ", "(",round(ci.yr[1,],2), ",", round(ci.yr[2,],2),")") # Macro-scale temporal variation: proportion of variance in slope explained by population-year r.slpopyr = s.marginal.slpopyr/totvar r.slpopyr.mean = mean(r.slpopyr[,1]) ci.popyr = as.data.frame(quantile(r.slpopyr[,1], prob = c(0.025, 0.975))) r.slpopyr.df = paste0(round(r.slpopyr.mean,2)," ", "(",round(ci.popyr[1,],2), ",", round(ci.popyr[2,],2),")") # Micro-scale temporal variation: proportion of variance in slope explained by plot-year r.slpy = s.marginal.slpy/totvar r.slpy.mean = mean(r.slpy[,1]) ci.py = as.data.frame(quantile(r.slpy[,1], prob = c(0.025, 0.975))) r.slpy.df = paste0(round(r.slpy.mean,2)," ", "(",round(ci.py[1,],2), ",", round(ci.py[2,],2),")") # Gather all proportions together into one dataframe r.all.df = as.data.frame(rbind(r.slpop.df, r.slp.df, r.slyr.df, r.slpopyr.df, r.slpy.df)) # Adult survival---- # Model formula_Survwodens_DSstdExp = Survived ~ DSstdsqrtExpScoreCor + I(DSstdsqrtExpScoreCor^2) + f(pop, model = "iid2d", n=N.Pop) + f(pop2, DSstdsqrtExpScoreCor, copy = "pop") + f(p, model = "iid2d", n=N.Plot) + f(p2, DSstdsqrtExpScoreCor, copy = "p") + f(yr, model = "iid2d", n=N.Year) + f(yr2, DSstdsqrtExpScoreCor, copy = "yr") + f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, DSstdsqrtExpScoreCor, copy = "popyr") + f(py, model = "iid2d", n=N.PY) + f(py2, DSstdsqrtExpScoreCor, copy = "py") modINLA_Survwodens_DSstdExp = inla (formula_Survwodens_DSstdExp, data = datPopComb, family = "binomial", control.compute=list(cpo = TRUE, config = TRUE), control.predictor=list(compute=T)) # Check if there are any model failures sum(modINLA_Survwodens_DSstdExp$cpo$failure) # recompute model failures manually modINLA_Survwodens_DSstdExp = inla.cpo(modINLA_Survwodens_DSstdExp) summary(modINLA_Survwodens_DSstdExp) # Table S1 (Survival): Linear and nonlinear unstandardized selection gradient fixefTable = modINLA_Survwodens_DSstdExp$summary.fixed # Table S1 (Survival): Variance (intercept and slope) and intercept-slope correlation estimates # For macro-spatial variation: among populations tau.int.Pop= modINLA_Survwodens_DSstdExp$internal.marginals.hyperpar$"Log precision for pop (component 1)" s.marginal.pop <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Pop) ranef.int.Pop = as.data.frame (inla.zmarginal(s.marginal.pop)) ranefPopmean = ranef.int.Pop$mean tau.slope.Pop= modINLA_Survwodens_DSstdExp$internal.marginals.hyperpar$"Log precision for pop (component 2)" s.marginal.slpop <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Pop) ranef.slope.Pop = as.data.frame (inla.zmarginal(s.marginal.slpop)) ranefslPopmean = ranef.slope.Pop$mean rho.int.Pop= modINLA_Survwodens_DSstdExp$internal.marginals.hyperpar$"Rho_internal1:2 for pop" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Pop) ranef.rho.int.Pop = as.data.frame (inla.zmarginal(s.marginal)) # For micro-spatial variation: among plots tau.int.Plot= modINLA_Survwodens_DSstdExp$internal.marginals.hyperpar$"Log precision for p (component 1)" s.marginal.p <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Plot) ranef.int.Plot = as.data.frame (inla.zmarginal(s.marginal.p)) ranefPlotmean = ranef.int.Plot$mean tau.slope.Plot= modINLA_Survwodens_DSstdExp$internal.marginals.hyperpar$"Log precision for p (component 2)" s.marginal.slp <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Plot) ranef.slope.Plot = as.data.frame (inla.zmarginal(s.marginal.slp)) ranefslPlotmean = ranef.slope.Plot$mean rho.int.Plot= modINLA_Survwodens_DSstdExp$internal.marginals.hyperpar$"Rho_internal1:2 for p" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Plot) ranef.rho.int.Plot = as.data.frame (inla.zmarginal(s.marginal)) # For temporal variation: among years tau.int.Yr = modINLA_Survwodens_DSstdExp$internal.marginals.hyperpar$"Log precision for yr (component 1)" s.marginal.yr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Yr) ranef.int.Yr = as.data.frame (inla.zmarginal(s.marginal.yr)) ranefYrmean = ranef.int.Yr$mean tau.slope.Yr= modINLA_Survwodens_DSstdExp$internal.marginals.hyperpar$"Log precision for yr (component 2)" s.marginal.slyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Yr) ranef.slope.Yr = as.data.frame (inla.zmarginal(s.marginal.slyr)) ranefslYrmean = ranef.slope.Yr$mean rho.int.Yr= modINLA_Survwodens_DSstdExp$internal.marginals.hyperpar$"Rho_internal1:2 for yr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Yr) ranef.rho.int.Yr = as.data.frame (inla.zmarginal(s.marginal)) # For macro-scale temporal variation: among population-years tau.int.PopYr= modINLA_Survwodens_DSstdExp$internal.marginals.hyperpar$"Log precision for popyr (component 1)" s.marginal.popyr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PopYr) ranef.int.PopYr = as.data.frame (inla.zmarginal(s.marginal.popyr)) ranefPopYrmean = ranef.int.PopYr$mean tau.slope.PopYr= modINLA_Survwodens_DSstdExp$internal.marginals.hyperpar$"Log precision for popyr (component 2)" s.marginal.slpopyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PopYr) ranef.slope.PopYr = as.data.frame (inla.zmarginal(s.marginal.slpopyr)) ranefslPopYrmean = ranef.slope.PopYr$mean rho.int.PopYr= modINLA_Survwodens_DSstdExp$internal.marginals.hyperpar$"Rho_internal1:2 for popyr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PopYr) ranef.rho.int.PopYr = as.data.frame (inla.zmarginal(s.marginal)) # For micro-scale temporal variation: among plot-years tau.int.PY= modINLA_Survwodens_DSstdExp$internal.marginals.hyperpar$"Log precision for py (component 1)" s.marginal.py <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PY) ranef.int.PY = as.data.frame (inla.zmarginal(s.marginal.py)) ranefPYmean = ranef.int.PY$mean tau.slope.PY= modINLA_Survwodens_DSstdExp$internal.marginals.hyperpar$"Log precision for py (component 2)" s.marginal.slpy <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PY) ranef.slope.PY = as.data.frame (inla.zmarginal(s.marginal.slpy)) ranefslPYmean = ranef.slope.PY$mean rho.int.PY= modINLA_Survwodens_DSstdExp$internal.marginals.hyperpar$"Rho_internal1:2 for py" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PY) ranef.rho.int.PY = as.data.frame (inla.zmarginal(s.marginal)) # Table S1 (Survival) - Permutation tests (permuting labels of each level at which selection can vary) # to test whether slopes estimated at each level are significantly different from each other # Macro-spatial variation (among populations) significance: Population label permutation sformula_Survwodens_DSstdExp_popRnd = Survived ~ DSstdsqrtExpScoreCor + I(DSstdsqrtExpScoreCor^2) + f(spop, model = "iid2d", n=sN.Pop) + f(spop2, DSstdsqrtExpScoreCor, copy = "spop") + f(yr, model = "iid2d", n=N.Year) + f(yr2, DSstdsqrtExpScoreCor, copy = "yr")+ f(p, model = "iid2d", n=N.Plot) + f(p2, DSstdsqrtExpScoreCor, copy = "p") + f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, DSstdsqrtExpScoreCor, copy = "popyr") + f(py, model = "iid2d", n=N.PY) + f(py2, DSstdsqrtExpScoreCor, copy = "py") n.sim <- 1000 respop_surv_DSstdExp<-matrix(NA, n.sim,2) for(i in 1:n.sim){ datPopComb$sPopulationID<-sample(datPopComb$PopulationID,nrow(datPopComb))## randomize plot labels spop = as.numeric(datPopComb$sPopulationID) sn.pop = max(spop) sN.Pop = 2*sn.pop spop2 = as.numeric(datPopComb$sPopulationID) + sn.pop datPopComb$spop = spop datPopComb$spop2 = spop2 modSurv_DSstdExp_popS.Inla = inla (sformula_Survwodens_DSstdExp_popRnd, data = datPopComb, family = "binomial", control.predictor=list(compute=T)) tau.int.Pop.DSstdExp = modSurv_DSstdExp_popS.Inla$internal.marginals.hyperpar$"Log precision for spop (component 1)" s.marginal.DSstdExp <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Pop.DSstdExp) ranef1.DSstdExp = as.data.frame (inla.zmarginal(s.marginal.DSstdExp)) respop_surv_DSstdExp[i,1]<-ranef1.DSstdExp$mean tau.slope.Pop.DSstdExp= modSurv_DSstdExp_popS.Inla$internal.marginals.hyperpar$"Log precision for spop (component 2)" s.marginal.sl.DSstdExp <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Pop.DSstdExp) ranef2.DSstdExp = as.data.frame (inla.zmarginal(s.marginal.sl.DSstdExp)) respop_surv_DSstdExp[i,2]<-ranef2.DSstdExp$mean } # Permutation p-value: Amount of variation expected by chance larger than observed variation in selection table(respop_surv_DSstdExp[,1]>ranefPopmean)/n.sim table(respop_surv_DSstdExp[,2]>ranefslPopmean)/n.sim # Save model output saveRDS(respop_surv_DSstdExp[,1], file = "PopPermutation_DSstdExp_Intercept_Survival.rds") saveRDS(respop_surv_DSstdExp[,2], file = "PopPermutation_DSstdExp_Slope_Survival.rds") # Micro-spatial variation (among plots) significance: Plot label permutation sformula_Survwodens_DSstdExp_plotRnd = Survived ~ DSstdsqrtExpScoreCor + I(DSstdsqrtExpScoreCor^2) + f(pop, model = "iid2d", n=N.Pop) + f(pop2, DSstdsqrtExpScoreCor, copy = "pop") + f(sp, model = "iid2d", n=sN.Plot) + f(sp2, DSstdsqrtExpScoreCor, copy = "sp") + f(yr, model = "iid2d", n=N.Year) + f(yr2, DSstdsqrtExpScoreCor, copy = "yr") + f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, DSstdsqrtExpScoreCor, copy = "popyr") + f(py, model = "iid2d", n=N.PY) + f(py2, DSstdsqrtExpScoreCor, copy = "py") n.sim <- 1000 resplot_surv_DSstdExp<-matrix(NA, n.sim,2) for(i in 1:n.sim){ datPopComb$sPlot<-sample(datPopComb$Plot,nrow(datPopComb))## randomize plot labels sn.plot = nlevels(datPopComb$sPlot) sN.Plot = 2*sn.plot sp = as.numeric(datPopComb$sPlot) sp2 = as.numeric(datPopComb$sPlot) + sn.plot datPopComb$sp = sp datPopComb$sp2 = sp2 modSurv_DSstdExp_plotS.Inla = inla (sformula_Survwodens_DSstdExp_plotRnd, data = datPopComb, family = "binomial", control.predictor=list(compute=T)) tau.int.Plot.DSstdExp = modSurv_DSstdExp_plotS.Inla$internal.marginals.hyperpar$"Log precision for sp (component 1)" s.marginal.DSstdExp <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Plot.DSstdExp) ranef1.DSstdExp = as.data.frame (inla.zmarginal(s.marginal.DSstdExp)) resplot_surv_DSstdExp[i,1]<-ranef1.DSstdExp$mean tau.slope.Plot.DSstdExp= modSurv_DSstdExp_plotS.Inla$internal.marginals.hyperpar$"Log precision for sp (component 2)" s.marginal.sl.DSstdExp <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Plot.DSstdExp) ranef2.DSstdExp = as.data.frame (inla.zmarginal(s.marginal.sl.DSstdExp)) resplot_surv_DSstdExp[i,2]<-ranef2.DSstdExp$mean } # Permutation p-value: Amount of variation expected by chance larger than observed variation in selection table(resplot_surv_DSstdExp[,1]>ranefPlotmean)/n.sim table(resplot_surv_DSstdExp[,2]>ranefslPlotmean)/n.sim # Save model output saveRDS(resplot_surv_DSstdExp[,1], file = "PlotPermutation_DSstdExp_Intercept_Survival.rds") saveRDS(resplot_surv_DSstdExp[,2], file = "PlotPermutation_DSstdExp_Slope_Survival.rds") # Temporal variation (among years) significance: Year label permutation sformula_Survwodens_DSstdExp_yrRnd = Survived ~ DSstdsqrtExpScoreCor + I(DSstdsqrtExpScoreCor^2) + f(pop, model = "iid2d", n=N.Pop) + f(pop2, DSstdsqrtExpScoreCor, copy = "pop") + f(p, model = "iid2d", n=N.Plot) + f(p2, DSstdsqrtExpScoreCor, copy = "p") + f(syr, model = "iid2d", n=sN.Year) + f(syr2, DSstdsqrtExpScoreCor, copy = "syr") + f(popyr, model = "iid2d", n=N.PopYear) + f(spopyr2, DSstdsqrtExpScoreCor, copy = "popyr") + f(py, model = "iid2d", n=N.PY) + f(py2, DSstdsqrtExpScoreCor, copy = "py") n.sim <- 1000 resyr_surv_DSstdExp<-matrix(NA, n.sim,2) for(i in 1:n.sim){ datPopComb$sBroodYear<-sample(datPopComb$BroodYear,nrow(datPopComb))## randomize year labels sn.year = nlevels(datPopComb$sBroodYear) sN.Year = 2*sn.year syr = as.numeric(datPopComb$sBroodYear) syr2 = as.numeric(datPopComb$sBroodYear) + sn.year datPopComb$syr = syr datPopComb$syr2 = syr2 modSurv_DSstdExp_yrS.Inla = inla (sformula_Survwodens_DSstdExp_yrRnd, data = datPopComb, family = "binomial", control.predictor=list(compute=T)) tau.int.Yr.DSstdExp = modSurv_DSstdExp_yrS.Inla$internal.marginals.hyperpar$"Log precision for syr (component 1)" s.marginal.DSstdExp <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Yr.DSstdExp) ranef1.DSstdExp = as.data.frame (inla.zmarginal(s.marginal.DSstdExp)) resyr_surv_DSstdExp[i,1]<-ranef1.DSstdExp$mean tau.slope.Yr.DSstdExp = modSurv_DSstdExp_yrS.Inla$internal.marginals.hyperpar$"Log precision for syr (component 2)" s.marginal.sl.DSstdExp <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Yr.DSstdExp) ranef2.DSstdExp = as.data.frame (inla.zmarginal(s.marginal.sl.DSstdExp)) resyr_surv_DSstdExp[i,2]<-ranef2.DSstdExp$mean } # Permutation p-value: Amount of variation expected by chance larger than observed variation in selection table(resyr_surv_DSstdExp[,1]>ranefYrmean)/n.sim table(resyr_surv_DSstdExp[,2]>ranefslYrmean)/n.sim # Save model output saveRDS(resyr_surv_DSstdExp[,1], file = "AmgPopYearPermutation_DSstdExp_Intercept_Survival.rds") saveRDS(resyr_surv_DSstdExp[,2], file = "AmgPopYearPermutation_DSstdExp_Slope_Survival.rds") # Macro-scale temporal variation (among population-years) significance:Population-Year label permutation. sformula_Survwodens_DSstdExp_popyrRnd = Survived ~ DSstdsqrtExpScoreCor + I(DSstdsqrtExpScoreCor^2) + f(pop, model = "iid2d", n=N.Pop) + f(pop2, DSstdsqrtExpScoreCor, copy = "pop") + f(p, model = "iid2d", n=N.Plot) + f(p2, DSstdsqrtExpScoreCor, copy = "p") + f(yr, model = "iid2d", n=N.Year) + f(yr2, DSstdsqrtExpScoreCor, copy = "yr") + f(spopyr, model = "iid2d", n=sN.PopYear) + f(spopyr2, DSstdsqrtExpScoreCor, copy = "spopyr") + f(py, model = "iid2d", n=N.PY) + f(py2, DSstdsqrtExpScoreCor, copy = "py") n.sim <- 1000 respopyr_surv_DSstdExp<-matrix(NA, n.sim,2) for(i in 1:n.sim){ datPopComb$sPopYear<-sample(datPopComb$PopYear,nrow(datPopComb))## randomize year labels sn.popyear = nlevels(datPopComb$sPopYear) sN.PopYear = 2*sn.popyear spopyr = as.numeric(datPopComb$sPopYear) spopyr2 = as.numeric(datPopComb$sPopYear) + sn.popyear datPopComb$spopyr = spopyr datPopComb$spopyr2 = spopyr2 modSurv_DSstdExp_popyrS.Inla = inla (sformula_Survwodens_DSstdExp_popyrRnd, data = datPopComb, family = "binomial", control.predictor=list(compute=T)) tau.int.PopYear.DSstdExp = modSurv_DSstdExp_popyrS.Inla$internal.marginals.hyperpar$"Log precision for spopyr (component 1)" s.marginal.DSstdExp <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PopYear.DSstdExp) ranef1.DSstdExp = as.data.frame (inla.zmarginal(s.marginal.DSstdExp)) respopyr_surv_DSstdExp[i,1]<-ranef1.DSstdExp$mean tau.slope.PopYear.DSstdExp = modSurv_DSstdExp_popyrS.Inla$internal.marginals.hyperpar$"Log precision for spopyr (component 2)" s.marginal.sl.DSstdExp <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PopYear.DSstdExp) ranef2.DSstdExp = as.data.frame (inla.zmarginal(s.marginal.sl.DSstdExp)) respopyr_surv_DSstdExp[i,2]<-ranef2.DSstdExp$mean } # Permutation p-value: Amount of variation expected by chance larger than observed variation in selection table(respopyr_surv_DSstdExp[,1]>ranefPopYrmean)/n.sim table(respopyr_surv_DSstdExp[,2]>ranefslPopYrmean)/n.sim # Save model output saveRDS(respopyr_surv_DSstdExp[,1], file = "YearUniqPermutation_DSstdExp_Intercept_Survival.rds") saveRDS(respopyr_surv_DSstdExp[,2], file = "YearUniqPermutation_DSstdExp_Slope_Survival.rds") # Micro-scale temporal variation (among plot-years) significance: Plot-Year label permutation sformula_Survwodens_DSstdExp_pyRnd = Survived ~ DSstdsqrtExpScoreCor + I(DSstdsqrtExpScoreCor^2) + f(pop, model = "iid2d", n=N.Pop) + f(pop2, DSstdsqrtExpScoreCor, copy = "pop") + f(p, model = "iid2d", n=N.Plot) + f(p2, DSstdsqrtExpScoreCor, copy = "p") + f(yr, model = "iid2d", n=N.Year) + f(yr2, DSstdsqrtExpScoreCor, copy = "yr") + f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, DSstdsqrtExpScoreCor, copy = "popyr") + f(spy, model = "iid2d", n=sN.PY) + f(spy2, DSstdsqrtExpScoreCor, copy = "spy") n.sim <- 1000 respy_surv_DSstdExp<-matrix(NA, n.sim,2) for(i in 1:n.sim){ datPopComb$sPlotYear<-sample(datPopComb$PlotYear,nrow(datPopComb))## randomize plot labels sn.PY = nlevels(datPopComb$sPlotYear) sN.PY = 2*sn.PY spy = as.numeric(datPopComb$sPlotYear) spy2 = as.numeric(datPopComb$sPlotYear) + sn.PY datPopComb$spy = spy datPopComb$spy2 = spy2 modSurv_DSstdExp_pyS.Inla = inla (sformula_Survwodens_DSstdExp_pyRnd, data = datPopComb, family = "binomial", control.predictor=list(compute=T)) tau.int.PY.DSstdExp = modSurv_DSstdExp_pyS.Inla$internal.marginals.hyperpar$"Log precision for spy (component 1)" s.marginal.DSstdExp <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PY.DSstdExp) ranef1.DSstdExp = as.data.frame (inla.zmarginal(s.marginal.DSstdExp)) respy_surv_DSstdExp[i,1]<-ranef1.DSstdExp$mean tau.slope.PY.DSstdExp = modSurv_DSstdExp_pyS.Inla$internal.marginals.hyperpar$"Log precision for spy (component 2)" s.marginal.sl.DSstdExp <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PY.DSstdExp) ranef2.DSstdExp = as.data.frame (inla.zmarginal(s.marginal.DSstdExp)) respy_surv_DSstdExp[i,2]<-ranef2.DSstdExp$mean } # Permutation p-value: Amount of variation expected by chance larger than observed variation in selection table(respy_surv_DSstdExp[,1]>ranefPYmean)/n.sim table(respy_surv_DSstdExp[,2]>ranefslPYmean)/n.sim # Save model output saveRDS(respy_surv_DSstdExp[,1], file = "PlotYearPermutation_DSstdExp_Intercept_Survival.rds") saveRDS(respy_surv_DSstdExp[,2], file = "PlotYearPermutation_DSstdExp_Slope_Survival.rds") # Table 2 (Survival): Proportion of variance attributable to each hierarchical level (i.e., random factor) # total variance slopes totvar = s.marginal.slpy + s.marginal.slp + s.marginal.slpopyr + s.marginal.slpop + s.marginal.slyr # Macro-spatial variation: proportion of variance in slope explained by population r.slpop = s.marginal.slpop/totvar r.slpop.mean = mean(r.slpop[,1]) ci.pop = as.data.frame(quantile(r.slpop[,1], prob = c(0.025, 0.975))) r.slpop.df = paste0(round(r.slpop.mean,2)," ", "(",round(ci.pop[1,],2), ",", round(ci.pop[2,],2),")") # Micro-spatial variation: proportion of variance in slope explained by plot r.slp = s.marginal.slp/totvar r.slp.mean = mean(r.slp[,1]) ci.p = as.data.frame(quantile(r.slp[,1], prob = c(0.025, 0.975))) r.slp.df = paste0(round(r.slp.mean,2)," ", "(",round(ci.p[1,],2), ",", round(ci.p[2,],2),")") # Temporal variation: proportion of variance in slope explained by year r.slyr = s.marginal.slyr/totvar r.slyr.mean = mean(r.slyr[,1]) ci.yr = as.data.frame(quantile(r.slyr[,1], prob = c(0.025, 0.975))) r.slyr.df = paste0(round(r.slyr.mean,2)," ", "(",round(ci.yr[1,],2), ",", round(ci.yr[2,],2),")") # Macro-scale temporal variation: proportion of variance in slope explained by population-year r.slpopyr = s.marginal.slpopyr/totvar r.slpopyr.mean = mean(r.slpopyr[,1]) ci.popyr = as.data.frame(quantile(r.slpopyr[,1], prob = c(0.025, 0.975))) r.slpopyr.df = paste0(round(r.slpopyr.mean,2)," ", "(",round(ci.popyr[1,],2), ",", round(ci.popyr[2,],2),")") # Micro-scale temporal variation: proportion of variance in slope explained by plot-year r.slpy = s.marginal.slpy/totvar r.slpy.mean = mean(r.slpy[,1]) ci.py = as.data.frame(quantile(r.slpy[,1], prob = c(0.025, 0.975))) r.slpy.df = paste0(round(r.slpy.mean,2)," ", "(",round(ci.py[1,],2), ",", round(ci.py[2,],2),")") # Gather all proportions together into one dataframe r.all.df = as.data.frame(rbind(r.slpop.df, r.slp.df, r.slyr.df, r.slpopyr.df, r.slpy.df)) # Table S3: Linear and nonlinear personality-related density-dependent unstandardized selection gradients---- # Integrative fitness---- # Model formula_densInt3_DSstd_IndFit = AnnualFitness ~ DSstdsqrtExpScoreCor + I(DSstdsqrtExpScoreCor^2) + densCenDSMean + PopDensCen + densCenDSMean:PopDensCen + DSstdsqrtExpScoreCor:densCenDSMean + I(DSstdsqrtExpScoreCor^2):densCenDSMean + DSstdsqrtExpScoreCor:PopDensCen + I(DSstdsqrtExpScoreCor^2):PopDensCen + DSstdsqrtExpScoreCor:densCenDSMean:PopDensCen + I(DSstdsqrtExpScoreCor^2):densCenDSMean:PopDensCen + f(pop, model = "iid2d", n=N.Pop) + f(pop2, DSstdsqrtExpScoreCor, copy = "pop") + f(p, model = "iid2d", n=N.Plot) + f(p2, DSstdsqrtExpScoreCor, copy = "p") + f(yr, model = "iid2d", n=N.Year) + f(yr2, DSstdsqrtExpScoreCor, copy = "yr") + f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, DSstdsqrtExpScoreCor, copy = "popyr") + f(py, model = "iid2d", n = N.PY) + f(py2, DSstdsqrtExpScoreCor, copy = "py") modINLA_densInt3_DSstd_IndFit = inla (formula_densInt3_DSstd_IndFit, data = datPopComb, family = "gaussian", control.compute=list(config = TRUE), control.predictor=list(compute=T)) summary(modINLA_densInt3_DSstd_IndFit) saveRDS(modINLA_densInt3_DSstd_IndFit, file = "INLAmodel_DensDepInt3_IndFitness_DSstd.rds") # Table S3 (Integrative fitness): Linear and nonlinear unstandardized selection gradient fixefTable = modINLA_densInt3_DSstd_IndFit$summary.fixed # Table S3 (Integrative fitness): Variance (intercept and slope) and intercept-slope correlation estimates # For macro-spatial variation: among populations tau.int.Pop= modINLA_densInt3_DSstd_IndFit$internal.marginals.hyperpar$"Log precision for pop (component 1)" s.marginal.pop <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Pop) ranef.int.Pop = as.data.frame (inla.zmarginal(s.marginal.pop)) ranefPopmean = ranef.int.Pop$mean tau.slope.Pop= modINLA_densInt3_DSstd_IndFit$internal.marginals.hyperpar$"Log precision for pop (component 2)" s.marginal.slpop <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Pop) ranef.slope.Pop = as.data.frame (inla.zmarginal(s.marginal.slpop)) ranefslPopmean = ranef.slope.Pop$mean rho.int.Pop= modINLA_densInt3_DSstd_IndFit$internal.marginals.hyperpar$"Rho_internal1:2 for pop" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Pop) ranef.rho.int.Pop = as.data.frame (inla.zmarginal(s.marginal)) # For micro-spatial variation: among plots tau.int.Plot= modINLA_densInt3_DSstd_IndFit$internal.marginals.hyperpar$"Log precision for p (component 1)" s.marginal.p <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Plot) ranef.int.Plot = as.data.frame (inla.zmarginal(s.marginal.p)) ranefPlotmean = ranef.int.Plot$mean tau.slope.Plot= modINLA_densInt3_DSstd_IndFit$internal.marginals.hyperpar$"Log precision for p (component 2)" s.marginal.slp <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Plot) ranef.slope.Plot = as.data.frame (inla.zmarginal(s.marginal.slp)) ranefslPlotmean = ranef.slope.Plot$mean rho.int.Plot= modINLA_densInt3_DSstd_IndFit$internal.marginals.hyperpar$"Rho_internal1:2 for p" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Plot) ranef.rho.int.Plot = as.data.frame (inla.zmarginal(s.marginal)) # For temporal variation: among years tau.int.Yr= modINLA_densInt3_DSstd_IndFit$internal.marginals.hyperpar$"Log precision for yr (component 1)" s.marginal.yr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Yr) ranef.int.Yr = as.data.frame (inla.zmarginal(s.marginal.yr)) ranefYrmean = ranef.int.Yr$mean tau.slope.Yr= modINLA_densInt3_DSstd_IndFit$internal.marginals.hyperpar$"Log precision for yr (component 2)" s.marginal.slyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Yr) ranef.slope.Yr = as.data.frame (inla.zmarginal(s.marginal.slyr)) ranefslYrmean = ranef.slope.Yr$mean rho.int.Yr= modINLA_densInt3_DSstd_IndFit$internal.marginals.hyperpar$"Rho_internal1:2 for yr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Yr) ranef.rho.int.Yr = as.data.frame (inla.zmarginal(s.marginal)) # For macro-scale temporal variation: among population-years tau.int.PopYr= modINLA_densInt3_DSstd_IndFit$internal.marginals.hyperpar$"Log precision for popyr (component 1)" s.marginal.popyr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PopYr) ranef.int.PopYr = as.data.frame (inla.zmarginal(s.marginal.popyr)) ranefPopYrmean = ranef.int.PopYr$mean tau.slope.PopYr= modINLA_densInt3_DSstd_IndFit$internal.marginals.hyperpar$"Log precision for popyr (component 2)" s.marginal.slpopyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PopYr) ranef.slope.PopYr = as.data.frame (inla.zmarginal(s.marginal.slpopyr)) ranefslPopYrmean = ranef.slope.PopYr$mean rho.int.PopYr= modINLA_densInt3_DSstd_IndFit$internal.marginals.hyperpar$"Rho_internal1:2 for popyr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PopYr) ranef.rho.int.PopYr = as.data.frame (inla.zmarginal(s.marginal)) # For micro-scale temporal variation: among plot-years tau.int.PY= modINLA_densInt3_DSstd_IndFit$internal.marginals.hyperpar$"Log precision for py (component 1)" s.marginal.py <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PY) ranef.int.PY = as.data.frame (inla.zmarginal(s.marginal.py)) ranefPYmean = ranef.int.PY$mean tau.slope.PY= modINLA_densInt3_DSstd_IndFit$internal.marginals.hyperpar$"Log precision for py (component 2)" s.marginal.slpy <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PY) ranef.slope.PY = as.data.frame (inla.zmarginal(s.marginal.slpy)) ranefslPYmean = ranef.slope.PY$mean rho.int.PY= modINLA_densInt3_DSstd_IndFit$internal.marginals.hyperpar$"Rho_internal1:2 for py" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PY) ranef.rho.int.PY = as.data.frame (inla.zmarginal(s.marginal)) # Residuals tau.int.resi= modINLA_densInt3_DSstd_IndFit$internal.marginals.hyperpar$"Log precision for Gaussian observations" s.marginal.resi <- inla.tmarginal(function(x) (1/exp(x)), tau.int.resi) ranef.int.resi = as.data.frame (inla.zmarginal(s.marginal.resi)) # Offspring recruitment---- # Define alpha and beta parameters for poisson model apar <- 0.5 bpar = apar*var(datPopComb$AnnualRecruits) lgprior <- list(prec = list(prior="loggamma", param = c(apar,bpar))) # Model formula_densInt3_DSstd = AnnualRecruits ~ DSstdsqrtExpScoreCor + I(DSstdsqrtExpScoreCor^2) + densCenDSMean + PopDensCen + densCenDSMean:PopDensCen + DSstdsqrtExpScoreCor:densCenDSMean + I(DSstdsqrtExpScoreCor^2):densCenDSMean + DSstdsqrtExpScoreCor:PopDensCen + I(DSstdsqrtExpScoreCor^2):PopDensCen + DSstdsqrtExpScoreCor:densCenDSMean:PopDensCen + I(DSstdsqrtExpScoreCor^2):densCenDSMean:PopDensCen + f(pop, model = "iid2d", n=N.Pop) + f(pop2, DSstdsqrtExpScoreCor, copy = "pop") + f(p, model = "iid2d", n=N.Plot) + f(p2, DSstdsqrtExpScoreCor, copy = "p") + f(yr, model = "iid2d", n=N.Year) + f(yr2, DSstdsqrtExpScoreCor, copy = "yr") + f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, DSstdsqrtExpScoreCor, copy = "popyr") + f(py, model = "iid2d", n = N.PY) + f(py2, DSstdsqrtExpScoreCor, copy = "py") + f(obsid, model = "iid", hyper = lgprior) modINLA_densInt3_DSstd = inla (formula_densInt3_DSstd, data = datPopComb, family = "poisson", control.compute=list(cpo=TRUE, config = TRUE), control.predictor=list(compute=T)) # Check if there are any failures of the model (data point for which estimation poorly performed) sum(modINLA_densInt3_DSstd$cpo$failure) # recompute model failures manually modINLA_densInt3_DSstd = inla.cpo(modINLA_densInt3_DSstd) summary(modINLA_densInt3_DSstd) saveRDS(modINLA_densInt3_DSstd, file = "INLAmodel_DensDepInt3_Recruit_DSstd.rds") # Table S3 (Recruitment): Linear and nonlinear unstandardized selection gradient fixefTable = modINLA_densInt3_DSstd_IndFit$summary.fixed # Table S3 (Recruitment): Variance (intercept and slope) and intercept-slope correlation estimates # For macro-spatial variation: among populations tau.int.Pop= modINLA_densInt3_DSstd$internal.marginals.hyperpar$"Log precision for pop (component 1)" s.marginal.pop <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Pop) ranef.int.Pop = as.data.frame (inla.zmarginal(s.marginal.pop)) ranefPopmean = ranef.int.Pop$mean tau.slope.Pop= modINLA_densInt3_DSstd$internal.marginals.hyperpar$"Log precision for pop (component 2)" s.marginal.slpop <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Pop) ranef.slope.Pop = as.data.frame (inla.zmarginal(s.marginal.slpop)) ranefslPopmean = ranef.slope.Pop$mean rho.int.Pop= modINLA_densInt3_DSstd$internal.marginals.hyperpar$"Rho_internal1:2 for pop" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Pop) ranef.rho.int.Pop = as.data.frame (inla.zmarginal(s.marginal)) # For micro-spatial variation: among plots tau.int.Plot= modINLA_densInt3_DSstd$internal.marginals.hyperpar$"Log precision for p (component 1)" s.marginal.p <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Plot) ranef.int.Plot = as.data.frame (inla.zmarginal(s.marginal.p)) ranefPlotmean = ranef.int.Plot$mean tau.slope.Plot= modINLA_densInt3_DSstd$internal.marginals.hyperpar$"Log precision for p (component 2)" s.marginal.slp <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Plot) ranef.slope.Plot = as.data.frame (inla.zmarginal(s.marginal.slp)) ranefslPlotmean = ranef.slope.Plot$mean rho.int.Plot= modINLA_densInt3_DSstd$internal.marginals.hyperpar$"Rho_internal1:2 for p" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Plot) ranef.rho.int.Plot = as.data.frame (inla.zmarginal(s.marginal)) # For temporal variation: among years tau.int.Yr= modINLA_densInt3_DSstd$internal.marginals.hyperpar$"Log precision for yr (component 1)" s.marginal.yr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Yr) ranef.int.Yr = as.data.frame (inla.zmarginal(s.marginal.yr)) ranefYrmean = ranef.int.Yr$mean tau.slope.Yr= modINLA_densInt3_DSstd$internal.marginals.hyperpar$"Log precision for yr (component 2)" s.marginal.slyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Yr) ranef.slope.Yr = as.data.frame (inla.zmarginal(s.marginal.slyr)) ranefslYrmean = ranef.slope.Yr$mean rho.int.Yr= modINLA_densInt3_DSstd$internal.marginals.hyperpar$"Rho_internal1:2 for yr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Yr) ranef.rho.int.Yr = as.data.frame (inla.zmarginal(s.marginal)) # For macro-scale temporal variation: among population-years tau.int.PopYr= modINLA_densInt3_DSstd$internal.marginals.hyperpar$"Log precision for popyr (component 1)" s.marginal.popyr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PopYr) ranef.int.PopYr = as.data.frame (inla.zmarginal(s.marginal.popyr)) ranefPopYrmean = ranef.int.PopYr$mean tau.slope.PopYr= modINLA_densInt3_DSstd$internal.marginals.hyperpar$"Log precision for popyr (component 2)" s.marginal.slpopyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PopYr) ranef.slope.PopYr = as.data.frame (inla.zmarginal(s.marginal.slpopyr)) ranefslPopYrmean = ranef.slope.PopYr$mean rho.int.PopYr= modINLA_densInt3_DSstd$internal.marginals.hyperpar$"Rho_internal1:2 for popyr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PopYr) ranef.rho.int.PopYr = as.data.frame (inla.zmarginal(s.marginal)) # For micro-scale temporal variation: among plot-years tau.int.PY= modINLA_densInt3_DSstd$internal.marginals.hyperpar$"Log precision for py (component 1)" s.marginal.py <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PY) ranef.int.PY = as.data.frame (inla.zmarginal(s.marginal.py)) ranefPYmean = ranef.int.PY$mean tau.slope.PY= modINLA_densInt3_DSstd$internal.marginals.hyperpar$"Log precision for py (component 2)" s.marginal.slpy <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PY) ranef.slope.PY = as.data.frame (inla.zmarginal(s.marginal.slpy)) ranefslPYmean = ranef.slope.PY$mean rho.int.PY= modINLA_densInt3_DSstd$internal.marginals.hyperpar$"Rho_internal1:2 for py" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PY) ranef.rho.int.PY = as.data.frame (inla.zmarginal(s.marginal)) # Overdispersion estimate tau.int.obsid= modINLA_densInt3_DSstd$internal.marginals.hyperpar$"Log precision for obsid" s.marginal.obsid <- inla.tmarginal(function(x) (1/exp(x)), tau.int.obsid) ranef.int.obsid = as.data.frame (inla.zmarginal(s.marginal.obsid)) # Adult dispersal---- formula_densInt3_DSstd_surv = Survived ~ DSstdsqrtExpScoreCor + I(DSstdsqrtExpScoreCor^2) + densCenDSMean + PopDensCen + densCenDSMean:PopDensCen + DSstdsqrtExpScoreCor:densCenDSMean + I(DSstdsqrtExpScoreCor^2):densCenDSMean + DSstdsqrtExpScoreCor:PopDensCen + I(DSstdsqrtExpScoreCor^2):PopDensCen + DSstdsqrtExpScoreCor:densCenDSMean:PopDensCen + I(DSstdsqrtExpScoreCor^2):densCenDSMean:PopDensCen + f(pop, model = "iid2d", n=N.Pop) + f(pop2, DSstdsqrtExpScoreCor, copy = "pop") + f(p, model = "iid2d", n=N.Plot) + f(p2, DSstdsqrtExpScoreCor, copy = "p") + f(yr, model = "iid2d", n=N.Year) + f(yr2, DSstdsqrtExpScoreCor, copy = "yr") + f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, DSstdsqrtExpScoreCor, copy = "popyr") + f(py, model = "iid2d", n = N.PY) + f(py2, DSstdsqrtExpScoreCor, copy = "py") modINLA_densInt3_DSstd_surv = inla (formula_densInt3_DSstd_surv, data = datPopComb, family = "binomial", control.family = list(link = "logit"), control.compute=list(cpo=TRUE, config=TRUE), control.inla = list(strategy = "laplace", npoints = 21), control.predictor=list(compute=T)) # Check if there are any failures of the model (data point for which estimation poorly performed) sum(modINLA_densInt3_DSstd_surv$cpo$failure) # recompute model failures manually modINLA_densInt3_DSstd_surv = inla.cpo(modINLA_densInt3_DSstd_surv) summary(modINLA_densInt3_DSstd_surv) saveRDS(modINLA_densInt3_DSstd_surv, file = "INLAmodel_DensDepInt3_Survival_DSstd.rds") # Table S3 (Survival): Linear and nonlinear unstandardized selection gradient fixefTable = modINLA_densInt3_DSstd_surv$summary.fixed # Table S3 (Survival): Variance (intercept and slope) and intercept-slope correlation estimates # For macro-spatial variation: among populations tau.int.Pop= modINLA_densInt3_DSstd_surv$internal.marginals.hyperpar$"Log precision for pop (component 1)" s.marginal.pop <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Pop) ranef.int.Pop = as.data.frame (inla.zmarginal(s.marginal.pop)) ranefPopmean = ranef.int.Pop$mean tau.slope.Pop= modINLA_densInt3_DSstd_surv$internal.marginals.hyperpar$"Log precision for pop (component 2)" s.marginal.slpop <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Pop) ranef.slope.Pop = as.data.frame (inla.zmarginal(s.marginal.slpop)) ranefslPopmean = ranef.slope.Pop$mean rho.int.Pop= modINLA_densInt3_DSstd_surv$internal.marginals.hyperpar$"Rho_internal1:2 for pop" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Pop) ranef.rho.int.Pop = as.data.frame (inla.zmarginal(s.marginal)) # For micro-spatial variation: among plots tau.int.Plot= modINLA_densInt3_DSstd_surv$internal.marginals.hyperpar$"Log precision for p (component 1)" s.marginal.p <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Plot) ranef.int.Plot = as.data.frame (inla.zmarginal(s.marginal.p)) ranefPlotmean = ranef.int.Plot$mean tau.slope.Plot= modINLA_densInt3_DSstd_surv$internal.marginals.hyperpar$"Log precision for p (component 2)" s.marginal.slp <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Plot) ranef.slope.Plot = as.data.frame (inla.zmarginal(s.marginal.slp)) ranefslPlotmean = ranef.slope.Plot$mean rho.int.Plot= modINLA_densInt3_DSstd_surv$internal.marginals.hyperpar$"Rho_internal1:2 for p" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Plot) ranef.rho.int.Plot = as.data.frame (inla.zmarginal(s.marginal)) # For temporal variation: among years tau.int.Yr= modINLA_densInt3_DSstd_surv$internal.marginals.hyperpar$"Log precision for yr (component 1)" s.marginal.yr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Yr) ranef.int.Yr = as.data.frame (inla.zmarginal(s.marginal.yr)) ranefYrmean = ranef.int.Yr$mean tau.slope.Yr= modINLA_densInt3_DSstd_surv$internal.marginals.hyperpar$"Log precision for yr (component 2)" s.marginal.slyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Yr) ranef.slope.Yr = as.data.frame (inla.zmarginal(s.marginal.slyr)) ranefslYrmean = ranef.slope.Yr$mean rho.int.Yr= modINLA_densInt3_DSstd_surv$internal.marginals.hyperpar$"Rho_internal1:2 for yr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Yr) ranef.rho.int.Yr = as.data.frame (inla.zmarginal(s.marginal)) # For macro-scale temporal variation: among population-years tau.int.PopYr= modINLA_densInt3_DSstd_surv$internal.marginals.hyperpar$"Log precision for popyr (component 1)" s.marginal.popyr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PopYr) ranef.int.PopYr = as.data.frame (inla.zmarginal(s.marginal.popyr)) ranefPopYrmean = ranef.int.PopYr$mean tau.slope.PopYr= modINLA_densInt3_DSstd_surv$internal.marginals.hyperpar$"Log precision for popyr (component 2)" s.marginal.slpopyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PopYr) ranef.slope.PopYr = as.data.frame (inla.zmarginal(s.marginal.slpopyr)) ranefslPopYrmean = ranef.slope.PopYr$mean rho.int.PopYr= modINLA_densInt3_DSstd_surv$internal.marginals.hyperpar$"Rho_internal1:2 for popyr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PopYr) ranef.rho.int.PopYr = as.data.frame (inla.zmarginal(s.marginal)) # For micro-scale temporal variation: among plot-years tau.int.PY= modINLA_densInt3_DSstd_surv$internal.marginals.hyperpar$"Log precision for py (component 1)" s.marginal.py <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PY) ranef.int.PY = as.data.frame (inla.zmarginal(s.marginal.py)) ranefPYmean = ranef.int.PY$mean tau.slope.PY= modINLA_densInt3_DSstd_surv$internal.marginals.hyperpar$"Log precision for py (component 2)" s.marginal.slpy <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PY) ranef.slope.PY = as.data.frame (inla.zmarginal(s.marginal.slpy)) ranefslPYmean = ranef.slope.PY$mean rho.int.PY= modINLA_densInt3_DSstd_surv$internal.marginals.hyperpar$"Rho_internal1:2 for py" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PY) ranef.rho.int.PY = as.data.frame (inla.zmarginal(s.marginal)) ## Table S4: Model controlling for variation in phenotypes among hierarchical levels (i.e., randam factors)---- ## = Include BLUPs for all levels as fixed effects ## BLUP (= mode of the posterior distribution) estimation and extraction # Prior of the model sdExp = sd(datPopComb$DSstdsqrtExpScoreCor) pcprior <- list(prec = list(prior="pc.prec", param = c((sdExp/0.31), 0.01))) # Model to estimate BLUPs---- formula_BLUP = DSstdsqrtExpScoreCor ~ f(pop, model = "iid") + f(p, model = "iid", hyper = pcprior) + f(yr, model = "iid", hyper = pcprior) + f(popyr, model = "iid", hyper = pcprior) + f(py, model = "iid", hyper = pcprior) modINLA_BLUP = inla (formula_BLUP, data = datPopComb, family = "gaussian", control.compute=list(cpo = TRUE, config = TRUE), control.predictor=list(compute=T)) summary(modINLA_BLUP) # Extract BLUP for each level ID # BLUP for population blup.pop = modINLA_BLUP$summary.random$pop blup.pop2 = data.frame(cbind(pop = blup.pop$ID, BlupMode.pop = blup.pop$mode)) # BLUP for plot blup.p = modINLA_BLUP$summary.random$p blup.p2 = as.data.frame(cbind(p = blup.p$ID, BlupMode.p = blup.p$mode)) # BLUP for year blup.yr = modINLA_BLUP$summary.random$yr blup.yr2 = data.frame(cbind(yr = blup.yr$ID, BlupMode.yr = blup.yr$mode)) # BLUP for Population-Year blup.popyr = modINLA_BLUP$summary.random$popyr blup.popyr2 = data.frame(cbind(popyr = blup.popyr$ID, BlupMode.popyr = blup.popyr$mode)) # BLUP for plot-year blup.py = modINLA_BLUP$summary.random$py blup.py2 = as.data.frame(cbind(py = blup.py$ID, BlupMode.py = blup.py$mode)) # Add BLUPs to dataset datPopComb2 = merge(datPopComb, blup.py2, by = "py", all.x =TRUE) datPopComb2 = merge(datPopComb2, blup.p2, by = "p", all.x =TRUE) datPopComb2 = merge(datPopComb2, blup.popyr2, by = "popyr", all.x =TRUE) datPopComb2 = merge(datPopComb2, blup.pop2, by = "pop", all.x =TRUE) datPopComb2 = merge(datPopComb2, blup.yr2, by = "yr", all.x =TRUE) # Models with BLUPs---- # Integrative fitness---- # Model formula_IndFit_BLUP = AnnualFitness ~ BlupMode.pop + BlupMode.p + BlupMode.yr + BlupMode.popyr + BlupMode.py + f(pop, model = "iid2d", n=N.Pop) + f(pop2, BlupMode.pop, copy = "pop") + f(p, model = "iid2d", n=N.Plot) + f(p2, BlupMode.p, copy = "p") + f(yr, model = "iid2d", n=N.Year) + f(yr2, BlupMode.yr, copy = "yr") + f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, BlupMode.popyr, copy = "popyr") + f(py, model = "iid2d", n=N.PY) + f(py2, BlupMode.py, copy = "py") modINLA_IndFit_BLUP = inla (formula_IndFit_BLUP, data = datPopComb2, family = "gaussian", control.compute=list(cpo = TRUE, config = TRUE), control.predictor=list(compute=T)) summary(modINLA_IndFit_BLUP) saveRDS(modINLA_IndFit_BLUP, file = "INLAmodel_BLUP_IndFitness.rds") # Table S4 (Integrative fitness): Linear and nonlinear selection gradient fixefTable = modINLA_IndFit_BLUP$summary.fixed # Table S4 (Integrative fitness): Variance (intercept and slope) and intercept-slope correlation estimates # For macro-spatial variation: among populations tau.int.Pop= modINLA_IndFit_BLUP$internal.marginals.hyperpar$"Log precision for pop (component 1)" s.marginal.pop <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Pop) ranef.int.Pop = as.data.frame (inla.zmarginal(s.marginal.pop)) ranefPopmean = ranef.int.Pop$mean tau.slope.Pop= modINLA_IndFit_BLUP$internal.marginals.hyperpar$"Log precision for pop (component 2)" s.marginal.slpop <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Pop) ranef.slope.Pop = as.data.frame (inla.zmarginal(s.marginal.slpop)) ranefslPopmean = ranef.slope.Pop$mean rho.int.Pop= modINLA_IndFit_BLUP$internal.marginals.hyperpar$"Rho_internal1:2 for pop" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Pop) ranef.rho.int.Pop = as.data.frame (inla.zmarginal(s.marginal)) # For micro-spatial variation: among plots tau.int.Plot= modINLA_IndFit_BLUP$internal.marginals.hyperpar$"Log precision for p (component 1)" s.marginal.p <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Plot) ranef.int.Plot = as.data.frame (inla.zmarginal(s.marginal.p)) ranefPlotmean = ranef.int.Plot$mean tau.slope.Plot= modINLA_IndFit_BLUP$internal.marginals.hyperpar$"Log precision for p (component 2)" s.marginal.slp <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Plot) ranef.slope.Plot = as.data.frame (inla.zmarginal(s.marginal.slp)) ranefslPlotmean = ranef.slope.Plot$mean rho.int.Plot= modINLA_IndFit_BLUP$internal.marginals.hyperpar$"Rho_internal1:2 for p" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Plot) ranef.rho.int.Plot = as.data.frame (inla.zmarginal(s.marginal)) tau.int.PopYr= modINLA_IndFit_BLUP$internal.marginals.hyperpar$"Log precision for popyr (component 1)" s.marginal.popyr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PopYr) ranef.int.PopYr = as.data.frame (inla.zmarginal(s.marginal.popyr)) ranefPopYrmean = ranef.int.PopYr$mean tau.slope.PopYr= modINLA_IndFit_BLUP$internal.marginals.hyperpar$"Log precision for popyr (component 2)" s.marginal.slpopyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PopYr) ranef.slope.PopYr = as.data.frame (inla.zmarginal(s.marginal.slpopyr)) ranefslPopYrmean = ranef.slope.PopYr$mean rho.int.PopYr= modINLA_IndFit_BLUP$internal.marginals.hyperpar$"Rho_internal1:2 for popyr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PopYr) ranef.rho.int.PopYr = as.data.frame (inla.zmarginal(s.marginal)) tau.int.Yr = modINLA_IndFit_BLUP$internal.marginals.hyperpar$"Log precision for yr (component 1)" s.marginal.yr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Yr) ranef.int.Yr = as.data.frame (inla.zmarginal(s.marginal.yr)) ranefYrmean = ranef.int.Yr$mean tau.slope.Yr= modINLA_IndFit_BLUP$internal.marginals.hyperpar$"Log precision for yr (component 2)" s.marginal.slyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Yr) ranef.slope.Yr = as.data.frame (inla.zmarginal(s.marginal.slyr)) ranefslYrmean = ranef.slope.Yr$mean rho.int.Yr= modINLA_IndFit_BLUP$internal.marginals.hyperpar$"Rho_internal1:2 for yr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Yr) ranef.rho.int.Yr = as.data.frame (inla.zmarginal(s.marginal)) tau.int.PY= modINLA_IndFit_BLUP$internal.marginals.hyperpar$"Log precision for py (component 1)" s.marginal.py <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PY) ranef.int.PY = as.data.frame (inla.zmarginal(s.marginal.py)) ranefPYmean = ranef.int.PY$mean tau.slope.PY= modINLA_IndFit_BLUP$internal.marginals.hyperpar$"Log precision for py (component 2)" s.marginal.slpy <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PY) ranef.slope.PY = as.data.frame (inla.zmarginal(s.marginal.slpy)) ranefslPYmean = ranef.slope.PY$mean rho.int.PY= modINLA_IndFit_BLUP$internal.marginals.hyperpar$"Rho_internal1:2 for py" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PY) ranef.rho.int.PY = as.data.frame (inla.zmarginal(s.marginal)) # Residual tau.gausobs= modINLA_IndFit_BLUP$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)) # Offspring recruitment---- # Define alpha and beta parameters of poisson model apar <- 0.5 bpar = apar*var(datPopComb2$AnnualRecruits) lgprior <- list(prec = list(prior="loggamma", param = c(apar,bpar))) # Model formula_Recruitwodens_BLUP = AnnualRecruits ~ DSstdsqrtExpScoreCor + I(DSstdsqrtExpScoreCor^2) + BlupMode.py + BlupMode.p + BlupMode.popyr + BlupMode.pop + BlupMode.yr + I(BlupMode.py^2) + I(BlupMode.p^2) + I(BlupMode.popyr^2) + I(BlupMode.pop^2) + I(BlupMode.yr^2) + f(pop, model = "iid2d", n=N.Pop) + f(pop2, DSstdsqrtExpScoreCor, copy = "pop") + f(p, model = "iid2d", n=N.Plot) + f(p2, DSstdsqrtExpScoreCor, copy = "p") + f(yr, model = "iid2d", n=N.Year) + f(yr2, DSstdsqrtExpScoreCor, copy = "yr") f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, DSstdsqrtExpScoreCor, copy = "popyr") + f(py, model = "iid2d", n=N.PY) + f(py2, DSstdsqrtExpScoreCor, copy = "py") + f(obsid, model = "iid", hyper = lgprior) modINLA_Recruitwodens_BLUP = inla (formula_Recruitwodens_BLUP, data = datPopComb2, family = "poisson", control.inla = list(strategy = "laplace", npoints = 21)) # Check if there are any failures of the model (data point for which estimation poorly performed) sum(modINLA_Recruitwodens_BLUP$cpo$failure) summary(modINLA_Recruitwodens_BLUP) saveRDS(modINLA_Recruitwodens_BLUP, file = "INLAmodel_BLUP_Recruit.rds") # Table S4 (Recruitment): Linear and nonlinear selection gradient fixefTable = modINLA_Recruitwodens_BLUP$summary.fixed # Table S4 (Recruitment): Variance (intercept and slope) and intercept-slope correlation estimates # For macro-spatial variation: among populations tau.int.Pop= modINLA_Recruitwodens_BLUP$internal.marginals.hyperpar$"Log precision for pop (component 1)" s.marginal.pop <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Pop) ranef.int.Pop = as.data.frame (inla.zmarginal(s.marginal.pop)) ranefPopmean = ranef.int.Pop$mean tau.slope.Pop= modINLA_Recruitwodens_BLUP$internal.marginals.hyperpar$"Log precision for pop (component 2)" s.marginal.slpop <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Pop) ranef.slope.Pop = as.data.frame (inla.zmarginal(s.marginal.slpop)) ranefslPopmean = ranef.slope.Pop$mean rho.int.Pop= modINLA_Recruitwodens_BLUP$internal.marginals.hyperpar$"Rho_internal1:2 for pop" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Pop) ranef.rho.int.Pop = as.data.frame (inla.zmarginal(s.marginal)) # For micro-spatial variation: among plots tau.int.Plot= modINLA_Recruitwodens_BLUP$internal.marginals.hyperpar$"Log precision for p (component 1)" s.marginal.p <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Plot) ranef.int.Plot = as.data.frame (inla.zmarginal(s.marginal.p)) ranefPlotmean = ranef.int.Plot$mean tau.slope.Plot= modINLA_Recruitwodens_BLUP$internal.marginals.hyperpar$"Log precision for p (component 2)" s.marginal.slp <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Plot) ranef.slope.Plot = as.data.frame (inla.zmarginal(s.marginal.slp)) ranefslPlotmean = ranef.slope.Plot$mean rho.int.Plot= modINLA_Recruitwodens_BLUP$internal.marginals.hyperpar$"Rho_internal1:2 for p" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Plot) ranef.rho.int.Plot = as.data.frame (inla.zmarginal(s.marginal)) # For temporal variation: among years tau.int.Yr = modINLA_Recruitwodens_BLUP$internal.marginals.hyperpar$"Log precision for yr (component 1)" s.marginal.yr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Yr) ranef.int.Yr = as.data.frame (inla.zmarginal(s.marginal.yr)) ranefYrmean = ranef.int.Yr$mean tau.slope.Yr= modINLA_Recruitwodens_BLUP$internal.marginals.hyperpar$"Log precision for yr (component 2)" s.marginal.slyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Yr) ranef.slope.Yr = as.data.frame (inla.zmarginal(s.marginal.slyr)) ranefslYrmean = ranef.slope.Yr$mean rho.int.Yr= modINLA_Recruitwodens_BLUP$internal.marginals.hyperpar$"Rho_internal1:2 for yr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Yr) ranef.rho.int.Yr = as.data.frame (inla.zmarginal(s.marginal)) # For macro-scale temporal variation: among population-years tau.int.PopYr= modINLA_Recruitwodens_BLUP$internal.marginals.hyperpar$"Log precision for popyr (component 1)" s.marginal.popyr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PopYr) ranef.int.PopYr = as.data.frame (inla.zmarginal(s.marginal.popyr)) ranefPopYrmean = ranef.int.PopYr$mean tau.slope.PopYr= modINLA_Recruitwodens_BLUP$internal.marginals.hyperpar$"Log precision for popyr (component 2)" s.marginal.slpopyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PopYr) ranef.slope.PopYr = as.data.frame (inla.zmarginal(s.marginal.slpopyr)) ranefslPopYrmean = ranef.slope.PopYr$mean rho.int.PopYr= modINLA_Recruitwodens_BLUP$internal.marginals.hyperpar$"Rho_internal1:2 for popyr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PopYr) ranef.rho.int.PopYr = as.data.frame (inla.zmarginal(s.marginal)) # For micro-scale temporal variation: among plot-years tau.int.PY= modINLA_Recruitwodens_BLUP$internal.marginals.hyperpar$"Log precision for py (component 1)" s.marginal.py <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PY) ranef.int.PY = as.data.frame (inla.zmarginal(s.marginal.py)) ranefPYmean = ranef.int.PY$mean tau.slope.PY= modINLA_Recruitwodens_BLUP$internal.marginals.hyperpar$"Log precision for py (component 2)" s.marginal.slpy <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PY) ranef.slope.PY = as.data.frame (inla.zmarginal(s.marginal.slpy)) ranefslPYmean = ranef.slope.PY$mean mode.slPY = inla.mmarginal(modINLA_Recruitwodens_BLUP$internal.marginals.hyperpar$"Log precision for py (component 2)") rho.int.PY= modINLA_Recruitwodens_BLUP$internal.marginals.hyperpar$"Rho_internal1:2 for py" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PY) ranef.rho.int.PY = as.data.frame (inla.zmarginal(s.marginal)) # Adult survival---- # Model formula_Survwodens_BLUP = Survived ~ DSstdsqrtExpScoreCor + BlupMode.py + BlupMode.p + BlupMode.popyr + BlupMode.pop + BlupMode.yr + I(DSstdsqrtExpScoreCor^2) + I(BlupMode.py^2) + I(BlupMode.p^2) + I(BlupMode.popyr^2) + I(BlupMode.pop^2) + I(BlupMode.yr^2) + f(pop, model = "iid2d", n=N.Pop) + f(pop2, BlupMode.pop, copy = "pop") + f(p, model = "iid2d", n=N.Plot) + f(p2, BlupMode.p, copy = "p") + f(yr, model = "iid2d", n=N.Year) + f(yr2, BlupMode.yr, copy = "yr") + f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, BlupMode.popyr, copy = "popyr") + f(py, model = "iid2d", n=N.PY) + f(py2, BlupMode.py, copy = "py") modINLA_Survwodens_BLUP = inla (formula_Survwodens_BLUP, data = datPopComb2, family = "binomial", control.family=list(link = "logit")) summary(modINLA_Survwodens_BLUP) # Create string with both estimate and CIs in one column fixefTable = modINLA_Survwodens_BLUP$summary.fixed # Table S4 (Survival): Variance (intercept and slope) and intercept-slope correlation estimates # For macro-spatial variation: among populations tau.int.Pop= modINLA_Survwodens_BLUP$internal.marginals.hyperpar$"Log precision for pop (component 1)" s.marginal.pop <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Pop) ranef.int.Pop = as.data.frame (inla.zmarginal(s.marginal.pop)) ranefPopmean = ranef.int.Pop$mean tau.slope.Pop= modINLA_Survwodens_BLUP$internal.marginals.hyperpar$"Log precision for pop (component 2)" s.marginal.slpop <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Pop) ranef.slope.Pop = as.data.frame (inla.zmarginal(s.marginal.slpop)) ranefslPopmean = ranef.slope.Pop$mean rho.int.Pop= modINLA_Survwodens_BLUP$internal.marginals.hyperpar$"Rho_internal1:2 for pop" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Pop) ranef.rho.int.Pop = as.data.frame (inla.zmarginal(s.marginal)) # For micro-spatial variation: among plots tau.int.Plot= modINLA_Survwodens_BLUP$internal.marginals.hyperpar$"Log precision for p (component 1)" s.marginal.p <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Plot) ranef.int.Plot = as.data.frame (inla.zmarginal(s.marginal.p)) ranefPlotmean = ranef.int.Plot$mean tau.slope.Plot= modINLA_Survwodens_BLUP$internal.marginals.hyperpar$"Log precision for p (component 2)" s.marginal.slp <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Plot) ranef.slope.Plot = as.data.frame (inla.zmarginal(s.marginal.slp)) ranefslPlotmean = ranef.slope.Plot$mean rho.int.Plot= modINLA_Survwodens_BLUP$internal.marginals.hyperpar$"Rho_internal1:2 for p" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Plot) ranef.rho.int.Plot = as.data.frame (inla.zmarginal(s.marginal)) # For temporal variation: among years tau.int.Yr = modINLA_Survwodens_BLUP$internal.marginals.hyperpar$"Log precision for yr (component 1)" s.marginal.yr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Yr) ranef.int.Yr = as.data.frame (inla.zmarginal(s.marginal.yr)) ranefYrmean = ranef.int.Yr$mean tau.slope.Yr= modINLA_Survwodens_BLUP$internal.marginals.hyperpar$"Log precision for yr (component 2)" s.marginal.slyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Yr) ranef.slope.Yr = as.data.frame (inla.zmarginal(s.marginal.slyr)) ranefslYrmean = ranef.slope.Yr$mean rho.int.Yr= modINLA_Survwodens_BLUP$internal.marginals.hyperpar$"Rho_internal1:2 for yr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Yr) ranef.rho.int.Yr = as.data.frame (inla.zmarginal(s.marginal)) # For macro-scale temporal variation: among population-years tau.int.PopYr= modINLA_Survwodens_BLUP$internal.marginals.hyperpar$"Log precision for popyr (component 1)" s.marginal.popyr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PopYr) ranef.int.PopYr = as.data.frame (inla.zmarginal(s.marginal.popyr)) ranefPopYrmean = ranef.int.PopYr$mean tau.slope.PopYr= modINLA_Survwodens_BLUP$internal.marginals.hyperpar$"Log precision for popyr (component 2)" s.marginal.slpopyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PopYr) ranef.slope.PopYr = as.data.frame (inla.zmarginal(s.marginal.slpopyr)) ranefslPopYrmean = ranef.slope.PopYr$mean rho.int.PopYr= modINLA_Survwodens_BLUP$internal.marginals.hyperpar$"Rho_internal1:2 for popyr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PopYr) ranef.rho.int.PopYr = as.data.frame (inla.zmarginal(s.marginal)) # For micro-scale temporal variation: among plot-years tau.int.PY= modINLA_Survwodens_BLUP$internal.marginals.hyperpar$"Log precision for py (component 1)" s.marginal.py <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PY) ranef.int.PY = as.data.frame (inla.zmarginal(s.marginal.py)) ranefPYmean = ranef.int.PY$mean tau.slope.PY= modINLA_Survwodens_BLUP$internal.marginals.hyperpar$"Log precision for py (component 2)" s.marginal.slpy <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PY) ranef.slope.PY = as.data.frame (inla.zmarginal(s.marginal.slpy)) ranefslPYmean = ranef.slope.PY$mean rho.int.PY= modINLA_Survwodens_BLUP$internal.marginals.hyperpar$"Rho_internal1:2 for py" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PY) ranef.rho.int.PY = as.data.frame (inla.zmarginal(s.marginal)) # Table S5: Linear and nonlinear selection gradients standardized at plot-year level---- # Integrative fitness---- formula_IndFit_PYstd = AnnualFitness ~ PYstdsqrtExpScoreCor + I(PYstdsqrtExpScoreCor^2) + f(pop, model = "iid2d", n=N.Pop) + f(pop2, PYstdsqrtExpScoreCor, copy = "pop") + f(p, model = "iid2d", n=N.Plot) + f(p2, PYstdsqrtExpScoreCor, copy = "p") + f(yr, model = "iid2d", n=N.Year) + f(yr2, PYstdsqrtExpScoreCor, copy = "yr") + f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, PYstdsqrtExpScoreCor, copy = "popyr") + f(py, model = "iid2d", n=N.PY) + f(py2, PYstdsqrtExpScoreCor, copy = "py") modINLA_IndFit_PYstd = inla (formula_IndFit_PYstd, data = datPopComb, family = "gaussian", control.compute=list(cpo = TRUE, config = TRUE), control.predictor=list(compute=T)) summary(modINLA_IndFit_PYstd) saveRDS(modINLA_IndFit_PYstd, file = "INLAmodel_PYStd_IndFitness.rds") # Table S4 (Integrative fitness): Linear and nonlinear selection gradient fixefTable = modINLA_IndFit_PYstd$summary.fixed # Table S4 (Integrative fitness): Variance (intercept and slope) and intercept-slope correlation estimates # For macro-spatial variation: among populations tau.int.Pop= modINLA_IndFit_PYstd$internal.marginals.hyperpar$"Log precision for pop (component 1)" s.marginal.pop <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Pop) ranef.int.Pop = as.data.frame (inla.zmarginal(s.marginal.pop)) ranefPopmean = ranef.int.Pop$mean tau.slope.Pop= modINLA_IndFit_PYstd$internal.marginals.hyperpar$"Log precision for pop (component 2)" s.marginal.slpop <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Pop) ranef.slope.Pop = as.data.frame (inla.zmarginal(s.marginal.slpop)) ranefslPopmean = ranef.slope.Pop$mean rho.int.Pop= modINLA_IndFit_PYstd$internal.marginals.hyperpar$"Rho_internal1:2 for pop" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Pop) ranef.rho.int.Pop = as.data.frame (inla.zmarginal(s.marginal)) # For micro-spatial variation: among plots tau.int.Plot= modINLA_IndFit_PYstd$internal.marginals.hyperpar$"Log precision for p (component 1)" s.marginal.p <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Plot) ranef.int.Plot = as.data.frame (inla.zmarginal(s.marginal.p)) ranefPlotmean = ranef.int.Plot$mean tau.slope.Plot= modINLA_IndFit_PYstd$internal.marginals.hyperpar$"Log precision for p (component 2)" s.marginal.slp <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Plot) ranef.slope.Plot = as.data.frame (inla.zmarginal(s.marginal.slp)) ranefslPlotmean = ranef.slope.Plot$mean rho.int.Plot= modINLA_IndFit_PYstd$internal.marginals.hyperpar$"Rho_internal1:2 for p" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Plot) ranef.rho.int.Plot = as.data.frame (inla.zmarginal(s.marginal)) # For temporal variation: among years tau.int.Yr = modINLA_IndFit_PYstd$internal.marginals.hyperpar$"Log precision for yr (component 1)" s.marginal.yr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Yr) ranef.int.Yr = as.data.frame (inla.zmarginal(s.marginal.yr)) ranefYrmean = ranef.int.Yr$mean tau.slope.Yr= modINLA_IndFit_PYstd$internal.marginals.hyperpar$"Log precision for yr (component 2)" s.marginal.slyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Yr) ranef.slope.Yr = as.data.frame (inla.zmarginal(s.marginal.slyr)) ranefslYrmean = ranef.slope.Yr$mean rho.int.Yr= modINLA_IndFit_PYstd$internal.marginals.hyperpar$"Rho_internal1:2 for yr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Yr) ranef.rho.int.Yr = as.data.frame (inla.zmarginal(s.marginal)) # For macro-scale temporal variation: among population-years tau.int.PopYr= modINLA_IndFit_PYstd$internal.marginals.hyperpar$"Log precision for popyr (component 1)" s.marginal.popyr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PopYr) ranef.int.PopYr = as.data.frame (inla.zmarginal(s.marginal.popyr)) ranefPopYrmean = ranef.int.PopYr$mean tau.slope.PopYr= modINLA_IndFit_PYstd$internal.marginals.hyperpar$"Log precision for popyr (component 2)" s.marginal.slpopyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PopYr) ranef.slope.PopYr = as.data.frame (inla.zmarginal(s.marginal.slpopyr)) ranefslPopYrmean = ranef.slope.PopYr$mean rho.int.PopYr= modINLA_IndFit_PYstd$internal.marginals.hyperpar$"Rho_internal1:2 for popyr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PopYr) ranef.rho.int.PopYr = as.data.frame (inla.zmarginal(s.marginal)) # For micro-scale temporal variation: among plot-years tau.int.PY= modINLA_IndFit_PYstd$internal.marginals.hyperpar$"Log precision for py (component 1)" s.marginal.py <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PY) ranef.int.PY = as.data.frame (inla.zmarginal(s.marginal.py)) ranefPYmean = ranef.int.PY$mean tau.slope.PY= modINLA_IndFit_PYstd$internal.marginals.hyperpar$"Log precision for py (component 2)" s.marginal.slpy <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PY) ranef.slope.PY = as.data.frame (inla.zmarginal(s.marginal.slpy)) ranefslPYmean = ranef.slope.PY$mean rho.int.PY= modINLA_IndFit_PYstd$internal.marginals.hyperpar$"Rho_internal1:2 for py" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PY) ranef.rho.int.PY = as.data.frame (inla.zmarginal(s.marginal)) # Residual tau.resi = modINLA_IndFit_PYstd$internal.marginals.hyperpar$"Log precision for the Gaussian observations" s.marginal.resi <- inla.tmarginal(function(x) (1/exp(x)), tau.resi) ranef.resi = as.data.frame (inla.zmarginal(s.marginal.resi)) # Offspring recruitment---- ## Model with exploration standardised at the Plot-Year level # Define alpha and beta parameters for poisson model apar <- 0.5 bpar = apar*var(datPopComb$AnnualRecruits) lgprior <- list(prec = list(prior="loggamma", param = c(apar,bpar))) # Model formula_rndquad_wodens = AnnualRecruits ~ PYstdsqrtExpScoreCor + I(PYstdsqrtExpScoreCor^2) + f(pop, model = "iid2d", n=N.Pop) + f(pop2, PYstdsqrtExpScoreCor, copy = "pop") + f(p, model = "iid2d", n=N.Plot) + f(p2, PYstdsqrtExpScoreCor, copy = "p") + f(syr, model = "iid2d", n=sN.Year) + f(syr2, PYstdsqrtExpScoreCor, copy = "syr")+ f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, PYstdsqrtExpScoreCor, copy = "popyr") + f(py, model = "iid2d", n = N.PY) + f(py2, PYstdsqrtExpScoreCor, copy = "py") + f(obsid, model = "iid", hyper = lgprior) modINLArndquad_wodens = inla (formula_rndquad_wodens, data = datPopComb, family = "poisson", control.compute=list(cpo=TRUE, config=TRUE), control.inla = list(strategy = "laplace", npoints = 21), control.predictor=list(compute=T)) # Check if there are any failures of the model (data point for which estimation poorly performed) sum(modINLArndquad_wodens$cpo$failure) # recompute model failures manually modINLArndquad_wodens = inla.cpo(modINLArndquad_wodens) summary(modINLArndquad_wodens) # Table S5 (Recruitment): Linear and nonlinear selection gradients fixefTable = modINLA_wRec_PYstd$summary.fixed # Table S5 (Recruitment): Variance (intercept and slope) and intercept-slope correlation estimates # For macro-spatial variation: among populations tau.int.Pop= modINLArndquad_wodens$internal.marginals.hyperpar$"Log precision for pop (component 1)" s.marginal.pop <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Pop) ranef.int.Pop = as.data.frame (inla.zmarginal(s.marginal.pop)) tau.slope.Pop= modINLArndquad_wodens$internal.marginals.hyperpar$"Log precision for pop (component 2)" s.marginal.slpop <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Pop) ranef.slope.Pop = as.data.frame (inla.zmarginal(s.marginal.slpop)) rho.int.Pop= modINLArndquad_wodens$internal.marginals.hyperpar$"Rho_internal1:2 for pop" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Pop) ranef.rho.int.Pop = as.data.frame (inla.zmarginal(s.marginal)) tau.int.Plot= modINLArndquad_wodens$internal.marginals.hyperpar$"Log precision for p (component 1)" s.marginal.p <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Plot) ranef.int.Plot = as.data.frame (inla.zmarginal(s.marginal.p)) tau.slope.Plot= modINLArndquad_wodens$internal.marginals.hyperpar$"Log precision for p (component 2)" s.marginal.slp <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Plot) ranef.slope.Plot = as.data.frame (inla.zmarginal(s.marginal.slp)) rho.int.Plot= modINLArndquad_wodens$internal.marginals.hyperpar$"Rho_internal1:2 for p" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Plot) ranef.rho.int.Plot = as.data.frame (inla.zmarginal(s.marginal)) tau.int.PopYr= modINLArndquad_wodens$internal.marginals.hyperpar$"Log precision for popyr (component 1)" s.marginal.popyr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PopYr) ranef.int.PopYr = as.data.frame (inla.zmarginal(s.marginal.popyr)) tau.slope.PopYr= modINLArndquad_wodens$internal.marginals.hyperpar$"Log precision for popyr (component 2)" s.marginal.slpopyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PopYr) ranef.slope.PopYr = as.data.frame (inla.zmarginal(s.marginal.slpopyr)) rho.int.PopYr= modINLArndquad_wodens$internal.marginals.hyperpar$"Rho_internal1:2 for popyr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PopYr) ranef.rho.int.PopYr = as.data.frame (inla.zmarginal(s.marginal)) tau.int.Yr= modINLArndquad_wodens$internal.marginals.hyperpar$"Log precision for yr (component 1)" s.marginal.yr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Yr) ranef.int.Yr = as.data.frame (inla.zmarginal(s.marginal.yr)) tau.slope.Yr= modINLArndquad_wodens$internal.marginals.hyperpar$"Log precision for yr (component 2)" s.marginal.slyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Yr) ranef.slope.Yr = as.data.frame (inla.zmarginal(s.marginal.slyr)) rho.int.Yr= modINLArndquad_wodens$internal.marginals.hyperpar$"Rho_internal1:2 for yr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Yr) ranef.rho.int.Yr = as.data.frame (inla.zmarginal(s.marginal)) # For micro-scale temporal variation: among plot-years tau.int.PY= modINLArndquad_wodens$internal.marginals.hyperpar$"Log precision for py (component 1)" s.marginal.py <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PY) ranef.int.PY = as.data.frame (inla.zmarginal(s.marginal.py)) tau.slope.PY= modINLArndquad_wodens$internal.marginals.hyperpar$"Log precision for py (component 2)" s.marginal.slpy <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PY) ranef.slope.PY = as.data.frame (inla.zmarginal(s.marginal.slpy)) rho.int.PY= modINLArndquad_wodens$internal.marginals.hyperpar$"Rho_internal1:2 for py" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PY) ranef.rho.int.PY = as.data.frame (inla.zmarginal(s.marginal)) # Overdispersion estimate tau.int.obsid= modINLArndquad_wodens$internal.marginals.hyperpar$"Log precision for obsid" s.marginal.obsid <- inla.tmarginal(function(x) (1/exp(x)), tau.int.obsid) ranef.int.obsid = as.data.frame (inla.zmarginal(s.marginal.obsid)) # Adult survival---- formula_PYstd_Survival = Survived ~ PYstdsqrtExpScoreCor + I(PYstdsqrtExpScoreCor^2) + f(pop, model = "iid2d", n=N.Pop) + f(pop2, PYstdsqrtExpScoreCor, copy = "pop") + f(p, model = "iid2d", n=N.Plot) + f(p2, PYstdsqrtExpScoreCor, copy = "p") + f(yr, model = "iid2d", n=N.Year) + f(yr2, PYstdsqrtExpScoreCor, copy = "yr") + f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, PYstdsqrtExpScoreCor, copy = "popyr") + f(py, model = "iid2d", n = N.PY) + f(py2, PYstdsqrtExpScoreCor, copy = "py") modINLA_PYstd_Survival = inla (formula_PYstd_Survival, data = datPopComb, family = "binomial", control.family = list(link = "logit"), control.compute=list(cpo=TRUE, config=TRUE), control.inla = list(strategy = "laplace", npoints = 21), control.predictor=list(compute=T)) # Check if there are any failures of the model (data point for which estimation poorly performed) sum(modINLA_PYstd_Survival$cpo$failure) # recompute model failures manually modINLA_PYstd_Survival = inla.cpo(modINLA_PYstd_Survival) summary(modINLA_PYstd_Survival) # Table S5 (Survival): Linear and nonlinear selection gradients fixefTable = modINLA_amgPop_woDens_quad$summary.fixed # Table S5 (Survival): Variance (intercept and slope) and intercept-slope correlation estimates # For macro-spatial variation: among populations tau.int.Pop= modINLA_PYstd_Survival$internal.marginals.hyperpar$"Log precision for pop (component 1)" s.marginal.pop <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Pop) ranef.int.Pop = as.data.frame (inla.zmarginal(s.marginal.pop)) ranefPopmean.PYstd = ranef.int.Pop$mean tau.slope.Pop= modINLA_PYstd_Survival$internal.marginals.hyperpar$"Log precision for pop (component 2)" s.marginal.slpop <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Pop) ranef.slope.Pop = as.data.frame (inla.zmarginal(s.marginal.slpop)) ranefslPopmean.PYstd = ranef.slope.Pop$mean rho.int.Pop= modINLA_PYstd_Survival$internal.marginals.hyperpar$"Rho_internal1:2 for pop" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Pop) ranef.rho.int.Pop = as.data.frame (inla.zmarginal(s.marginal)) # For micro-spatial variation: among plots tau.int.Plot= modINLA_PYstd_Survival$internal.marginals.hyperpar$"Log precision for p (component 1)" s.marginal.p <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Plot) ranef.int.Plot = as.data.frame (inla.zmarginal(s.marginal.p)) ranefPlotmean.PYstd = ranef.int.Plot$mean tau.slope.Plot= modINLA_PYstd_Survival$internal.marginals.hyperpar$"Log precision for p (component 2)" s.marginal.slp <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Plot) ranef.slope.Plot = as.data.frame (inla.zmarginal(s.marginal.slp)) ranefslPlotmean.PYstd = ranef.slope.Plot$mean rho.int.Plot= modINLA_PYstd_Survival$internal.marginals.hyperpar$"Rho_internal1:2 for p" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Plot) ranef.rho.int.Plot = as.data.frame (inla.zmarginal(s.marginal)) # For temporal variation: among years tau.int.Yr = modINLA_PYstd_Survival$internal.marginals.hyperpar$"Log precision for yr (component 1)" s.marginal.yr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Yr) ranef.int.Yr = as.data.frame (inla.zmarginal(s.marginal.yr)) ranefYrmean.PYstd = ranef.int.Yr$mean tau.slope.Yr = modINLA_PYstd_Survival$internal.marginals.hyperpar$"Log precision for yr (component 2)" s.marginal.slyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Yr) ranef.slope.Yr = as.data.frame (inla.zmarginal(s.marginal.slyr)) ranefslYrmean.PYstd = ranef.slope.Yr$mean rho.int.Yr = modINLA_PYstd_Survival$internal.marginals.hyperpar$"Rho_internal1:2 for yr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Yr) ranef.rho.int.Yr = as.data.frame (inla.zmarginal(s.marginal)) # For macro-scale temporal variation: among population-years tau.int.PopYr = modINLA_PYstd_Survival$internal.marginals.hyperpar$"Log precision for popyr (component 1)" s.marginal.popyr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PopYr) ranef.int.PopYr = as.data.frame (inla.zmarginal(s.marginal.popyr)) ranefPopYrmean.PYstd = ranef.int.PopYr$mean tau.slope.PopYr = modINLA_PYstd_Survival$internal.marginals.hyperpar$"Log precision for popyr (component 2)" s.marginal.slpopyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PopYr) ranef.slope.PopYr = as.data.frame (inla.zmarginal(s.marginal.slpopyr)) ranefslPopYrmean.PYstd = ranef.slope.PopYr$mean rho.int.PopYr = modINLA_PYstd_Survival$internal.marginals.hyperpar$"Rho_internal1:2 for popyr" s.marginal.popyr <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PopYr) ranef.rho.int.PopYr = as.data.frame (inla.zmarginal(s.marginal.popyr)) # For micro-scale temporal variation: among plot-years tau.int.PY= modINLA_PYstd_Survival$internal.marginals.hyperpar$"Log precision for py (component 1)" s.marginal.py <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PY) ranef.int.PY = as.data.frame (inla.zmarginal(s.marginal.py)) ranefPYmean.PYstd = ranef.int.PY$mean tau.slope.PY= modINLA_PYstd_Survival$internal.marginals.hyperpar$"Log precision for py (component 2)" s.marginal.slpy <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PY) ranef.slope.PY = as.data.frame (inla.zmarginal(s.marginal.slpy)) ranefslPYmean.PYstd = ranef.slope.PY$mean rho.int.PY= modINLA_PYstd_Survival$internal.marginals.hyperpar$"Rho_internal1:2 for py" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PY) ranef.rho.int.PY = as.data.frame (inla.zmarginal(s.marginal)) # Table S6 (Fledglings): Variance (intercept and slope) and intercept-slope correlation estimates---- # Model for integrative fitness with annual number fledglings ---- # Exploration standardized over whole DS # Model formula_FitFled_DSstdExp = AnnualFitnessFled ~ DSstdsqrtExpScoreCor + I(DSstdsqrtExpScoreCor^2) + f(pop, model = "iid2d", n=N.Pop) + f(pop2, DSstdsqrtExpScoreCor, copy = "pop") + f(p, model = "iid2d", n=N.Plot) + f(p2, DSstdsqrtExpScoreCor, copy = "p") + f(yr, model = "iid2d", n=N.Year) + f(yr2, DSstdsqrtExpScoreCor, copy = "yr") + f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, DSstdsqrtExpScoreCor, copy = "popyr") + f(py, model = "iid2d", n=N.PY) + f(py2, DSstdsqrtExpScoreCor, copy = "py") modINLA_FitFled_DSstdExp = inla (formula_FitFled_DSstdExp, data = datPopComb, family = "gaussian", control.compute=list(dic = TRUE, cpo = TRUE, waic = T, config = TRUE), control.predictor=list(compute=T)) summary(modINLA_FitFled_DSstdExp) saveRDS(modINLA_FitFled_DSstdExp, file = "INLAmodel_DSStd_IndFitnessFled.rds") # Retrieve variances and correlations of random factors # For macro-spatial variation: among populations tau.int.Pop= modINLA_FitFled_DSstdExp$internal.marginals.hyperpar$"Log precision for pop (component 1)" s.marginal.pop <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Pop) ranef.int.Pop = as.data.frame (inla.zmarginal(s.marginal.pop)) ranefPopmean = ranef.int.Pop$mean tau.slope.Pop= modINLA_FitFled_DSstdExp$internal.marginals.hyperpar$"Log precision for pop (component 2)" s.marginal.slpop <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Pop) ranef.slope.Pop = as.data.frame (inla.zmarginal(s.marginal.slpop)) ranefslPopmean = ranef.slope.Pop$mean rho.int.Pop= modINLA_FitFled_DSstdExp$internal.marginals.hyperpar$"Rho_internal1:2 for pop" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Pop) ranef.rho.int.Pop = as.data.frame (inla.zmarginal(s.marginal)) # For micro-spatial variation: among plots tau.int.Plot= modINLA_FitFled_DSstdExp$internal.marginals.hyperpar$"Log precision for p (component 1)" s.marginal.p <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Plot) ranef.int.Plot = as.data.frame (inla.zmarginal(s.marginal.p)) ranefPlotmean = ranef.int.Plot$mean tau.slope.Plot= modINLA_FitFled_DSstdExp$internal.marginals.hyperpar$"Log precision for p (component 2)" s.marginal.slp <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Plot) ranef.slope.Plot = as.data.frame (inla.zmarginal(s.marginal.slp)) ranefslPlotmean = ranef.slope.Plot$mean rho.int.Plot= modINLA_FitFled_DSstdExp$internal.marginals.hyperpar$"Rho_internal1:2 for p" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Plot) ranef.rho.int.Plot = as.data.frame (inla.zmarginal(s.marginal)) # For temporal variation: among years tau.int.Yr = modINLA_FitFled_DSstdExp$internal.marginals.hyperpar$"Log precision for yr (component 1)" s.marginal.yr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Yr) ranef.int.Yr = as.data.frame (inla.zmarginal(s.marginal.yr)) ranefYrmean = ranef.int.Yr$mean tau.slope.Yr= modINLA_FitFled_DSstdExp$internal.marginals.hyperpar$"Log precision for yr (component 2)" s.marginal.slyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Yr) ranef.slope.Yr = as.data.frame (inla.zmarginal(s.marginal.slyr)) ranefslYrmean = ranef.slope.Yr$mean rho.int.Yr= modINLA_FitFled_DSstdExp$internal.marginals.hyperpar$"Rho_internal1:2 for yr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Yr) ranef.rho.int.Yr = as.data.frame (inla.zmarginal(s.marginal)) # For macro-scale temporal variation: among population-years tau.int.PopYr= modINLA_FitFled_DSstdExp$internal.marginals.hyperpar$"Log precision for popyr (component 1)" s.marginal.popyr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PopYr) ranef.int.PopYr = as.data.frame (inla.zmarginal(s.marginal.popyr)) ranefPopYrmean = ranef.int.PopYr$mean tau.slope.PopYr= modINLA_FitFled_DSstdExp$internal.marginals.hyperpar$"Log precision for popyr (component 2)" s.marginal.slpopyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PopYr) ranef.slope.PopYr = as.data.frame (inla.zmarginal(s.marginal.slpopyr)) ranefslPopYrmean = ranef.slope.PopYr$mean rho.int.PopYr= modINLA_FitFled_DSstdExp$internal.marginals.hyperpar$"Rho_internal1:2 for popyr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PopYr) ranef.rho.int.PopYr = as.data.frame (inla.zmarginal(s.marginal)) # For micro-scale temporal variation: among plot-years tau.int.PY= modINLA_FitFled_DSstdExp$internal.marginals.hyperpar$"Log precision for py (component 1)" s.marginal.py <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PY) ranef.int.PY = as.data.frame (inla.zmarginal(s.marginal.py)) ranefPYmean = ranef.int.PY$mean tau.slope.PY= modINLA_FitFled_DSstdExp$internal.marginals.hyperpar$"Log precision for py (component 2)" s.marginal.slpy <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PY) ranef.slope.PY = as.data.frame (inla.zmarginal(s.marginal.slpy)) ranefslPYmean = ranef.slope.PY$mean rho.int.PY= modINLA_FitFled_DSstdExp$internal.marginals.hyperpar$"Rho_internal1:2 for py" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PY) ranef.rho.int.PY = as.data.frame (inla.zmarginal(s.marginal)) # Residual tau.resi= modINLA_FitFled_DSstdExp$internal.marginals.hyperpar$"Log precision for the Gaussian observations" s.marginal.resi <- inla.tmarginal(function(x) (1/exp(x)), tau.resi) ranef.resi = as.data.frame (inla.zmarginal(s.marginal.resi)) # Estimates for fixed effects fixefTable = modINLA_FitwodensFled_DSstdExp$summary.fixed # Proportion in variance for all ecological levels---- # total variance for slopes totvar = s.marginal.slpop + s.marginal.slp + s.marginal.slyr + s.marginal.slpopyr + s.marginal.slpy # for population r.slpop = s.marginal.slpop/totvar r.slpop.mean = mean(r.slpop[,1]) ci.pop = as.data.frame(quantile(r.slpop[,1], prob = c(0.025, 0.975))) r.slpop.df = paste0(round(r.slpop.mean,2)," ", "(",round(ci.pop[1,],2), ",", round(ci.pop[2,],2),")") # for plot r.slp = s.marginal.slp/totvar r.slp.mean = mean(r.slp[,1]) ci.p = as.data.frame(quantile(r.slp[,1], prob = c(0.025, 0.975))) r.slp.df = paste0(round(r.slp.mean,2)," ", "(",round(ci.p[1,],2), ",", round(ci.p[2,],2),")") # for year r.slyr = s.marginal.slyr/totvar r.slyr.mean = mean(r.slyr[,1]) ci.yr = as.data.frame(quantile(r.slyr[,1], prob = c(0.025, 0.975))) r.slyr.df = paste0(round(r.slyr.mean,2)," ", "(",round(ci.yr[1,],2), ",", round(ci.yr[2,],2),")") # for population-year r.slpopyr = s.marginal.slpopyr/totvar r.slpopyr.mean = mean(r.slpopyr[,1]) ci.popyr = as.data.frame(quantile(r.slpopyr[,1], prob = c(0.025, 0.975))) r.slpopyr.df = paste0(round(r.slpopyr.mean,2)," ", "(",round(ci.popyr[1,],2), ",", round(ci.popyr[2,],2),")") # for plot-year r.slpy = s.marginal.slpy/totvar r.slpy.mean = mean(r.slpy[,1]) ci.py = as.data.frame(quantile(r.slpy[,1], prob = c(0.025, 0.975))) r.slpy.df = paste0(round(r.slpy.mean,2)," ", "(",round(ci.py[1,],2), ",", round(ci.py[2,],2),")") # Gather all proportions of variances in one data frame r.all.df = as.data.frame(rbind(r.slpop.df, r.slp.df, r.slyr.df, r.slpopyr.df, r.slpy.df)) # Model for number of fledglings as fitness metric---- # Exploration standardized over whole DS # Define alpha and beta parameters for poisson model apar <- 0.5 bpar = apar*var(datPopComb$AnnualFledged) lgprior <- list(prec = list(prior="loggamma", param = c(apar,bpar))) # Model formula_Fled_DSstd = AnnualFledged ~ DSstdsqrtExpScoreCor + I(DSstdsqrtExpScoreCor^2) + f(pop, model = "iid2d", n=N.Pop) + f(pop2, DSstdsqrtExpScoreCor, copy = "pop") + f(p, model = "iid2d", n=N.Plot) + f(p2, DSstdsqrtExpScoreCor, copy = "p") + f(yr, model = "iid2d", n=N.Year) + f(yr2, DSstdsqrtExpScoreCor, copy = "yr") + f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, DSstdsqrtExpScoreCor, copy = "popyr") + f(py, model = "iid2d", n = N.PY) + f(py2, DSstdsqrtExpScoreCor, copy = "py") f(obsid, model = "iid", hyper = lgprior) modINLA_Fled_DSstd = inla (formula_Fled_DSstd, data = datPopComb, family = "poisson", control.compute=list(cpo=TRUE), control.predictor=list(compute=T)) # Check if there are any failures of the model (data point for which estimation poorly performed) sum(modINLA_Fled_DSstd$cpo$failure) # recompute model failures manually modINLA_Fled_DSstd = inla.cpo(modINLA_Fled_DSstd) summary(modINLA_Fled_DSstd) saveRDS(modINLA_Fled_DSstd, file = "INLAmodel_DSStd_Fled.rds") # Retrieve variances and correlations of random factors # For macro-spatial variation: among populations tau.int.Pop= modINLA_Fled_DSstd$internal.marginals.hyperpar$"Log precision for pop (component 1)" s.marginal.pop <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Pop) ranef.int.Pop = as.data.frame (inla.zmarginal(s.marginal.pop)) ranefPopmean = ranef.int.Pop$mean tau.slope.Pop= modINLA_Fled_DSstd$internal.marginals.hyperpar$"Log precision for pop (component 2)" s.marginal.slpop <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Pop) ranef.slope.Pop = as.data.frame (inla.zmarginal(s.marginal.slpop)) ranefslPopmean = ranef.slope.Pop$mean rho.int.Pop= modINLA_Fled_DSstd$internal.marginals.hyperpar$"Rho_internal1:2 for pop" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Pop) ranef.rho.int.Pop = as.data.frame (inla.zmarginal(s.marginal)) # For micro-spatial variation: among plots tau.int.Plot= modINLA_Fled_DSstd$internal.marginals.hyperpar$"Log precision for p (component 1)" s.marginal.p <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Plot) ranef.int.Plot = as.data.frame (inla.zmarginal(s.marginal.p)) ranefPlotmean = ranef.int.Plot$mean tau.slope.Plot= modINLA_Fled_DSstd$internal.marginals.hyperpar$"Log precision for p (component 2)" s.marginal.slp <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Plot) ranef.slope.Plot = as.data.frame (inla.zmarginal(s.marginal.slp)) ranefslPlotmean = ranef.slope.Plot$mean rho.int.Plot= modINLA_Fled_DSstd$internal.marginals.hyperpar$"Rho_internal1:2 for p" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Plot) ranef.rho.int.Plot = as.data.frame (inla.zmarginal(s.marginal)) # For temporal variation: among years tau.int.Yr= modINLA_Fled_DSstd$internal.marginals.hyperpar$"Log precision for yr (component 1)" s.marginal.yr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.Yr) ranef.int.Yr = as.data.frame (inla.zmarginal(s.marginal.yr)) ranefYrmean = ranef.int.Yr$mean tau.slope.Yr= modINLA_Fled_DSstd$internal.marginals.hyperpar$"Log precision for yr (component 2)" s.marginal.slyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.Yr) ranef.slope.Yr = as.data.frame (inla.zmarginal(s.marginal.slyr)) ranefslYrmean = ranef.slope.Yr$mean rho.int.Yr= modINLA_Fled_DSstd$internal.marginals.hyperpar$"Rho_internal1:2 for yr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.Yr) ranef.rho.int.Yr = as.data.frame (inla.zmarginal(s.marginal)) # For macro-scale temporal variation: among population-years tau.int.PopYr= modINLA_Fled_DSstd$internal.marginals.hyperpar$"Log precision for popyr (component 1)" s.marginal.popyr <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PopYr) ranef.int.PopYr = as.data.frame (inla.zmarginal(s.marginal.popyr)) ranefPopYrmean = ranef.int.PopYr$mean tau.slope.PopYr= modINLA_Fled_DSstd$internal.marginals.hyperpar$"Log precision for popyr (component 2)" s.marginal.slpopyr <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PopYr) ranef.slope.PopYr = as.data.frame (inla.zmarginal(s.marginal.slpopyr)) ranefslPopYrmean = ranef.slope.PopYr$mean rho.int.PopYr= modINLA_Fled_DSstd$internal.marginals.hyperpar$"Rho_internal1:2 for popyr" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PopYr) ranef.rho.int.PopYr = as.data.frame (inla.zmarginal(s.marginal)) # For micro-scale temporal variation: among plot-years tau.int.PY= modINLA_Fled_DSstd$internal.marginals.hyperpar$"Log precision for py (component 1)" s.marginal.py <- inla.tmarginal(function(x) (1/exp(x)), tau.int.PY) ranef.int.PY = as.data.frame (inla.zmarginal(s.marginal.py)) ranefPYmean = ranef.int.PY$mean tau.slope.PY= modINLA_Fled_DSstd$internal.marginals.hyperpar$"Log precision for py (component 2)" s.marginal.slpy <- inla.tmarginal(function(x) (1/exp(x)), tau.slope.PY) ranef.slope.PY = as.data.frame (inla.zmarginal(s.marginal.slpy)) ranefslPYmean = ranef.slope.PY$mean rho.int.PY= modINLA_Fled_DSstd$internal.marginals.hyperpar$"Rho_internal1:2 for py" s.marginal <- inla.tmarginal(function(x) (2*exp(x)/(1+exp(x))-1), rho.int.PY) ranef.rho.int.PY = as.data.frame (inla.zmarginal(s.marginal)) # Overdispersion tau.int.obsid= modINLA_Fled_DSstd$internal.marginals.hyperpar$"Log precision for obsid" s.marginal.obsid <- inla.tmarginal(function(x) (1/exp(x)), tau.int.obsid) ranef.int.obsid = as.data.frame (inla.zmarginal(s.marginal.obsid)) # Proportion of variances for all ecological levels---- # total variance slopes totvar = s.marginal.slpop + s.marginal.slp + s.marginal.slyr + s.marginal.slpopyr + s.marginal.slpy # for population r.slpop = s.marginal.slpop/totvar r.slpop.mean = mean(r.slpop[,1]) ci.pop = as.data.frame(quantile(r.slpop[,1], prob = c(0.025, 0.975))) r.slpop.df = paste0(round(r.slpop.mean,2)," ", "(",round(ci.pop[1,],2), ",", round(ci.pop[2,],2),")") # for plot r.slp = s.marginal.slp/totvar r.slp.mean = mean(r.slp[,1]) ci.p = as.data.frame(quantile(r.slp[,1], prob = c(0.025, 0.975))) r.slp.df = paste0(round(r.slp.mean,2)," ", "(",round(ci.p[1,],2), ",", round(ci.p[2,],2),")") # for years r.slyr = s.marginal.slyr/totvar r.slyr.mean = mean(r.slyr[,1]) ci.yr = as.data.frame(quantile(r.slyr[,1], prob = c(0.025, 0.975))) r.slyr.df = paste0(round(r.slyr.mean,2)," ", "(",round(ci.yr[1,],2), ",", round(ci.yr[2,],2),")") # for population-year r.slpopyr = s.marginal.slpopyr/totvar r.slpopyr.mean = mean(r.slpopyr[,1]) ci.popyr = as.data.frame(quantile(r.slpopyr[,1], prob = c(0.025, 0.975))) r.slpopyr.df = paste0(round(r.slpopyr.mean,2)," ", "(",round(ci.popyr[1,],2), ",", round(ci.popyr[2,],2),")") # for plot-year r.slpy = s.marginal.slpy/totvar r.slpy.mean = mean(r.slpy[,1]) ci.py = as.data.frame(quantile(r.slpy[,1], prob = c(0.025, 0.975))) r.slpy.df = paste0(round(r.slpy.mean,2)," ", "(",round(ci.py[1,],2), ",", round(ci.py[2,],2),")") # Gather all proportions of variance in one data frame r.all.df = as.data.frame(rbind(r.slpop.df, r.slp.df, r.slyr.df, r.slpopyr.df, r.slpy.df)) # Figure 2---- # Extract selection gradients from model with exploration standardized at plot-year level # to represent gradients estimated a plot-year level and grouped by plot formula_IndFit_PYstd_w = wAnnualFitness ~ PYstdsqrtExpScoreCor + I(PYstdsqrtExpScoreCor^2) + f(pop, model = "iid2d", n=N.Pop) + f(pop2, PYstdsqrtExpScoreCor, copy = "pop") + f(p, model = "iid2d", n=N.Plot) + f(p2, PYstdsqrtExpScoreCor, copy = "p") + f(yr, model = "iid2d", n=N.Year) + f(yr2, PYstdsqrtExpScoreCor, copy = "yr") + f(popyr, model = "iid2d", n=N.PopYear) + f(popyr2, PYstdsqrtExpScoreCor, copy = "popyr") + f(py, model = "iid2d", n=N.PY) + f(py2, PYstdsqrtExpScoreCor, copy = "py") modINLA_IndFit_PYstd_w = inla (formula_IndFit_PYstd_w, data = datPopComb, family = "gaussian", control.compute=list(cpo = TRUE, config = TRUE), control.predictor=list(compute=T)) saveRDS(modINLA_IndFit_PYstd_w, file = "modINLA_IndFit_PYstd_w.rds") modINLA_IndFit_PYstd_w = readRDS("modINLA_IndFit_PYstd_w.rds") fixefTable = modINLA_IndFit_PYstd_w$summary.fixed # Selection gradient for each plot-year selgrd.py = subset(modINLA_IndFit_PYstd_w$summary.random$py2, ID>200) # select estimates for slopes selgrd.py$PYID = 1:nrow(selgrd.py) # rename IDs starting from 1 colnames(selgrd.py) = c("py2", "SelGrd", "sd", "lwrCI", "0.5CI", "uprCI", "mode", "kld", "PYID") # select variables from full dataset to merge with gradients linkPop.py = distinct(dplyr::select(datPopComb, PopulationID, Plot, py2), py2, .keep_all = T) selgrd.py = merge(selgrd.py, linkPop.py, by = "py2") # Remove selection gradient estimate at plot-year level for populations with only one plot selgrd.py$SelGrd[selgrd.py$PopulationID =="4"] = "NA" selgrd.py$SelGrd[selgrd.py$PopulationID =="5"] = "NA" # Create new Plot IDs to be more meaningful and consistent in the article # Data frame with each original Plot ID being unique uniqPlot_selgrd = distinct(dplyr::select(datPopComb, Plot, PopulationID), Plot, .keep_all = T) # Switch Population IDs for they appear in the wanted order in the graph uniqPlot_selgrd$PopulationID = as.factor (uniqPlot_selgrd$PopulationID) levels (uniqPlot_selgrd$PopulationID)[levels(uniqPlot_selgrd$PopulationID) ==c("1","2", "3", "4", "5")] = c("2","1","3", "4", "5") uniqPlot_selgrd$PopulationID = factor (uniqPlot_selgrd$PopulationID, levels = c("1", "2", "3", "4", "5")) uniqPlot_selgrd = uniqPlot_selgrd[order(uniqPlot_selgrd$PopulationID),] # New Plot IDs uniqPlot_selgrd$PlotID = 1:nrow(uniqPlot_selgrd) uniqPlot_selgrd$PlotID = as.factor(uniqPlot_selgrd$PlotID) # Select only columns needed uniqPlot_selgrd = dplyr::select(uniqPlot_selgrd, Plot, PlotID) # Merge new Plot IDs with selection gradient estimates and CIs selgrd.py_plotgraph = merge(uniqPlot_selgrd, selgrd.py, by = "Plot") # Add mean selection gradient estimate (beta exploration from the model) to # plot selection gradients relative to zero rather than to mean gradient selgrd.py_plotgraph$SelGrd = as.numeric(selgrd.py_plotgraph$SelGrd) selgrd.py_plotgraph$SelGrdTrue = selgrd.py_plotgraph$SelGrd + fixefTable$mean[2] # Selection gradients for each population-year for populations with only one plot selgrd.popyr = subset(modINLA_IndFit_PYstd_w$summary.random$popyr2, ID>51) # select slope estimates colnames(selgrd.popyr)[1] = "popyr2" # select variables from full dataset to merge with gradients linkPop.popyr = distinct(dplyr::select(datPopComb, BroodYear, PopYear, PopulationID, Plot, popyr2), popyr2, .keep_all = T) selgrd.popyr = merge(selgrd.popyr, linkPop.popyr, by = "popyr2") # data frame with only populations with one plot selgrd.popyr.whww = selgrd.popyr %>% filter(PopulationID == "4"| PopulationID=="5") colnames(selgrd.popyr.whww) = c("popyr2", "SelGrd", "sd", "lwrCI", "0.5CI", "uprCI", "mode", "kld","BroodYear", "PopYear", "PopulationID", "Plot") # Rename Plot IDs selgrd.popyr.whww$PlotID = c(rep("34", 18), rep("35", 12)) # Add mean selection gradient estimate (beta exploration from the model) to # plot selection gradients relative to zero rather than to mean gradient selgrd.popyr.whww$SelGrdTrue = selgrd.popyr.whww$SelGrd + fixefTable$mean[2] # Plot Figure 2 # Create a vector of x-axis labels (i.e., Plot IDs) for only every other ID appears # Create a vector of length of Plot IDs breaks <- seq(1, 37, 1) # Replace labels of even IDs with blank labels <- as.character(breaks) labels[!(breaks %% 2 != 0)] <- '' breaks = as.factor(breaks) # Plot ggplot(data = selgrd.py_plotgraph, aes(PlotID, SelGrdTrue, colour = PopulationID)) + theme(rect = element_blank(), axis.line = element_line(colour = "black"), legend.position = "none", axis.ticks = element_line(colour = "black"), axis.text.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.title = element_text(size = 12)) + scale_y_continuous (breaks = seq (-0.40, 0.50, 0.1)) + scale_x_discrete(breaks = breaks, labels = labels) + ylab(expression(paste("Linear selection gradient ", " ", (beta)))) + xlab("Plot ID") + geom_boxplot(alpha =0.4) + geom_jitter(alpha =0.3, size = 2, width = 0.2) + geom_boxplot(data = selgrd.popyr.whww,aes(PlotID, SelGrd, colour = PopulationID), # boxplots for populations with only one plot inherit.aes =FALSE, alpha =0.4) + geom_jitter(data = selgrd.popyr.whww,aes(PlotID, SelGrd, colour = PopulationID), inherit.aes =FALSE, alpha =0.3, size = 2, width = 0.2) + scale_colour_manual(values = c("#1f77b4" , "#d62728", "#ff7f0e", "#9467bd", "#e377c2" )) + expand_limits(x = c(1,37)) + annotate("text", label = "Belgium", x = 5, y = 0.35, size = 3) + annotate("text", label = "Germany", x = 15, y = 0.35, size = 3) + annotate("text", label = "Netherlands", x = 28, y = 0.35, size = 3) + annotate("text", label = "UK", x = 35.6, y = 0.35, size = 3) + geom_segment(aes(x = 1, y = 0.33, xend = 9.5, yend = 0.33), colour = "black") + geom_segment(aes(x = 10, y = 0.33, xend = 21.5, yend = 0.33), colour = "black") + geom_segment(aes(x = 22, y = 0.33, xend = 34.5, yend = 0.33), colour = "black") + geom_segment(aes(x = 35, y = 0.33, xend = 36.5, yend = 0.33), colour = "black")