#Packages# library(readr) library(tidyverse) library(plotrix) library(lme4) library(MuMIn) library(gridExtra) library(ggplot2) #age at first reproduction - females# #data_cont_fem = data, subset age_at_first_reproduction for females mod_cont_fem <- glm(age ~ PrimOri + PopSize + F + PartExp + I(PopSize^2), data = data_cont_fem , na.action = "na.fail", family=poisson) summary(mod_cont_fem) dfun <- function(object) { with(object,sum((weights * residuals^2)[weights > 0])/df.residual) } x.quasipoisson <- function(...) { res <- quasipoisson(...) res$aic <- poisson(...)$aic res } test <- update(mod_cont_fem, family="x.quasipoisson", na.action = na.fail) (gg <- dredge(test,rank="QAIC", chat=dfun(mod_cont_fem),subset = dc(PopSize, I(PopSize^2)))) gg_subset <- subset(gg, delta<=2.0) summary(model.avg(gg_subset)) model_cont_fem <- model.avg(get.models(gg, subset=delta<=2.0), fit = TRUE) new_data_cont_fem_nb_terr <- expand.grid(PopSize = seq(0,max(data_cont_fem$PopSize)), F = mean(data_cont_fem$F), PartExp = c("no"), PrimOri = c("yes")) predict_cont_fem_nb_terr <- predict(model_cont_fem , new_data_cont_fem_nb_terr, type = "response", se.fit = T, FULL = TRUE) new_data_cont_fem_nb_terr$Fit <- predict_cont_fem_nb_terr$fit new_data_cont_fem_nb_terr$se <- predict_cont_fem_nb_terr$se.fit Figure_cont_fem_nb_terr <- ggplot() + scale_x_continuous(breaks = seq(0,80,10), limits = c(0,80)) + scale_y_continuous(breaks = seq(1,7,1), limits = c(1,7)) + scale_color_manual("",values=c("#000000","#000000"), guide=F) + geom_ribbon(data = new_data_cont_fem_nb_terr, aes(ymax=Fit+se,ymin=Fit-se, x=PopSize),alpha=0.5, fill=c("grey")) + geom_jitter(data = data_cont_fem, aes(x=PopSize, y=age), size=1, height = 0.1, width = 0) + geom_line(data = new_data_cont_fem_nb_terr, aes(x=PopSize, y=Fit), colour=c("#000000"), size=0.5) + labs(y = "", x = "",size=0.2) + theme(axis.ticks = element_line(colour = "#000000", size = 0.2), axis.line = element_line(colour = "#000000", size = 0.2), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), panel.background = element_blank(), axis.text.x = element_text(size=12), axis.text.y = element_text(size=12), plot.title =element_text(size=12)) + ggtitle("a)") Figure_cont_fem_nb_terr #age at first reproduction - males# #data_cont_male = data, subset age_at_first_reproduction for males mod_cont_male <- glm(age ~ PrimOri + PopSize + F + PartExp + I(PopSize^2), data = data_cont_male , na.action = "na.fail", family=poisson) summary(mod_cont_male) dfun <- function(object) { with(object,sum((weights * residuals^2)[weights > 0])/df.residual) } x.quasipoisson <- function(...) { res <- quasipoisson(...) res$aic <- poisson(...)$aic res } test <- update(mod_cont_male, family="x.quasipoisson", na.action = na.fail) (gg <- dredge(test,rank="QAIC", chat=dfun(mod_cont_male),subset = dc(PopSize, I(PopSize^2)))) gg_subset <- subset(gg, delta<=2.0) summary(model.avg(gg_subset)) model_cont_male <- model.avg(get.models(gg, subset=delta<=2.0), fit = TRUE) new_data_cont_male_nb_terr <- expand.grid(PopSize = seq(0,max(data_cont_male$PopSize)), F = mean(data_cont_male$F), PartExp = c("no"), PrimOri = c("yes")) predict_cont_male_nb_terr <- predict(model_cont_male , new_data_cont_male_nb_terr, type = "response", se.fit = T, FULL = TRUE) new_data_cont_male_nb_terr$Fit <- predict_cont_male_nb_terr$fit new_data_cont_male_nb_terr$se <- predict_cont_male_nb_terr$se.fit Figure_cont_male_nb_terr <- ggplot() + scale_x_continuous(breaks = seq(0,80,10), limits = c(0,80)) + scale_y_continuous(breaks = seq(1,7,1), limits = c(1,7)) + scale_color_manual("",values=c("#000000","#000000"), guide=F) + geom_ribbon(data = new_data_cont_male_nb_terr, aes(ymax=Fit+se,ymin=Fit-se, x=PopSize),alpha=0.5, fill=c("grey")) + geom_jitter(data = data_cont_male, aes(x=PopSize, y=age), size=1, height = 0.1, width = 0) + geom_line(data = new_data_cont_male_nb_terr, aes(x=PopSize, y=Fit), colour=c("#000000"), size=0.5) + labs(y = "", x = "",size=0.2) + theme(axis.ticks = element_line(colour = "#000000", size = 0.2), axis.line = element_line(colour = "#000000", size = 0.2), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), panel.background = element_blank(), axis.text.x = element_text(size=12), axis.text.y = element_text(size=12), plot.title =element_text(size=12)) + ggtitle("b)") Figure_cont_male_nb_terr #Figure 2# Figure_2 <- grid.arrange(Figure_cont_fem_nb_terr, Figure_cont_male_nb_terr, ncol = 2, left = "Age at first reproduction (years)", bottom = "Population size (number of wolf territories)") plot(Figure_2) ggsave("Figure2.jpg",Figure_2,width = 16,height = 7,dpi = 600, units = "cm") #reproductive age category - females# #data_le_fem = data, subset reproductive_age_category for females mod_le_fem <- glm(earlier_first_reproduction ~ PartExp + PopSize + I(PopSize^2) + F, data = data_le_fem, na.action = "na.fail",family=binomial) dred <- dredge(mod_le_fem, subset = dc(PopSize, I(PopSize^2))) summary(model.avg(dred)) subset(dred, delta <= 2.0) model_le_fem = model.avg(get.models(dred, subset=delta <= 2.0)) summary(model_le_fem) new_le_fem <- expand.grid(PartExp="no", PopSize=seq(0,max(data_le_fem$PopSize)), F=mean(data_le_fem$F)) predictl_le_fem <- predict(model_le_fem,new_le_fem,type="response",se.fit=T,full=TRUE) new_le_fem$Fit <- predictl_le_fem$fit new_le_fem$se <- predictl_le_fem$se.fit Figure_le_fem_nb_terr <- ggplot() + xlim(0,80) + ylim(-0.05, 1.05) + geom_jitter(data=data_le_fem, height=0.01, width=0.5,aes(x=PopSize, y=earlier_first_reproduction), size=0.5) + scale_color_manual("",values=c("#000000","#000000"), guide=F) + geom_ribbon(data = new_le_fem, aes(ymax=Fit+se,ymin=Fit-se, x=PopSize),alpha=0.5, fill=c("grey")) + geom_line(data = new_le_fem, aes(x=PopSize, y=Fit), colour=c("#000000"), size=0.5) + labs(y = "", x = "",size=12) + theme(axis.ticks = element_line(colour = "#000000", size = 0.2), axis.line = element_line(colour = "#000000", size = 0.2), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x = element_text(size=12), axis.title.y = element_text(size=12), axis.text.x = element_text(size=12), axis.text.y = element_text(size=12), plot.title =element_text(size=12)) + ggtitle("a)") Figure_le_fem_nb_terr #reproductive age category - males# #data_le_male = data, subset reproductive_age_category for males mod_le_male <- glm(earlier_first_reproduction ~ PartExp + PopSize + I(PopSize^2) + F, data = data_le_male, na.action = "na.fail",family=binomial) dred <- dredge(mod_le_male, subset = dc(PopSize, I(PopSize^2))) summary(model.avg(dred)) subset(dred, delta <= 2.0) model_le_male = model.avg(get.models(dred, subset=delta <= 2.0)) summary(model_le_male) new_le_male <- expand.grid(PartExp="no", PopSize=seq(0,max(data_le_male$PopSize)), F=mean(data_le_male$F)) predictl_le_male <- predict(model_le_male,new_le_male,type="response",se.fit=T,full=TRUE) new_le_male$Fit <- predictl_le_male$fit new_le_male$se <- predictl_le_male$se.fit Figure_le_male_nb_terr <- ggplot() + xlim(0,80) + ylim(-0.05, 1.05) + geom_jitter(data=data_le_male, height=0.01, width=0.5,aes(x=PopSize, y=earlier_first_reproduction), size=0.5) + scale_color_manual("",values=c("#000000","#000000"), guide=F) + geom_ribbon(data = new_le_male, aes(ymax=Fit+se,ymin=Fit-se, x=PopSize),alpha=0.5, fill=c("grey")) + geom_line(data = new_le_male, aes(x=PopSize, y=Fit), colour=c("#000000"), size=0.5) + labs(y = "", x = "",size=12) + theme(axis.ticks = element_line(colour = "#000000", size = 0.2), axis.line = element_line(colour = "#000000", size = 0.2), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x = element_text(size=12), axis.title.y = element_text(size=12), axis.text.x = element_text(size=12), axis.text.y = element_text(size=12), plot.title =element_text(size=12)) + ggtitle("b)") Figure_le_male_nb_terr #Figure 3# Figure_3 <- grid.arrange(Figure_le_fem_nb_terr, Figure_le_male_nb_terr, ncol = 2, left = "Probability to reproduce later", bottom = "Population size (number of wolf territories)") plot(Figure_3) ggsave("Figure_3.jpg",Figure_3,width = 16,height = 7,dpi = 600, units = "cm")