#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #2.1 fit lm of N swans ~ time (2017-2019) to get estimate of lambda (slope of N swans ~ time) #values are from counts <- data.frame(time=c(2017, 2018, 2019), AHY=c(4976.5, 5668, 6180.25), HY=c(3035, 1548, 2102.75)) # enter data counts$total <- counts$AHY+counts$HY # calc total swans (hatch-year + after-hatch-year birds) WETcounts <- data.frame(time=seq(2012, 2019,1), total=c(8598, 7473, 5528, 7186, 5898, 9531, 8185, 9971)) # note I've subtracted 1 year from each data point (e.g. data begins with 2012 instead of 2013) bc counts were conducted in Feb, but I consider them part of the breeding season that occured the prior year. This makes the data match with ours. # LM based on our aerial count data m1 <- lm(log(total) ~ time, data=counts) # fit lm m1.slope <- summary(m1)$coefficients[2] # extract slope for total ~ time m1.slope.se <- summary(m1)$coefficients[2,2] # extract se for slope # LM based on WET count data m2 <- lm(log(total) ~ time, data=WETcounts) m2.slope <- summary(m2)$coefficients[2] # extract slope for total ~ time m2.slope.se <- summary(m2)$coefficients[2,2] # extract se for slope # predictions + CIs for m1 and m2 CIs newx1 <- seq(min(counts$time), max(counts$time), length.out = 100) preds1 <- predict.lm(m1, newdata = data.frame(time=newx1), se.fit = T) newx2 <- seq(min(WETcounts$time), max(WETcounts$time), length.out = 100) preds2 <- predict.lm(m2, newdata = data.frame(time=newx2), se.fit = T) #plot m1 vs m2 #png("Plot - trend.png", res = 300, width = 4, height = 4, units = 'in') #par(mgp=c(2,0.75,0), mar=c(5,4,4,2)) plot(WETcounts$time, log(WETcounts$total), ylab="Log population size", xlab="Breeding season", ylim=c(8.2,9.5), col="red", pch=16, cex.axis=0.7,las=1,tck=-0.02) points(counts$time, log(counts$total), pch=16, col="blue") clip(2017,2019,0,10); abline(m1, col='blue', lwd=2) clip(2012,2019,0,10); abline(m2, col='red', lwd=2) lines(newx2, preds2$fit+1.96*preds2$se.fit, lty = 'dashed', col = 'red') lines(newx2, preds2$fit-1.96*preds2$se.fit, lty = 'dashed', col = 'red') lines(newx1, preds1$fit+1.96*preds1$se.fit, lty = 'dashed', col = 'blue') lines(newx1, preds1$fit-1.96*preds1$se.fit, lty = 'dashed', col = 'blue') legend(2012, 9.5,legend=c("Aerial counts", "WET counts"),col=c("blue","red"),box.lty=0,lty=1,cex=0.7) text(2012, 9.2, "slope = 0.035 (SE = 0.032)", col="red", cex=0.55, pos=4) text(2015.5, 8.4, "slope = 0.017 (SE = 0.070)", col="blue", cex=0.55, pos=4) #dev.off() #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #2.2 generate sampling distribution for lambda following normal distribution based on LM slope above # lambda dist from m1 (our data) m1.slope.dist <- rnorm(10000, m1.slope, m1.slope.se) m1.lambda.dist <- exp(m1.slope.dist) mean(m1.lambda.dist) lambda<-m1.lambda.dist # save another copy for using later in back-calculations quantile(lambda,probs=c(.025,.975)) # clean up rm(counts,WETcounts,preds1,preds2,m1.lambda.dist,m1.slope,m1.slope.dist,m1.slope.se,m2.slope,m2.slope.se,newx1,newx2) rm(m1,m2,m3,m4,m5,m6,fKM) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #2.3 generate sampling distributions for other parameters (except s.hy=hatch year survival, which we will back calculate) # annual survival after hatch year s.ahy<-rbinom(10000,17,14/17)/17 # based on n = 17 GPS collars with known fate from Te Waihora plot(density(s.ahy,adjust=3)) hist(s.ahy) mean(s.ahy) quantile(s.ahy,probs=c(.025,.975)) #se<-sqrt((48*mean(s.ahy)*(1-mean(s.ahy)))/length(s.ahy)) # for calculating SE of mean of n independent samples from Binomial(k,p) distribution, use sqrt(kpq/n) # abundance of after-hatch-year swans ahy<-data.frame(ahy=c(5395,5538,5756,5983)) # based on 2 counts by each of 2 independent observers in 2018 n.ahy<-rpois(10000,mean(ahy$ahy)) plot(density(n.ahy,adjust=2)) hist(n.ahy) mean(n.ahy) quantile(n.ahy,probs=c(.025,.975)) # number of nests nests<-data.frame(nests=c(840,871)) # based on 18 Sept 2018 census (116 nests counted outside photos area + 724 confirmed active nests (N=840), and 31 additional potentially active nests (N=871). n.nests<-rpois(10000,mean(nests$nests)) # using rpois plot(density(n.nests,adjust=2)) hist(n.nests) mean(n.nests) quantile(n.nests,probs=c(.025,.975)) # recruitment per unharvested nest (i.e. f.0 in the paper; fertility per unharvested nest) recruit = read.csv("newnestdata_for_backcalc_SHY.csv") # load data recruit2 <- subset(recruit, Camera.initial == "N") # remove nests with cameras, bc we found previously that cameras reduce recruitment AT <- subset(recruit2, Colony == "AT") #subset by Ataahua colony KA <- subset(recruit2, Colony == "KA") #subset by Kaituna colony KA <- subset(KA, Nest.ID != "KA16") # make sure KA16 is not included bc all eggs had unknown fate YA <- subset(recruit2, Colony == "YA") #subset by Yarrs colony # calc r mean for each colony with no harvest r.AT <- mean(AT$Recruitment.attempt1[AT$Eggs.removed=="0"]) # r for AT with no harvest r.KA <- mean(KA$Recruitment.attempt1[KA$Eggs.removed=="0"]) # r for KA with no harvest r.YA <- mean(YA$Recruitment.attempt1[YA$Eggs.removed=="0"]) # r for YA with no harvest # calc weights based on colony size # w.AT <- 130/647 # AT weight based on 130 nests divided by 647 (total for AT, KA, and YA) w.KA <- 393/647 # KA weight based on 393 nests divided by 647 (total for AT, KA, and YA) w.YA <- 124/647 # YA weight based on 124 nests divided by 647 (total for AT, KA, and YA) #wt.r <- ((w.AT*r.AT)+(w.KA*r.KA)+(w.YA*r.YA))/(w.AT+w.KA+w.YA) # calc weighted mean r #wt.sd <- sqrt(((w.AT*(r.AT-wt.r)^2)+(w.KA*(r.KA-wt.r)^2)+(w.YA*(r.YA-wt.r)^2))/((3-1)/3)*(w.AT+w.KA+w.YA)) # calc weighted sd # n nests in each colony sample n.AT <- 7 n.KA <- 7 n.YA <- 8 # sampling distribution for r, weighted by colony size, with uncertainty based on the sample size of nests per colony r <- w.AT*rpois(n=10000, lambda=r.AT*n.AT)/n.AT + w.KA*rpois(n=10000, lambda=r.KA*n.KA)/n.KA + w.YA*rpois(n=10000, lambda=r.YA*n.YA)/n.YA # plot r plot(density(r,adjust=3)) hist(r) mean(r) quantile(r,probs=c(.025,.975)) # survival of hatchling from day 0 to 40 s.0to40<-rbinom(10000,50,0.679)/50 # based on n = 50 vhf tags, KM model, with no events (confirmed mortalities) after day 40 plot(density(s.0to40,adjust=2)) hist(s.0to40) mean(s.0to40) quantile(s.0to40,probs=c(.025,.975)) #define sex ratio sex.ra<-0.5 # clean up rm(AT,KA,YA,r.AT,r.KA,r.YA,w.AT,w.KA,w.YA,n.AT,n.KA,n.YA,nests,recruit,recruit2) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #2.4 plots of the original sampling distributions used for back-calculating s.hy library(ggplot2) # annual survival after hatch year df.sahy <- as.data.frame(s.ahy) # calculate age-distribution-weighted mean annual survival after hatch year for M. Williams data to use in plot # Post-Wahine age structure from Williams 1979 Prop.Sj1 <- 0.193 # proportion of AHY swans that are age 1-2 (in early 1970s post-wahine; Williams 1979) Prop.Sj2 <- 0.079 # proportion of AHY swans that are age 2-3 (in early 1970s post-wahine; Williams 1979) Prop.Sj3 <- 0.117 # proportion of AHY swans that are age 3-4 (in early 1970s post-wahine; Williams 1979) Prop.A <- 0.611 # proportion of AHY swans that are age 4+ (in early 1970s post-wahine; Williams 1979) # age-specific annual survival from 1956-74 Sj1 <- 0.661 # survival from age 1-2 (from Williams 1979) Sj2 <- 0.784 # survival from age 2-3 (from Williams 1979) Sj3 <- 0.764 # survival from age 3-4 (from Williams 1979) Sa <- 0.836 # annual survival for ages 4+ (from Williams 1981) # age-structure-weighted mean annual survival based on Williams 1979, 1981 Age.wt.S.ahy <- (Sj1*Prop.Sj1)+(Sj2*Prop.Sj2)+(Sj3*Prop.Sj3)+(Sa*Prop.A) # 0.789693 # histogram with density plot p1.sahy <- ggplot(df.sahy, aes(x=s.ahy)) + geom_histogram(aes(y=..density..), colour="black", fill="white", binwidth = 0.05, boundary=1, size=0.6)+ geom_density(alpha=.4, adjust=3, color="black", fill="#00798c", size=0.6) + geom_vline(aes(xintercept=mean(s.ahy)), color="black",linetype="dashed", size=0.6)+ xlab("Annual survival after hatch year") + ylab("Density") + theme_classic() + theme(aspect.ratio=1) + annotate("segment", x = 0.84, xend = 0.84, y = 1, yend = 0, arrow=arrow(angle=15, type = "closed"), size=0.7) + # from Barker and Buchanan 1993 annotate("segment", x = 0.79, xend = 0.79, y = 1, yend = 0, arrow=arrow(angle=15, type = "closed"), size=0.7) + # from Williams 1979 geom_text(x=0.84, y=1.3, label="b", size=5) + geom_text(x=0.79, y=1.3, label="a", size=5) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), axis.title.x = element_text(face = "plain", colour = "black", size = 18), axis.title.y = element_text(face = "plain", colour = "black", size = 18), axis.ticks.x = element_line(color="black", size = 1.2), axis.ticks.y = element_line(color= "black", size = 1.2), axis.ticks.length = unit(0.2, "cm"), axis.text.x = element_text(family = "sans", face = "plain", colour = "black", size = 18, margin = margin(8,0,10,0)), axis.text.y = element_text(family = "sans", face = "plain", colour = "black", size = 18, margin = margin(0,8,0,10)), axis.line.x = element_line(color="black", size = 1), axis.line.y = element_line(color="black", size = 1)) p1.sahy #ggsave(p1.sahy, file="hist_sahy_originaldata.png", dpi = 300) # number of nests df.nests <- as.data.frame(n.nests) # histogram with density plot p1.nests <- ggplot(df.nests, aes(x=n.nests)) + geom_histogram(aes(y=..density..), colour="black", fill="white", binwidth = 25, boundary=800, size=0.6)+ geom_density(alpha=.4, adjust=3, color="black", fill="#00798c", size=0.6) + geom_vline(aes(xintercept=mean(n.nests)), color="black",linetype="dashed", size=0.6)+ xlab("Nest abundance") + ylab("Density") + scale_x_continuous(breaks=c(seq(750,1000,50)), limits = c(750,1010)) + theme_classic() + theme(aspect.ratio=1) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), axis.title.x = element_text(face = "plain", colour = "black", size = 18), axis.title.y = element_text(face = "plain", colour = "black", size = 18), axis.ticks.x = element_line(color="black", size = 1.2), axis.ticks.y = element_line(color= "black", size = 1.2), axis.ticks.length = unit(0.2, "cm"), axis.text.x = element_text(family = "sans", face = "plain", colour = "black", size = 18, margin = margin(8,0,10,0)), axis.text.y = element_text(family = "sans", face = "plain", colour = "black", size = 18, margin = margin(0,8,0,10)), axis.line.x = element_line(color="black", size = 1), axis.line.y = element_line(color="black", size = 1)) p1.nests #ggsave(p1.nests, file="hist_nests_originaldata.png", dpi = 300) # number of AHY swans df.nahy <- as.data.frame(n.ahy) # histogram with density plot p1.nahy <- ggplot(df.nahy, aes(x=n.ahy)) + geom_histogram(aes(y=..density..), colour="black", fill="white", boundary=5600, binwidth = 50, size=0.6)+ geom_density(alpha=.4, adjust=3, color="black", fill="#00798c", size=0.6) + geom_vline(aes(xintercept=mean(n.ahy)), color="black",linetype="dashed", size=0.6)+ xlab("After-hatch-year abundance") + ylab("Density") + scale_x_continuous(breaks=c(seq(5400,6000,200)), limits = c(5350,6050)) + theme_classic() + theme(aspect.ratio=1) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), axis.title.x = element_text(face = "plain", colour = "black", size = 18), axis.title.y = element_text(face = "plain", colour = "black", size = 18), axis.ticks.x = element_line(color="black", size = 1.2), axis.ticks.y = element_line(color= "black", size = 1.2), axis.ticks.length = unit(0.2, "cm"), axis.text.x = element_text(family = "sans", face = "plain", colour = "black", size = 18, margin = margin(8,0,10,0)), axis.text.y = element_text(family = "sans", face = "plain", colour = "black", size = 18, margin = margin(0,8,0,10)), axis.line.x = element_line(color="black", size = 1), axis.line.y = element_line(color="black", size = 1)) p1.nahy #ggsave(p1.nahy, file="hist_nahy_originaldata.png", dpi = 300) # recruitment df.r <- as.data.frame(r) # histogram with density plot p1.r <- ggplot(df.r, aes(x=r)) + geom_histogram(aes(y=..density..), colour="black", fill="white", boundary=4.0, binwidth = 0.4, size=0.6)+ geom_density(alpha=.4, adjust=3, color="black", fill="#00798c", size=0.6) + geom_vline(aes(xintercept=mean(r)), color="black",linetype="dashed", size=0.6)+ xlab("Fertility per unharvested nest") + ylab("Density") + #scale_x_continuous(breaks=c(seq(3,5,0.5)), limits = c(3,5)) + theme_classic() + theme(aspect.ratio=1) + annotate("segment", x = 3.16, xend = 3.16, y = .1, yend = 0, arrow=arrow(angle=15, type = "closed"), size=0.7) + annotate("segment", x = 3.87, xend = 3.87, y = .1, yend = 0, arrow=arrow(angle=15, type = "closed"), size=0.7) + annotate("segment", x = 3.93, xend = 3.93, y = .1, yend = 0, arrow=arrow(angle=15, type = "closed"), size=0.7) + geom_text(x=3.16, y=.14, label="a", size=5) + geom_text(x=3.8, y=.14, label="b", size=5) + geom_text(x=4, y=.14, label="c", size=5) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), axis.title.x = element_text(face = "plain", colour = "black", size = 18), axis.title.y = element_text(face = "plain", colour = "black", size = 18), axis.ticks.x = element_line(color="black", size = 1.2), axis.ticks.y = element_line(color= "black", size = 1.2), axis.ticks.length = unit(0.2, "cm"), axis.text.x = element_text(family = "sans", face = "plain", colour = "black", size = 18, margin = margin(8,0,10,0)), axis.text.y = element_text(family = "sans", face = "plain", colour = "black", size = 18, margin = margin(0,8,0,10)), axis.line.x = element_line(color="black", size = 1), axis.line.y = element_line(color="black", size = 1)) p1.r #ggsave(p1.r, file="hist_r_originaldata.png", dpi = 300) #lambda df.lam <- as.data.frame(lambda) # histogram with density plot p1.lam <- ggplot(df.lam, aes(x=lambda)) + geom_histogram(aes(y=..density..), colour="black", fill="white", boundary=0, binwidth = 0.05, size=0.6)+ geom_density(alpha=.4, adjust=3, color="black", fill="#00798c", size=0.6) + geom_vline(aes(xintercept=mean(lambda)), color="black",linetype="dashed", size=0.6)+ xlab("Population growth rate") + ylab("Density") + theme_classic() + theme(aspect.ratio=1) + #geom_text(x=1.2, y=5, label="mean = 1.019", size=5) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), axis.title.x = element_text(face = "plain", colour = "black", size = 18), axis.title.y = element_text(face = "plain", colour = "black", size = 18), axis.ticks.x = element_line(color="black", size = 1.2), axis.ticks.y = element_line(color= "black", size = 1.2), axis.ticks.length = unit(0.2, "cm"), axis.text.x = element_text(family = "sans", face = "plain", colour = "black", size = 18, margin = margin(8,0,10,0)), axis.text.y = element_text(family = "sans", face = "plain", colour = "black", size = 18, margin = margin(0,8,0,10)), axis.line.x = element_line(color="black", size = 1), axis.line.y = element_line(color="black", size = 1)) p1.lam #ggsave(p1.lam, file="hist_lambda_originaldata_aerialcounts.png", dpi = 300) #survival of VHF cygnets over first 40 days df.s0to40 <- as.data.frame(s.0to40) # histogram with density plot p1.s0to40 <- ggplot(df.lam, aes(x=s.0to40)) + geom_histogram(aes(y=..density..), colour="black", fill="white", boundary=0, binwidth = 0.04, size=0.6)+ geom_density(alpha=.4, adjust=3, color="black", fill="#00798c", size=0.6) + geom_vline(aes(xintercept=mean(s.0to40)), color="black",linetype="dashed", size=0.6)+ xlab("40-day hatchling survival probability") + ylab("Density") + theme_classic() + theme(aspect.ratio=1) + geom_text(x=1.2, y=5, label="mean = 1.019", size=5) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), axis.title.x = element_text(face = "plain", colour = "black", size = 16), axis.title.y = element_text(face = "plain", colour = "black", size = 16), axis.ticks.x = element_line(color="black", size = 1.2), axis.ticks.y = element_line(color= "black", size = 1.2), axis.ticks.length = unit(0.2, "cm"), axis.text.x = element_text(family = "sans", face = "plain", colour = "black", size = 18, margin = margin(8,0,10,0)), axis.text.y = element_text(family = "sans", face = "plain", colour = "black", size = 18, margin = margin(0,8,0,10)), axis.line.x = element_line(color="black", size = 1), axis.line.y = element_line(color="black", size = 1)) p1.s0to40 #ggsave(p1.s0to40, file="hist_s0to40.png", dpi = 300) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #2.5 back-calculate hatch year survival # 2.5.1 randomly draw from lambda, n.ahy, n.nests, r, s.ahy, s.0to40 # 2.5.2 solve s.hy = (lambda^2-lambda*s.ahy)/((2*n.nests/n.ahy)*r*sex.ra*s.ahy) # 2.5.3 make sure s.hy is < than s.hy0to40; if not, re-draw everything # 2.5.4 make sure s.hy is >= 0; if not, re-draw everything # 2.5.5 make sure s.40to356, i.e. (s.hy/s.0to40))^(1/325), is >= (s.0to40)^(1/40); if not, re-draw everything # 2.5.6 put all values in table or df # 2.5.7 repeat until you have 10000 samples df <- data.frame(lambda=lambda, n.ahy=n.ahy, n.nests=n.nests, r=r, s.ahy=s.ahy, s.0to40=s.0to40) blah <- matrix(ncol=ncol(df), nrow=1) blah <- as.data.frame(blah) colnames(blah) <- c("lambda", "n.ahy", "n.nests", "r", "s.ahy", "s.0to40") output <- rep(list(),length=10000) # create list for output == 10000 (desired length) output_s.hy_unconditiional <- rep(list(),length=100000) # create list for all s.hy values created, even when conditions are not met (use length longer than we need) output_sdaily_unconditional <- rep(list(),length=100000) # create list for all sdaily values created when first two conditions are met but third conditions are not met for (j in 1 : 100000) { # tell the outside loop to run more times than we need for (i in 1 : ncol(df)){ # i varies from 1 to n of col in df blah[,i] <- sample(df[,i],1) # randomly sample 1 value from each col in df; we need a new blah for each run of the outside loop } s.hy <- ((blah$lambda^2)-(blah$lambda*blah$s.ahy))/(((2*blah$n.nests)/blah$n.ahy)*blah$r*sex.ra*blah$s.ahy) # use randomly drawn values to calculate s.hy blah2 <- c(blah, s.hy) # add result of equation (s.hy) into object blah2 with the values used to calculate it blah2 <- as.data.frame(blah2) colnames(blah2)[7] <- "s.hy" # naming s.hy column output_s.hy_unconditiional[[j]] <- s.hy # saving the value of s.hy to the list if (s.hy <= blah2$s.0to40) { # first condition - s.hy must be <= to randomly drawn s0to40; only run next part of code if condition is met if (s.hy >= 0) { # second condition - if first condition is met, check that s.hy is >= 0 sdaily40to365 <- (s.hy/blah2$s.0to40)^(1/325) # defining equation for daily survival from day 40 to 365 sdaily0to40 <- (blah2$s.0to40)^(1/40) # defining equation for daily survival from day 0 to 40 sdaily <- c(sdaily0to40,sdaily40to365,s.hy) # concatenate sdaily values and s.hy sdaily <- t(as.data.frame(sdaily)) # save sdaily values and s.hy to df colnames(sdaily) <- c("sdaily0to40","sdaily40to365","s.hy") # name columns in df output_sdaily_unconditional[[j]] <- sdaily # save sdaily df to list if (sdaily40to365 >= sdaily0to40) { # third condition (this might need to be > rather than >= ?) parameters <- as.data.frame(c(blah2, sdaily0to40, sdaily40to365)) # combining blah2 with the two new objects (daily survival) colnames(parameters)[c(8,9)] <- c("sdaily0to40","sdaily40to365") output[[j]] <- parameters output_final <- output[-which(lapply(output, is.null)==T)] # remove the null elements in the output list, so the condition below (length==10000) isn't met based with null values if(length(output_final)==10000){ # if output has 10000 elements, stop running the outside loop break } } } } } # code for post-processing, to determine how many runs failed, and identify those runs: # create lists of unconditional data; the first list includes elements for which conditions 1-2 were not met (and therefore includes elements for which condition 3 was not met); the second list includes elements for which conditions 1-2 were met, and condition 3 were not met #output_s.hy_unconditiional2 <- output_s.hy_unconditiional[-which(lapply(output_s.hy_unconditiional, is.null)==T)] # remove all remaining null elements from list after finishing #output_sdaily_unconditional2 <- output_sdaily_unconditional[-which(lapply(output_sdaily_unconditional, is.null)==T)] # remove all remaining null elements from list after finishing # function for identifying the difference between s40to365 and s0to40 from columns 2 and 1, respectively, in sdaily list defined above #function_mark <- function(x) (x[,2]-x[,1]) #dif <- lapply(output_sdaily_unconditional2, function_mark) # use function for all elements of list #positions <- which(dif<=0) # positions in which the subtraction is < 0 #output_sdaily_unconditional2[positions] # sublist showing elements where condition 1 was not met # unlist the new distributions for s.hy and other parameters to use in the final model # hatch year survival s.hy <- unlist(lapply(output_final, '[', c("s.hy")), use.names = F) hist(s.hy) mean(s.hy) sd(s.hy) quantile(s.hy,probs=c(.025,.975)) # after hatch year survival s.ahy <- unlist(lapply(output_final, '[', c("s.ahy")), use.names = F) hist(s.ahy) plot(density(s.ahy,adjust=5)) mean(s.ahy) sd(s.ahy) quantile(s.ahy,probs=c(.025,.975)) # recruitment, i.e. f.0 r <- unlist(lapply(output_final, '[', c("r")), use.names = F) hist(r) plot(density(r,adjust=5)) mean(r) sd(r) quantile(r,probs=c(.025,.975)) # number of AHY swans n.ahy <- unlist(lapply(output_final, '[', c("n.ahy")), use.names = F) hist(n.ahy) plot(density(n.ahy,adjust=5)) mean(n.ahy) sd(n.ahy) quantile(n.ahy,probs=c(.025,.975)) # number of nests n.nests <- unlist(lapply(output_final, '[', c("n.nests")), use.names = F) hist(n.nests) plot(density(n.nests,adjust=5)) mean(n.nests) sd(n.nests) quantile(n.nests,probs=c(.025,.975)) # lambda lam <- unlist(lapply(output_final, '[', c("lambda")), use.names = F) mean(lam) median(lam) sd(lam) quantile(lam,probs=c(.025,.975)) # hachling survival s.hatchling <- unlist(lapply(output_final, '[', c("s.0to40")), use.names = F) mean(s.hatchling) median(s.hatchling) sd(s.hatchling) quantile(s.hatchling,probs=c(.025,.975)) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #2.6 plots of distributions resulting/retained in final model (these distributions have been adjusted) # s.hy s.hy.final <- as.data.frame(s.hy) # histogram with density plot p.s.hy.final <- ggplot(s.hy.final, aes(x=s.hy)) + #geom_histogram(aes(y=..density..), colour="black", fill="white", boundary=0, binwidth = 0.05, size=0.6)+ geom_density(alpha=.4, adjust=3, color="black", fill="#00798c", size=0.8) + geom_vline(aes(xintercept=mean(s.hy)), color="black",linetype="dashed", size=1)+ xlab("Hatch-year survival probability") + ylab("Density") + theme_classic() + theme(aspect.ratio=1) + scale_x_continuous(breaks = c(0,0.2,0.4,0.6,0.8,1), limits = c(0,1))+ #annotate("segment", x = 0.3, xend = 0.3, y = 0.4, yend = 0, arrow=arrow(angle=15, type = "closed"), size=0.7) + #geom_text(x=0.3, y=0.5, label="A", size=5) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), #panel.border = element_blank(), panel.border = element_rect(colour = "black", fill=NA, size=1), axis.title.x = element_text(face = "plain", colour = "black", size = 18), axis.title.y = element_text(face = "plain", colour = "black", size = 18), axis.ticks.x = element_line(color="black", size = 1), axis.ticks.y = element_line(color= "black", size = 1), axis.ticks.length = unit(0.2, "cm"), axis.text.x = element_text(family = "sans", face = "plain", colour = "black", size = 21, margin = margin(8,0,10,0)), axis.text.y = element_text(family = "sans", face = "plain", colour = "black", size = 21, margin = margin(0,8,0,10)), axis.line.x = element_line(color="black", size = 0.8), axis.line.y = element_line(color="black", size = 0.8))+ coord_flip() p.s.hy.final #ggsave(p.s.hy.final, file="hist_shy_final_forplottingwithtrend.png", dpi = 300) # sahy s.ahy.final <- unlist(lapply(output_final, '[', c("s.ahy")), use.names = F) s.ahy.final <- as.data.frame(s.ahy.final) # histogram with density plot p.s.ahy.final <- ggplot(s.ahy.final, aes(x=s.ahy)) + geom_histogram(aes(y=..density..), colour="black", fill="white", boundary=0, binwidth = 0.05, size=0.6)+ geom_density(alpha=.4, adjust=3, color="black", fill="#00798c", size=0.6) + geom_vline(aes(xintercept=mean(s.ahy.final)), color="black",linetype="dashed", size=0.6)+ xlab("Annual survival after hatch year") + ylab("Density") + scale_x_continuous(breaks=c(seq(0.6,1,0.2)), limits = c(0.5,1)) + theme_classic() + theme(aspect.ratio=1) + annotate("segment", x = 0.84, xend = 0.84, y = 1, yend = 0, arrow=arrow(angle=15, type = "closed"), size=0.7) + # from Barker and Buchanan 1993 annotate("segment", x = 0.79, xend = 0.79, y = 1, yend = 0, arrow=arrow(angle=15, type = "closed"), size=0.7) + # from Williams 1979 geom_text(x=0.84, y=1.3, label="b", size=5) + geom_text(x=0.79, y=1.3, label="a", size=5) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), axis.title.x = element_text(face = "plain", colour = "black", size = 18), axis.title.y = element_text(face = "plain", colour = "black", size = 18), axis.ticks.x = element_line(color="black", size = 1.2), axis.ticks.y = element_line(color= "black", size = 1.2), axis.ticks.length = unit(0.2, "cm"), axis.text.x = element_text(family = "sans", face = "plain", colour = "black", size = 18, margin = margin(8,0,10,0)), axis.text.y = element_text(family = "sans", face = "plain", colour = "black", size = 18, margin = margin(0,8,0,10)), axis.line.x = element_line(color="black", size = 1), axis.line.y = element_line(color="black", size = 1)) p.s.ahy.final #ggsave(p.s.ahy.final, file="hist_sahy_final.png", dpi = 300) # recruitment r.final <- unlist(lapply(output_final, '[', c("r")), use.names = F) r.final <- as.data.frame(r.final) # histogram with density plot p.r.final <- ggplot(r.final, aes(x=r)) + geom_histogram(aes(y=..density..), colour="black", fill="white", boundary=0, binwidth = 0.3, size=0.6)+ geom_density(alpha=.4, adjust=3, color="black", fill="#00798c", size=0.6) + geom_vline(aes(xintercept=mean(r.final)), color="black",linetype="dashed", size=0.6)+ xlab("No. eggs hatched per nest") + ylab("Density") + #scale_x_continuous(breaks=c(seq(3,5,0.5)), limits = c(3,5)) + theme_classic() + theme(aspect.ratio=1) + annotate("segment", x = 3.16, xend = 3.16, y = .1, yend = 0, arrow=arrow(angle=15, type = "closed"), size=0.7) + annotate("segment", x = 3.87, xend = 3.87, y = .1, yend = 0, arrow=arrow(angle=15, type = "closed"), size=0.7) + annotate("segment", x = 3.93, xend = 3.93, y = .1, yend = 0, arrow=arrow(angle=15, type = "closed"), size=0.7) + geom_text(x=3.16, y=.12, label="a", size=5) + geom_text(x=3.8, y=.12, label="b", size=5) + geom_text(x=4, y=.12, label="c", size=5) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), axis.title.x = element_text(face = "plain", colour = "black", size = 18), axis.title.y = element_text(face = "plain", colour = "black", size = 18), axis.ticks.x = element_line(color="black", size = 1.2), axis.ticks.y = element_line(color= "black", size = 1.2), axis.ticks.length = unit(0.2, "cm"), axis.text.x = element_text(family = "sans", face = "plain", colour = "black", size = 18, margin = margin(8,0,10,0)), axis.text.y = element_text(family = "sans", face = "plain", colour = "black", size = 18, margin = margin(0,8,0,10)), axis.line.x = element_line(color="black", size = 1), axis.line.y = element_line(color="black", size = 1)) p.r.final #ggsave(p.r.final, file="hist_r_final.png", dpi = 300) # N nests nests.final <- unlist(lapply(output_final, '[', c("n.nests")), use.names = F) nests.final <- as.data.frame(nests.final) # histogram with density plot p.nests.final <- ggplot(nests.final, aes(x=n.nests)) + geom_histogram(aes(y=..density..), colour="black", fill="white", binwidth = 20, boundary=800, size=0.6)+ geom_density(alpha=.4, adjust=3, color="black", fill="#00798c", size=0.6) + geom_vline(aes(xintercept=mean(n.nests)), color="black",linetype="dashed", size=0.6)+ xlab("Nest abundance") + ylab("Density") + scale_x_continuous(breaks=c(seq(750,950,50)), limits = c(750,1010)) + theme_classic() + theme(aspect.ratio=1) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), axis.title.x = element_text(face = "plain", colour = "black", size = 18), axis.title.y = element_text(face = "plain", colour = "black", size = 18), axis.ticks.x = element_line(color="black", size = 1.2), axis.ticks.y = element_line(color= "black", size = 1.2), axis.ticks.length = unit(0.2, "cm"), axis.text.x = element_text(family = "sans", face = "plain", colour = "black", size = 18, margin = margin(8,0,10,0)), axis.text.y = element_text(family = "sans", face = "plain", colour = "black", size = 18, margin = margin(0,8,0,10)), axis.line.x = element_line(color="black", size = 1), axis.line.y = element_line(color="black", size = 1)) p.nests.final #ggsave(p.nests.final, file="hist_nests_final.png", dpi = 300) # N ahy ahy.final <- unlist(lapply(output_final, '[', c("n.ahy")), use.names = F) ahy.final <- as.data.frame(ahy.final) # histogram with density plot p.nahy.final <- ggplot(ahy.final, aes(x=n.ahy)) + geom_histogram(aes(y=..density..), colour="black", fill="white", boundary=5600, binwidth = 50, size=0.6)+ geom_density(alpha=.4, adjust=3, color="black", fill="#00798c", size=0.6) + geom_vline(aes(xintercept=mean(n.ahy)), color="black",linetype="dashed", size=0.6)+ xlab("After-hatch-year abundance") + ylab("Density") + scale_x_continuous(breaks=c(seq(5400,6000,200)), limits = c(5350,6050)) + theme_classic() + theme(aspect.ratio=1) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), axis.title.x = element_text(face = "plain", colour = "black", size = 18), axis.title.y = element_text(face = "plain", colour = "black", size = 18), axis.ticks.x = element_line(color="black", size = 1.2), axis.ticks.y = element_line(color= "black", size = 1.2), axis.ticks.length = unit(0.2, "cm"), axis.text.x = element_text(family = "sans", face = "plain", colour = "black", size = 18, margin = margin(8,0,10,0)), axis.text.y = element_text(family = "sans", face = "plain", colour = "black", size = 18, margin = margin(0,8,0,10)), axis.line.x = element_line(color="black", size = 1), axis.line.y = element_line(color="black", size = 1)) p.nahy.final #ggsave(p.nahy.final, file="hist_nahy_final.png", dpi = 300)