In nwfsc-cb/srkw-status: Package for processing and exploring killer whale demographic data
knitr::opts_chunk$set(echo = FALSE, warning=FALSE, message=FALSE,
comment = "#>"
)
library(kwdemog)
library(ggplot2)
library(knitr)
library(dplyr)
library(mgcv)
library(rstanarm)
library(viridis)
library(ggridges)
data(orca)
year.end = 2024
refit_models = FALSE
report_start = year.end-5
# format whale names just for consistency
format_names = function(x) {
for(i in 1:2) {
indx = which(substr(x,2,2) == "0")
x[indx] = paste0(substr(x[indx],1,1), substr(x[indx],3,nchar(x[indx])))
}
return(x)
}
orca$animal = format_names(orca$animal)
orca$pod = format_names(orca$pod)
orca$matriline = format_names(orca$matriline)
orca$mom = format_names(orca$mom)
orca$dad = format_names(orca$dad)
report_dir = "projections/"
Abstract
This document is intended to provide and overview and status update of the Southern Resident population following the most recent summer census in r year.end. Many of the analyses and figures presented below are described in more detail elsewhere, particularly the most recent Status Review [@nmfs2016] and recent technical reports following the 2011 - 2012 bilateral workshops [@hilborn2012; @ward2013]. These analyses also formed the backbone of the Pacific Fishery Management Council's working group on salmon fisheries impacts to killer whales [@pfmc2020; @pfmc2020b]. In particular, these analyses are meant to describe recent changes in population size and age structure, change in demographic rates over time, and an updated projections of population viability.
Recent births and deaths
Since r report_start, the following animals have been born in the SRKW population and lived:
whaleData = orca
additions = whaleData[which(whaleData$pod%in%c("J1","K1","L1") &
whaleData$birth > report_start & is.na(whaleData$death)),]
additions = additions[,c("animal","birth","sexF1M2")]
additions$sexF1M2[additions$sexF1M2=="0"]="Unk"
additions$sexF1M2[additions$sexF1M2=="1"]="F"
additions$sexF1M2[additions$sexF1M2=="2"]="M"
names(additions)[3] = "sex"
additions = dplyr::arrange(additions, birth, animal)
additions = as.matrix(additions)
row.names(additions) = NULL
kable(as.matrix(additions))
The following whales have died since r report_start. Their deaths are also included in this table,
whaleData = orca
losses = dplyr::filter(whaleData, pod%in%c("J1","K1","L1"), death > report_start)
# subset columns
losses = losses[,c("animal","death","sexF1M2","birth")]
losses$age_death = losses$death-losses$birth
losses$sexF1M2[losses$sexF1M2=="0"]="Unk"
losses$sexF1M2[losses$sexF1M2=="1"]="F"
losses$sexF1M2[losses$sexF1M2=="2"]="M"
names(losses)[3] = "sex"
losses = as.matrix(losses[,c("animal","death","sex","age_death")])
row.names(losses) = NULL
# order
losses = arrange(as.data.frame(losses), death, animal) %>% as.matrix()
kable(losses)
Sex ratio at birth
During the 2011-2012 bilateral workshops [@hilborn2012; @ward2013], a comparison between NRKW and SRKW sex ratios at birth was presented, with calves being approximately 55\% female in the NRKW population and 45\% female in the SRKW population. This difference was assumed to be due to chance, and there was no evidence for a significant trend. As the proportion of males in the SRKW population has increased over time, it is worth re-examining the evidence supporting any trend.
sexr = dplyr::filter(whaleData, birth >= (year.end-3)) %>%
dplyr::summarize(males = length(which(sexF1M2 == 2)),
females = length(which(sexF1M2==1)))
# do a simple bayesian logistic regression
indx = which(orca$sexF1M2 != 0 & orca$birth >= 1976 & orca$birth <= year.end)
Y = orca$sexF1M2[indx] - 1
Xt = orca$birth[indx]
df = data.frame("year"=Xt,"sex"=Y)
g1 = ggplot(df, aes(year, sex)) + geom_hline(aes(yintercept=0.5),col="red",alpha=0.3) + geom_point(col="blue",alpha=0.3, position=position_dodge(width = 1)) + geom_smooth(method = "glm", method.args = list(family = "binomial")) + xlab("Year") + ylab("Sex (0 = F, 1 = M)") + theme_bw()
# do a quick stepwise model selection using AIC
#g.null = glm(Y ~ 1, family = "binomial")
g.t = stan_glm(Y ~ Xt, family = "binomial")
draws = as.matrix(g.t)
ppos = 100 * round(length(which(draws[,2] > 0)) / nrow(draws), 3)
df = data.frame("x" = draws[,2])
g2 = ggplot(data=df, aes(x=x)) +
geom_histogram(fill = "dodgerblue",alpha=0.6) +
xlab("Annual trend in male births") + ylab("Frequency") +
geom_vline(aes(xintercept=0), col="grey30") +
geom_vline(aes(xintercept=quantile(df$x,0.025)), col="red",alpha=0.3) +
geom_vline(aes(xintercept=quantile(df$x,0.975)), col="red",alpha=0.3) +
theme_bw()
To evaluate support for a trend, we fit Bayesian logistic regression models (GLMs with a binomial family and logit link function), to SRKW birth data over the period 1977-present. In recent years, since r year.end-3, there have been r sexr$males male births and r sexr$females female births. This analysis highlights that the probability of a positive trend is approximately r ppos\% (Fig. \ref{fig:srb}).
```r"}
gridExtra::grid.arrange(g1, g2, ncol=1)
# Changing Population Structure
One of the objectives of the recovery goals [@nmfs2008; @nmfs2011] was an age and sex distribution that is more similar to the Northern Resident population (at the time of ESA listing). The Southern Resident population has undergone a number of shifts in age and sex, and because the population is so small, the age or sex composition are more sensitive to individual births and deaths. Previous status reviews based these targets from 1973-1996 [@olesiuk2005]:
```r
m = matrix("", 4, 2)
m[,1] = c("Juveniles", "Reproductive females", "Post-reproductive females", "Adult males")
m[,2] = c("47 %", "24 %", "11 %", "18 %")
colnames(m) = c("Stage","")
kable(m)
We can re-evaluate these targets, both for the Northern and Southern Resident populations, using the most recent years of data available (2018 for Northern Resident, r year.end for Southern Resident). There are some small differences between life stages used in Olesiuk et al. 2005, and more recent work [@ward2013]. First, Olesiuk assumed animals to not be mature to 15.5 (wheras Ward et al. 2013 assumed females to be mature at age 10). Second, Olesiuk et al. (2005) defined post-reproductive animals to be animals who hadn't given birth in 10 years (Ward et al. 2013 used a cutoff of 42 years). For these calculations, we define age at maturity to be 10 years, and 42+ as the age of reproductive senescence.
whaleData = kwdemog::expand(orca, start_year = 1976,
current_year = year.end)
# Sexes need to be randomly filled in
whaleData$sexF1M2[which(whaleData$sexF1M2==0)] = sample(c(1,2), size=length(which(whaleData$sexF1M2==0)), replace=T)
whaleData$age1979 = 1979 - whaleData$birth
whaleData$age2010 = 2010 - whaleData$birth
whaleData$ageCurrent = year.end - whaleData$birth
NR = dplyr::filter(whaleData, year == 1979, pod %in% c("J1","K1","L1") == FALSE,
alive==1, age1979 > 0) %>%
dplyr::filter(death > 1979 | is.na(death))
NR_summary_1979 = dplyr::summarize(NR,
x1 = length(which(age1979 < 10)),
x2 = length(which(age1979 >= 10 & sexF1M2==2)),
x3 = length(which(age1979 >= 10 & age1979 <= 42 & sexF1M2==1)),
x4 = length(which(age1979 > 42 & sexF1M2==1))) %>% as.matrix()
NR = dplyr::filter(whaleData, year == 2018, pod %in% c("J1","K1","L1") == FALSE,
alive==1, ageCurrent > 0) %>%
dplyr::filter(death > 2010 | is.na(death))
NR_summary = dplyr::summarize(NR,
x1 = length(which(ageCurrent < 10)),
x2 = length(which(ageCurrent >= 10 & sexF1M2==2)),
x3 = length(which(ageCurrent >= 10 & age2010 <= 42 & sexF1M2==1)),
x4 = length(which(ageCurrent > 42 & sexF1M2==1))) %>% as.matrix()
SR = dplyr::filter(whaleData, year == 1979, pod %in% c("J1","K1","L1"),
alive==1, age1979 > 0) %>%
dplyr::filter(death > 1979 | is.na(death))
SR_summary_1979 = dplyr::summarize(SR,
x1 = length(which(age1979 < 10)),
x2 = length(which(age1979 >= 10 & sexF1M2==2)),
x3 = length(which(age1979 >= 10 & age1979 <= 42 & sexF1M2==1)),
x4 = length(which(age1979 > 42 & sexF1M2==1))) %>% as.matrix()
SR = dplyr::filter(whaleData, year == year.end, pod %in% c("J1","K1","L1"),
alive==1, ageCurrent > 0) %>%
dplyr::filter(death > year.end | is.na(death))
SR_summary = dplyr::summarize(SR,
x1 = length(which(ageCurrent < 10)),
x2 = length(which(ageCurrent >= 10 & sexF1M2==2)),
x3 = length(which(ageCurrent >= 10 & ageCurrent <= 42 & sexF1M2==1)),
x4 = length(which(ageCurrent > 42 & sexF1M2==1))) %>% as.matrix()
m = matrix(NA, 4, 4)
m[,1] = 100*round(c(SR_summary_1979/sum(SR_summary_1979)), 2)
m[,2] = 100*round(c(SR_summary/sum(SR_summary)), 2)
m[,3] = 100*round(c(NR_summary_1979/sum(NR_summary_1979)), 2)
m[,4] = 100*round(c(NR_summary/sum(NR_summary)), 2)
for(i in 1:nrow(m)) {
for(j in 1:ncol(m)) {
m[i,j] = paste0(m[i,j], " %")
}
}
rownames(m) = c("Juveniles (< 10)", "Adult males (10+)", "Adult females (10-42)", "Post-reproductive females (42+)")
colnames(m) = c("SRKW 1979", paste0("SRKW ",year.end), "NRKW 1979", "NRKW 2018")
kable(m)
Reproductive females
The number of reproductive aged females was at its lowest point in the late 1970s, in part because of the prior harvesting that occurred into the early 1970s (Fig. \ref{fig:ts-repro-females}). Though the overall number of reproductive females has been fluctuating between 25-35 for most of the last 40 years, there have been contrasting changes by pod, with declines in L pod females and increases in J pod (Fig. \ref{fig:ts-repro-females}). At the start of the survey in 1976, the distribution of females was skewed toward younger ages with few older, post-reproductive females. The distribution in recent years is more uniform across female ages (in other words, more females in their 30s, Fig. \ref{fig:plot-repro-females}).
```r"}
summarize repro females by year
repro = dplyr::filter(whaleData, pod %in% c("J1","K1","L1"), alive == 1, sexF1M2==1, age >=10, age<=42)
repro$pod[which(repro$pod=="J1")]="J"
repro$pod[which(repro$pod=="K1")]="K"
repro$pod[which(repro$pod=="L1")]="L"
repro$year = as.factor(repro$year)
ggplot(repro, aes(x = age, y = year, fill = ..x..)) +
geom_density_ridges_gradient(scale = 1, rel_min_height = 0.01) +
scale_fill_viridis(name = "Female age", option = "C") +
labs(title = 'Ages of reproductive - aged female SRKW') + xlab("Age") + ylab("Year") +
theme(axis.text.y = element_text(angle = 0, hjust = 1)) +
scale_y_discrete(breaks=c("1980","1985","1990","1995","2000","2005","2010","2015","2020"),
labels=c("1980","1985","1990","1995","2000","2005","2010","2015","2020")) +
coord_flip() +
theme_bw()
```r"}
g1 = group_by(repro, year) %>% summarize(m = length(unique(animal)),SRKW="Total") %>%
ggplot(aes(year, m, group=SRKW,color=SRKW)) + geom_line(col="black") + geom_point() +
xlab("Year") + ylab("Reproductive aged females") +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
theme_bw()
g2 = group_by(repro, year, pod) %>% summarize(m = length(unique(animal))) %>%
ggplot(aes(year, m, group=pod, color=pod)) + geom_line() + geom_point() +
xlab("Year") + ylab("Reproductive aged females") +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) + theme_bw() +
scale_color_viridis(discrete = TRUE, end = 0.8) +
scale_fill_viridis(discrete = TRUE, end = 0.8)
gridExtra::grid.arrange(g1, g2, ncol = 1)
For comparison, we can also look at the aggregate number of reproductive females in the NRKW population. This shows a nearly linear growth ovder time (Fig. \ref{fig:ts-reprofemales-nr}).
```r"}
repro = dplyr::filter(whaleData, pod %in% c("J1","K1","L1") ==FALSE, alive == 1, sexF1M2==1, age >=10, age<=42,
year<=2018)
g1 = group_by(repro, year) %>% summarize(m = length(unique(animal))) %>%
ggplot(aes(year, m, group=NA)) + geom_line(col="black") + geom_point() +
xlab("Year") + ylab("Reproductive aged females") +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
theme_bw()
g1
# Changing demographic rates
### Increased evidence of declining fecundity
There are several methodological changes from the projections done previously [@hilborn2012; @ward2013]. First, because we don't indices of salmon abundance available to whales (and none of the existing metrics of salmon abundance have been found to correlate with killer whale demography; PFMC Draft Report 2019), the estimation model was switched from a GLM to a GAM, which allows a smoother over year effects. In terms of interpretation, these smooth terms represent annual variation that may be driven by processes (prey, disease, etc) or changes in data quality or detectability (particularly for birth rates, which are partially a function of things like the fraction of time that particular social groups spend in inland waters where they're more likely to be seen). To evaluate changes in survival or fecundity rates over time, we can examine survival rates estimated using data through 2010, versus estimates using data through `r year.end`. For both instances, we can fit generalized additive models (GAMs) that include the effect of age (seperate splines by sex) and optionally a spline term over time. To not introduce bias, we did not allow animals born < 1970 to be included in the estimation (only the projection). For this sensitivity analysis, the SRKW and NRKW data are combined into a single model, with population as an estimated offset [@ward2013].
Results from these models indicate that of female survival, male survival, and female fecundity, fecundity rates have changed the most since 2010 and have declined (particularly for older ages, Fig. \ref{fig:demo-rates}). Time series of estimated fecundity rates also indicate low rates in recent years (this is not surprising, as there have been very few successful births). Since 1980, there have been several periods associated with high fecundity (late 1980s, mid 2000s, Fig. \ref{fig:ts-demo}). In contrast, estimated survival rates have been relatively flat (Fig. \ref{fig:ts-demo}), with the exception of low survival prior to the population being listed (SRKW reached a peak of 98 animals in 1995, and then dropped to 82 animals in 2003).
```r
# Have survival and fecundity rates changed since workshops?
whaleData = kwdemog::expand(orca)
whaleData = whaleData[whaleData$birth > 1970,]
# A handful of unknown sexes need to be randomly filled in
whaleData$sexF1M2[which(whaleData$sexF1M2==0)] = sample(c(1,2), size=length(which(whaleData$sexF1M2==0)), replace=T)
whaleData$region = ifelse(whaleData$pod %in% c("J1","K1","L1"), "SRKW", "NRKW")
whaleData$alive[which(whaleData$region=="NRKW" & whaleData$time > 2010)] = NA
survival_01 = gam(alive ~ s(age,by=as.factor(sexF1M2)) + s(year), family = "binomial", data = whaleData[which(whaleData$includeSurv == 1 & whaleData$year <= 2010), ])
survival_02 = gam(alive ~ s(age,by=as.factor(sexF1M2)) + s(year) + region, family = "binomial", data = whaleData[which(whaleData$includeSurv == 1), ])
survival_notime_01 = gam(alive ~ s(age,by=as.factor(sexF1M2)) + region, family = "binomial", data = whaleData[which(whaleData$includeSurv == 1 & whaleData$year <= 2010), ])
survival_notime_02 = gam(alive ~ s(age,by=as.factor(sexF1M2)) + region, family = "binomial", data = whaleData[which(whaleData$includeSurv == 1), ])
df = expand.grid("age"=1:40, "year"=2010, "region"="SRKW", "sexF1M2"=1:2, "model"=c("1970-2010",paste0("1970-",year.end)))
df$predict_surv = NA
df$predict_surv[df$model=="1970-2010"] = predict(survival_01, newdata=df[df$model=="1970-2010",], type="response")
df$predict_surv[df$model!="1970-2010"] = predict(survival_02, newdata=df[df$model!="1970-2010",], type="response")
df$predict_surv_notime = NA
df$predict_surv_notime[df$model=="1970-2010"] = predict(survival_notime_01, newdata=df[df$model=="1970-2010",], type="response")
df$predict_surv_notime[df$model!="1970-2010"] = predict(survival_notime_02, newdata=df[df$model!="1970-2010",], type="response")
g1 = ggplot(df[df$sexF1M2==1,], aes(age, predict_surv, group=model, color=model)) +
geom_line() + xlab("Age") + ylab("Female survivial") + ggtitle("Time included as predictor") +
theme_bw()
g1_notime = ggplot(df[df$sexF1M2==1,], aes(age, predict_surv_notime, group=model, color=model)) +
geom_line() + xlab("Age") + ylab("Female survivial") + ggtitle("Time not included as predictor") +
theme_bw()
g1_male = ggplot(df[df$sexF1M2==2 & df$age<=30,], aes(age, predict_surv, group=model, color=model)) +
geom_line() + xlab("Age") + ylab("Male survivial") + theme_bw()
g1_male_notime = ggplot(df[df$sexF1M2==2 & df$age<=30,], aes(age, predict_surv_notime, group=model, color=model)) +
geom_line() + xlab("Age") + ylab("Male survivial") + theme_bw()
# do same with fecundity
#whaleData$Birth[which(whaleData$region == "NRKW" & whaleData$year>2018)] = NA
whales_since76 = as.character(whaleData$animal[whaleData$birth > 1970])
fecundity_01 = gam(gave_birth ~ s(age,k=4) + s(year) + region, family = "binomial", data = whaleData[which(whaleData$animal%in%whales_since76 & whaleData$includeFec==1 & whaleData$sexF1M2=="1" & whaleData$age>=10 & whaleData$age< 43 & whaleData$year <= 2010), ])
fecundity_02 = gam(gave_birth ~ s(age,k=4) + s(year) + region, family = "binomial", data = whaleData[which(whaleData$animal%in%whales_since76 & whaleData$includeFec==1 & whaleData$sexF1M2=="1" & whaleData$age>=10 & whaleData$age< 43), ])
df$predict_fec=NA
df$predict_fec[df$model=="1970-2010"] = predict(fecundity_01, newdata=df[df$model=="1970-2010",], type="response")
df$predict_fec[df$model!="1970-2010"] = predict(fecundity_02, newdata=df[df$model!="1970-2010",], type="response")
g2 = ggplot(df, aes(age, predict_fec, group=model, color=model)) +
geom_line() + ylab("Predicted fecundity") + xlab("Age") + theme_bw()
fecundity_01_notime = gam(gave_birth ~ s(age,k=4) + region, family = "binomial", data = whaleData[which(whaleData$animal%in%whales_since76 & whaleData$includeFec==1 & whaleData$sexF1M2=="1" & whaleData$age>=10 & whaleData$age< 43 & whaleData$year <= 2010), ])
fecundity_02_notime = gam(gave_birth ~ s(age,k=4) + region, family = "binomial", data = whaleData[which(whaleData$animal%in%whales_since76 & whaleData$includeFec==1 & whaleData$sexF1M2=="1" & whaleData$age>=10 & whaleData$age< 43), ])
df$predict_fec=NA
df$predict_fec[df$model=="1970-2010"] = predict(fecundity_01_notime, newdata=df[df$model=="1970-2010",], type="response")
df$predict_fec[df$model!="1970-2010"] = predict(fecundity_02_notime, newdata=df[df$model!="1970-2010",], type="response")
g2_notime = ggplot(df, aes(age, predict_fec, group=model, color=model)) +
geom_line() + ylab("Predicted fecundity") + xlab("Age") + theme_bw()
```r"}
gridExtra::grid.arrange(g1, g1_notime, g1_male, g1_male_notime, g2, g2_notime, ncol=2)
```r",fig.height=8}
print(paste0(report_dir, "srkw_fec_age-year.rds"))
fecundity.model = readRDS(paste0(report_dir, "srkw_fec_age-year.rds"))
survival.model = readRDS(paste0(report_dir, "srkw_surv_age-year.rds"))
df = data.frame("age" = 20, "year"=1981:year.end, region="SRKW", "sexF1M2"=1)
df$pred = predict(fecundity.model, newdata=df, type="response")
df$low = predict(fecundity.model, newdata=df, type="response") - 2*predict(fecundity_02, newdata=df, type="response", se.fit=TRUE)$se.fit
df$hi = predict(fecundity.model, newdata=df, type="response") + 2*predict(fecundity_02, newdata=df, type="response", se.fit=TRUE)$se.fit
g = ggplot(df, aes(year, pred)) +
geom_ribbon(aes(ymin=low, ymax=hi), alpha=0.3, fill="dodgerblue") +
geom_line(color="dodgerblue") + ylab("Fecundity") +
xlab("Year") + theme_bw()
# same for survival
df = data.frame("age" = 20, "year"=1981:year.end, region="SRKW", "sexF1M2"=1,stage="young_female")
df$pred = predict(survival.model, newdata=df, type="response")
df$low = predict(survival.model, newdata=df, type="response") - 2*predict(survival_02, newdata=df, type="response", se.fit=TRUE)$se.fit
df$hi = predict(survival.model, newdata=df, type="response") + 2*predict(survival_02, newdata=df, type="response", se.fit=TRUE)$se.fit
g1 = ggplot(df, aes(year, pred)) +
geom_ribbon(aes(ymin=low, ymax=hi), alpha=0.3, fill="dodgerblue") +
geom_line(color="dodgerblue") + ylab("Female survival") +
xlab("Year") + theme_bw()
df = data.frame("age" = 20, "year"=1981:year.end, region="SRKW", "sexF1M2"=2,stage="young_male")
df$pred = predict(survival.model, newdata=df, type="response")
df$low = predict(survival.model, newdata=df, type="response") - 2*predict(survival.model, newdata=df, type="response", se.fit=TRUE)$se.fit
df$hi = predict(survival.model, newdata=df, type="response") + 2*predict(survival.model, newdata=df, type="response", se.fit=TRUE)$se.fit
g2 = ggplot(df, aes(year, pred)) +
geom_ribbon(aes(ymin=low, ymax=hi), alpha=0.3, fill="dodgerblue") +
geom_line(color="dodgerblue") + ylab("Male survival") +
xlab("Year") + theme_bw()
gridExtra::grid.arrange(g, g1, g2, ncol=1)
\break
Population projections
Given the current population age and sex structure, we performed a series of forecasts or projections, doing 1000 simulations of 25 years for each. Following previous work [@hilborn2012], projections beyond this time frame were not included, as longer term trajectories become negative, resulting in extinction or quasi-extinction. Following previous annual updates, we also used a sex ratio at birth of 55% male / 45% female, because of the historical skew in SRKW sex ratios at birth. Note that this differs from the previous reviews and published projectsions, which assumes a 50:50 sex ratio at birth [@hilborn2012; @ward2013].
The scenarios we considered were:
- Projections using fecundity rates and survival rates estimated over the entire time series (this is done with estimation models ignoring time)
- Projections using fecundity and survival rates estimated for the last 5 years,
r year.end-5 to r year.end
- Projections using the highest fecundity and survival rates estimated, in the period 1985-1989
Informative priors
Because this analysis was done in a Bayesian framework, we performed the analysis with and without informative prior distributions (scenarios without informative priors were based on Southern Resident data only). Data provided by DFO have updated the NRKW catalog from 2010-2011 to 2018. As such, we can (1) consider informative priors on the effect of age (for survival and fecundity rates), and (2) use informative priors on the effect of year. We included scenarios with just informative priors on age, or informative priors on the age and year terms (but not informative priors on just year alone).
Fecundity was modeled using a Bayesian GAM [@wood2016; @plummer2019], where
$$logit(p_{i,t})=B_{0} + f(year_{i,t}) + f(age_{i,t})$$
where $p_i$ represents the probability of animal $i$ giving birth at time $t$, $B_{0}$ is a population specific intercept, and smooth functions $f()$ are modeled over year and age. As in GAMs fit with maximum likelihood, the smooth functions can be written as $f(x)=\prod { j=1 }^{ J }{ { b }{ j }{ z }{ j }\left( x \right) }$, where $z{j}$ are basis functions for the smoother, $J$ is the dimensionality, and $b_{j}$ are estimated coefficients. Additional details are provided in [@wood2016].
Survival was also modeled using a logistic function, however we implemented a stage- rather than age-based model, following limitations of the data described in [@ward2013] and others. A primary concern for example is that ages of older animals at the start of the killer whale surveys were sometimes guessed based on reproductive histories, and potentially biased. The survival model predicting survival $s_{i,t}$ of animal $i$ at time $t$ can be written as
$$logit(s_{i,t})=B_{0} + B_{stage} + f(year_{i,t})$$
where $B_{0}$ is the population level offset, $B_{stage}$ is a stage-specific effect (stages as described in [@ward2013]) and $f(year_{i,t})$ the smooth term over years.
The fecundity and survival models were first fit to the NRKW data, using weakly informative priors, 50000 Markov chain Monte Carlo (MCMC) iterations, and 3 parallel MCMC chains. Output from the posterior distributions of these models were used to generate multivariate normal prior distributions for the SRKW population models. The effects of age or stage and year were separated, and the posterior mean and variance-covariance matrix of each was summarized.
if(refit_models == TRUE) {
# estimate models for NRKW
source(paste0(report_dir, "fit_nrkw.r"))
# use NRKW prior to do the same for SRKW
source(paste0(report_dir, "fit_srkw.r"))
# for each of the priors used for SRKW, do projections
#whaleData$Birth[which(whaleData$region== "NRKW" & whaleData$year>2018)] = NA
# (1) use age-year priors and status-quo projections
fecundity.model = readRDS(paste0(report_dir, "srkw_fec_age-year.rds"))
survival.model = readRDS(paste0(report_dir, "srkw_surv_age-year.rds"))
scenario = "status quo"
source(paste0(report_dir,"projections.R"))
save.image(paste0(report_dir,"projections_status-quo_age-year.Rdata"))
# (2) use age prior and status-quo projections
fecundity.model = readRDS(paste0(report_dir, "srkw_fec_age.rds"))
survival.model = readRDS(paste0(report_dir, "srkw_surv_age.rds"))
scenario = "status quo"
source(paste0(report_dir,"projections.R"))
save.image(paste0(report_dir,"projections_status-quo_age.Rdata"))
# (3) use age-year priors and last-5 projections
fecundity.model = readRDS(paste0(report_dir, "srkw_fec_age-year.rds"))
survival.model = readRDS(paste0(report_dir, "srkw_surv_age-year.rds"))
scenario = "last5"
source(paste0(report_dir,"projections.R"))
save.image(paste0(report_dir,"projections_last5_age-year.Rdata"))
# (4) use age prior and last-5 projections
fecundity.model = readRDS(paste0(report_dir, "srkw_fec_age.rds"))
survival.model = readRDS(paste0(report_dir, "srkw_surv_age.rds"))
scenario = "last5"
source(paste0(report_dir,"projections.R"))
save.image(paste0(report_dir,"projections_last5_age.Rdata"))
# (5) use age-year priors and last-5 projections
fecundity.model = readRDS(paste0(report_dir, "srkw_fec_age-year.rds"))
survival.model = readRDS(paste0(report_dir, "srkw_surv_age-year.rds"))
scenario = "best"
source(paste0(report_dir,"projections.R"))
save.image(paste0(report_dir,"projections_best_age-year.Rdata"))
# (6) use age prior and last-5 projections
fecundity.model = readRDS(paste0(report_dir, "srkw_fec_age.rds"))
survival.model = readRDS(paste0(report_dir, "srkw_surv_age.rds"))
scenario = "best"
source(paste0(report_dir,"projections.R"))
save.image(paste0(report_dir,"projections_best_age.Rdata"))
}
orca$age = year.end - orca$birth
orca_current = dplyr::filter(orca, !is.na(age), age < 25, pod %in% c("J001","K001","L001"), is.na(death), sexF1M2 == 1)
Projections done for the 2016 Status Review [@nmfs2016] also showed the population on a downward trajectory, and all three scenarios presented here provide further evidence of a declining population over the next 25 years (Fig. \ref{fig:proj1}). The most optimistic scenario, using demographic rates calculated from the 1985-1989 period, has a trajectory that increases and eventually declines after 2030. Additional runs for this scenario indicated a similar trajectory with a 50:50 sex ratio. Thus, this downward trend is likely driven by the current age and sex structure of young animals in the population, as well as the number of older animals. For example, there's currently r nrow(orca_current) females younger than 25 -- if all survive another 10 years, they will represent the majority of the reproductive females in the population and number fewer than the current population of reproductive females \ref{fig:plot-repro-females}.
A first set of projections represents the SRKW population projected using the NRKW age/stage priors, but not including priors on the year terms (essentially letting each population have different patterns over time). This
```r"}
f1 = paste0(report_dir,"projections_status-quo_age.Rdata")
f2 = paste0(report_dir,"projections_last5_age.Rdata")
f3 = paste0(report_dir,"projections_best_age.Rdata")
f4 = paste0(report_dir,"projections_status-quo_age-year.Rdata")
f5 = paste0(report_dir,"projections_last5_age-year.Rdata")
f6 = paste0(report_dir,"projections_best_age-year.Rdata")
load(f1)
d1 = data.frame("year"=year.end:(year.end+24), "mean"=apply(popSize,2,mean),
"low"=apply(popSize,2,quantile,0.05),
"hi"=apply(popSize,2,quantile,0.95),
"model"="All years")
load(f2)
y = paste0((year.end-5),"-",year.end)
d2 = data.frame("year"=year.end:(year.end+24), "mean"=apply(popSize,2,mean),
"low"=apply(popSize,2,quantile,0.05),
"hi"=apply(popSize,2,quantile,0.95),
"model"=paste0(year.end-4,"-",year.end))
load(f3)
d3 = data.frame("year"=year.end:(year.end+24), "mean"=apply(popSize,2,mean),
"low"=apply(popSize,2,quantile,0.05),
"hi"=apply(popSize,2,quantile,0.95),
"model"="1985-1989")
d = rbind(d1, d2, d3)
d = dplyr::rename(d, "Scenario"=model)
ggplot(d, aes(year, mean, group=Scenario, fill=Scenario,col=Scenario)) +
geom_ribbon(aes(ymin=low,ymax=hi),col=NA,alpha=0.3) +
geom_line() +
xlab("Year") +
ylab("Population size") +
theme_bw() +
scale_color_viridis(discrete = TRUE, end = 0.8) +
scale_fill_viridis(discrete = TRUE, end = 0.8)
The sensitivity comparing the SRKW projections under different prior distributions illustrated that for most scenarios, there appeared to be little effect of using the prior on year effects, in addition to age. Perhaps the scenario that showed the largest effect was the most optimistic scenario, using survival and birth rates from the 1985-1989 period (Fig. \ref{fig:proj2}).
```r"}
load(f4)
d1 = data.frame("year"=year.end:(year.end+24), "mean"=apply(popSize,2,mean),
"low"=apply(popSize,2,quantile,0.05),
"hi"=apply(popSize,2,quantile,0.95),
"model"="All years")
load(f5)
y = paste0((year.end-5),"-",year.end)
d2 = data.frame("year"=year.end:(year.end+24), "mean"=apply(popSize,2,mean),
"low"=apply(popSize,2,quantile,0.05),
"hi"=apply(popSize,2,quantile,0.95),
"model"=paste0(year.end-4,"-",year.end))
load(f6)
d3 = data.frame("year"=year.end:(year.end+24), "mean"=apply(popSize,2,mean),
"low"=apply(popSize,2,quantile,0.05),
"hi"=apply(popSize,2,quantile,0.95),
"model"="1985-1989")
dz = rbind(d1, d2, d3)
dz = dplyr::rename(dz, "Scenario"=model)
dz$prior = "age-year"
d$prior = "age"
dz = rbind(d, dz)
ggplot(dz, aes(year, mean, group=prior, fill=prior,col=prior)) +
geom_ribbon(aes(ymin=low,ymax=hi),col=NA,alpha=0.3) +
facet_wrap(~Scenario) +
geom_line() + xlab("Year") + ylab("Population size") + theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
scale_color_viridis(discrete = TRUE, end = 0.8) +
scale_fill_viridis(discrete = TRUE, end = 0.8)
It's also important to understand future changes to reproductive SRKW females (age 10-43). These plots use the same data as the projections for the whole population. The uncertainty bands around the projections are tightest for the next 10 years (those are based on animals currently alive in the population) -- but uncertainty increases in the future, as new simulated animals are added (Fig. \ref{fig:proj3}).
```r"}
load(f1)
d1 = data.frame("year"=year.end:(year.end+24), "mean"=apply(nFemales,2,mean),
"low"=apply(nFemales,2,quantile,0.05),
"hi"=apply(nFemales,2,quantile,0.95),
"model"="All years")
load(f2)
y = paste0((year.end-5),"-",year.end)
d2 = data.frame("year"=year.end:(year.end+24), "mean"=apply(nFemales,2,mean),
"low"=apply(nFemales,2,quantile,0.05),
"hi"=apply(nFemales,2,quantile,0.95),
"model"=paste0(year.end-4,"-",year.end))
load(f3)
d3 = data.frame("year"=year.end:(year.end+24), "mean"=apply(nFemales,2,mean),
"low"=apply(nFemales,2,quantile,0.05),
"hi"=apply(nFemales,2,quantile,0.95),
"model"="1985-1989")
d = rbind(d1, d2, d3)
d = dplyr::rename(d, "Scenario"=model)
ggplot(d, aes(year, mean, group=Scenario, fill=Scenario,col=Scenario)) +
geom_ribbon(aes(ymin=low,ymax=hi),col=NA,alpha=0.3) +
geom_line() +
xlab("Year") +
ylab("Reproductive females") +
theme_bw() +
scale_color_viridis(discrete = TRUE, end = 0.8) +
scale_fill_viridis(discrete = TRUE, end = 0.8)
### Variation in individual reproductive success
One of the factors that may contribute to SRKW declining faster than these projections represent is individual variation in reproductive success. All of the current modeling efforts assume that every female has the same probabilities of giving birth (adjusted for age and year). From an estimation standpoint, it's very difficult to fit more complex models that let individuals have unique rates, when reproductive histories are only partially observed. These generally result in estimates for young animals with wide confidence intervals.
```r
# Have survival and fecundity rates changed since workshops?
whaleData = kwdemog::expand(orca)
whaleData = whaleData[whaleData$birth > 1970,]
# A handful of unknown sexes need to be randomly filled in
whaleData$sexF1M2[which(whaleData$sexF1M2==0)] = sample(c(1,2), size=length(which(whaleData$sexF1M2==0)), replace=T)
whaleData$region = ifelse(whaleData$pod %in% c("J1","K1","L1"), "SRKW", "NRKW")
whales_since76 = as.character(whaleData$animal[whaleData$birth > 1970])
sub = whaleData[which(whaleData$animal%in%whales_since76 & whaleData$includeFec==1 & !is.na(whaleData$gave_birth) & whaleData$sexF1M2=="1" & whaleData$age>=10 & whaleData$age< 43), ]
sub$animal = as.factor(sub$animal)
fecundity_01 = gam(gave_birth ~ s(age,k=4) + s(year) + region + s(animal,bs="re"), family = "binomial", data = sub, method="REML")
extract_ranef(fecundity_01)
data(orca)
whaleData = kwdemog::expand(orca)
whaleData = whaleData[whaleData$birth > 1970,]
# A handful of unknown sexes need to be randomly filled in
whaleData$sexF1M2[which(whaleData$sexF1M2==0)] = sample(c(1,2), size=length(which(whaleData$sexF1M2==0)), replace=T)
whales_since76 = as.character(orca$animal[orca$birth > 1970])
sub = whaleData[which(whaleData$animal%in%whales_since76 & whaleData$population=="SRKW" & !is.na(whaleData$gave_birth) & whaleData$sexF1M2=="1" & whaleData$age>=10 & whaleData$age< 43), ]
sub$animal = as.factor(sub$animal)
# filter out females alive in current year
currently_alive = dplyr::filter(sub, year == year.end, is.na(death))
# find last birth for each of these
last_birth = group_by(dplyr::filter(sub, gave_birth==1), animal) %>%
dplyr::summarise(last_birth = max(year[which(gave_birth==1)]))
last_birth$last_birth[!is.finite(last_birth$last_birth)] = NA
currently_alive = dplyr::left_join(currently_alive, last_birth) %>%
dplyr::mutate(unlikely_future_mom="") %>%
dplyr::select(animal,age,last_birth,unlikely_future_mom) %>% data.frame()
currently_alive$unlikely_future_mom[which(currently_alive$last_birth <= (year.end-10))] = "*"
n_duds = length(which(currently_alive$unlikely_future_mom=="*"))
n_future = dim(dplyr::filter(currently_alive, age<=33, unlikely_future_mom!="*"))[1]
future_moms = dplyr::filter(orca, is.na(death), sexF1M2!=2, population=="SRKW", birth >= (year.end-10))
But we can look at the current potential moms in the population, and identify ones that haven't reproduced in the last decade. While there are r nrow(currently_alive) potential moms, r n_duds haven't given birth in the last decade, because of age or other reasons. Looking forward a decade and filtering these and older animals out, if the remaining all live there will be r n_future. Of the living SRKW animals who aren't yet of reproductive age, there are r length(which(future_moms$sexF1M2==1)) confirmed females and r length(which(future_moms$sexF1M2==0)) of unknown sex -- so if all these animals survive and are female, we'd expect to have r n_future+nrow(future_moms) SRKW moms in 10 years.
knitr::kable(currently_alive)
\break
Acknowledgments
Annual data for SRKW was collected by the Center for Whale Research. Recent updates to the NRKW catalog (initially provided by John Ford in 2011) and help interpreting data were provided by Thomas Doniol-Valcroze and Jared Towers.
References
nwfsc-cb/srkw-status documentation built on Jan. 16, 2025, 1 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = FALSE, warning=FALSE, message=FALSE, comment = "#>" ) library(kwdemog) library(ggplot2) library(knitr) library(dplyr) library(mgcv) library(rstanarm) library(viridis) library(ggridges) data(orca) year.end = 2024 refit_models = FALSE report_start = year.end-5 # format whale names just for consistency format_names = function(x) { for(i in 1:2) { indx = which(substr(x,2,2) == "0") x[indx] = paste0(substr(x[indx],1,1), substr(x[indx],3,nchar(x[indx]))) } return(x) } orca$animal = format_names(orca$animal) orca$pod = format_names(orca$pod) orca$matriline = format_names(orca$matriline) orca$mom = format_names(orca$mom) orca$dad = format_names(orca$dad) report_dir = "projections/"
Abstract
This document is intended to provide and overview and status update of the Southern Resident population following the most recent summer census in r year.end. Many of the analyses and figures presented below are described in more detail elsewhere, particularly the most recent Status Review [@nmfs2016] and recent technical reports following the 2011 - 2012 bilateral workshops [@hilborn2012; @ward2013]. These analyses also formed the backbone of the Pacific Fishery Management Council's working group on salmon fisheries impacts to killer whales [@pfmc2020; @pfmc2020b]. In particular, these analyses are meant to describe recent changes in population size and age structure, change in demographic rates over time, and an updated projections of population viability.
Recent births and deaths
Since r report_start, the following animals have been born in the SRKW population and lived:
whaleData = orca additions = whaleData[which(whaleData$pod%in%c("J1","K1","L1") & whaleData$birth > report_start & is.na(whaleData$death)),] additions = additions[,c("animal","birth","sexF1M2")] additions$sexF1M2[additions$sexF1M2=="0"]="Unk" additions$sexF1M2[additions$sexF1M2=="1"]="F" additions$sexF1M2[additions$sexF1M2=="2"]="M" names(additions)[3] = "sex" additions = dplyr::arrange(additions, birth, animal) additions = as.matrix(additions) row.names(additions) = NULL kable(as.matrix(additions))
The following whales have died since r report_start. Their deaths are also included in this table,
whaleData = orca losses = dplyr::filter(whaleData, pod%in%c("J1","K1","L1"), death > report_start) # subset columns losses = losses[,c("animal","death","sexF1M2","birth")] losses$age_death = losses$death-losses$birth losses$sexF1M2[losses$sexF1M2=="0"]="Unk" losses$sexF1M2[losses$sexF1M2=="1"]="F" losses$sexF1M2[losses$sexF1M2=="2"]="M" names(losses)[3] = "sex" losses = as.matrix(losses[,c("animal","death","sex","age_death")]) row.names(losses) = NULL # order losses = arrange(as.data.frame(losses), death, animal) %>% as.matrix() kable(losses)
Sex ratio at birth
During the 2011-2012 bilateral workshops [@hilborn2012; @ward2013], a comparison between NRKW and SRKW sex ratios at birth was presented, with calves being approximately 55\% female in the NRKW population and 45\% female in the SRKW population. This difference was assumed to be due to chance, and there was no evidence for a significant trend. As the proportion of males in the SRKW population has increased over time, it is worth re-examining the evidence supporting any trend.
sexr = dplyr::filter(whaleData, birth >= (year.end-3)) %>% dplyr::summarize(males = length(which(sexF1M2 == 2)), females = length(which(sexF1M2==1)))
# do a simple bayesian logistic regression indx = which(orca$sexF1M2 != 0 & orca$birth >= 1976 & orca$birth <= year.end) Y = orca$sexF1M2[indx] - 1 Xt = orca$birth[indx] df = data.frame("year"=Xt,"sex"=Y) g1 = ggplot(df, aes(year, sex)) + geom_hline(aes(yintercept=0.5),col="red",alpha=0.3) + geom_point(col="blue",alpha=0.3, position=position_dodge(width = 1)) + geom_smooth(method = "glm", method.args = list(family = "binomial")) + xlab("Year") + ylab("Sex (0 = F, 1 = M)") + theme_bw() # do a quick stepwise model selection using AIC #g.null = glm(Y ~ 1, family = "binomial") g.t = stan_glm(Y ~ Xt, family = "binomial") draws = as.matrix(g.t) ppos = 100 * round(length(which(draws[,2] > 0)) / nrow(draws), 3) df = data.frame("x" = draws[,2]) g2 = ggplot(data=df, aes(x=x)) + geom_histogram(fill = "dodgerblue",alpha=0.6) + xlab("Annual trend in male births") + ylab("Frequency") + geom_vline(aes(xintercept=0), col="grey30") + geom_vline(aes(xintercept=quantile(df$x,0.025)), col="red",alpha=0.3) + geom_vline(aes(xintercept=quantile(df$x,0.975)), col="red",alpha=0.3) + theme_bw()
To evaluate support for a trend, we fit Bayesian logistic regression models (GLMs with a binomial family and logit link function), to SRKW birth data over the period 1977-present. In recent years, since r year.end-3, there have been r sexr$males male births and r sexr$females female births. This analysis highlights that the probability of a positive trend is approximately r ppos\% (Fig. \ref{fig:srb}).
```r"} gridExtra::grid.arrange(g1, g2, ncol=1)
# Changing Population Structure
One of the objectives of the recovery goals [@nmfs2008; @nmfs2011] was an age and sex distribution that is more similar to the Northern Resident population (at the time of ESA listing). The Southern Resident population has undergone a number of shifts in age and sex, and because the population is so small, the age or sex composition are more sensitive to individual births and deaths. Previous status reviews based these targets from 1973-1996 [@olesiuk2005]:
```r
m = matrix("", 4, 2)
m[,1] = c("Juveniles", "Reproductive females", "Post-reproductive females", "Adult males")
m[,2] = c("47 %", "24 %", "11 %", "18 %")
colnames(m) = c("Stage","")
kable(m)
We can re-evaluate these targets, both for the Northern and Southern Resident populations, using the most recent years of data available (2018 for Northern Resident, r year.end for Southern Resident). There are some small differences between life stages used in Olesiuk et al. 2005, and more recent work [@ward2013]. First, Olesiuk assumed animals to not be mature to 15.5 (wheras Ward et al. 2013 assumed females to be mature at age 10). Second, Olesiuk et al. (2005) defined post-reproductive animals to be animals who hadn't given birth in 10 years (Ward et al. 2013 used a cutoff of 42 years). For these calculations, we define age at maturity to be 10 years, and 42+ as the age of reproductive senescence.
whaleData = kwdemog::expand(orca, start_year = 1976, current_year = year.end) # Sexes need to be randomly filled in whaleData$sexF1M2[which(whaleData$sexF1M2==0)] = sample(c(1,2), size=length(which(whaleData$sexF1M2==0)), replace=T) whaleData$age1979 = 1979 - whaleData$birth whaleData$age2010 = 2010 - whaleData$birth whaleData$ageCurrent = year.end - whaleData$birth NR = dplyr::filter(whaleData, year == 1979, pod %in% c("J1","K1","L1") == FALSE, alive==1, age1979 > 0) %>% dplyr::filter(death > 1979 | is.na(death)) NR_summary_1979 = dplyr::summarize(NR, x1 = length(which(age1979 < 10)), x2 = length(which(age1979 >= 10 & sexF1M2==2)), x3 = length(which(age1979 >= 10 & age1979 <= 42 & sexF1M2==1)), x4 = length(which(age1979 > 42 & sexF1M2==1))) %>% as.matrix() NR = dplyr::filter(whaleData, year == 2018, pod %in% c("J1","K1","L1") == FALSE, alive==1, ageCurrent > 0) %>% dplyr::filter(death > 2010 | is.na(death)) NR_summary = dplyr::summarize(NR, x1 = length(which(ageCurrent < 10)), x2 = length(which(ageCurrent >= 10 & sexF1M2==2)), x3 = length(which(ageCurrent >= 10 & age2010 <= 42 & sexF1M2==1)), x4 = length(which(ageCurrent > 42 & sexF1M2==1))) %>% as.matrix() SR = dplyr::filter(whaleData, year == 1979, pod %in% c("J1","K1","L1"), alive==1, age1979 > 0) %>% dplyr::filter(death > 1979 | is.na(death)) SR_summary_1979 = dplyr::summarize(SR, x1 = length(which(age1979 < 10)), x2 = length(which(age1979 >= 10 & sexF1M2==2)), x3 = length(which(age1979 >= 10 & age1979 <= 42 & sexF1M2==1)), x4 = length(which(age1979 > 42 & sexF1M2==1))) %>% as.matrix() SR = dplyr::filter(whaleData, year == year.end, pod %in% c("J1","K1","L1"), alive==1, ageCurrent > 0) %>% dplyr::filter(death > year.end | is.na(death)) SR_summary = dplyr::summarize(SR, x1 = length(which(ageCurrent < 10)), x2 = length(which(ageCurrent >= 10 & sexF1M2==2)), x3 = length(which(ageCurrent >= 10 & ageCurrent <= 42 & sexF1M2==1)), x4 = length(which(ageCurrent > 42 & sexF1M2==1))) %>% as.matrix() m = matrix(NA, 4, 4) m[,1] = 100*round(c(SR_summary_1979/sum(SR_summary_1979)), 2) m[,2] = 100*round(c(SR_summary/sum(SR_summary)), 2) m[,3] = 100*round(c(NR_summary_1979/sum(NR_summary_1979)), 2) m[,4] = 100*round(c(NR_summary/sum(NR_summary)), 2) for(i in 1:nrow(m)) { for(j in 1:ncol(m)) { m[i,j] = paste0(m[i,j], " %") } } rownames(m) = c("Juveniles (< 10)", "Adult males (10+)", "Adult females (10-42)", "Post-reproductive females (42+)") colnames(m) = c("SRKW 1979", paste0("SRKW ",year.end), "NRKW 1979", "NRKW 2018") kable(m)
Reproductive females
The number of reproductive aged females was at its lowest point in the late 1970s, in part because of the prior harvesting that occurred into the early 1970s (Fig. \ref{fig:ts-repro-females}). Though the overall number of reproductive females has been fluctuating between 25-35 for most of the last 40 years, there have been contrasting changes by pod, with declines in L pod females and increases in J pod (Fig. \ref{fig:ts-repro-females}). At the start of the survey in 1976, the distribution of females was skewed toward younger ages with few older, post-reproductive females. The distribution in recent years is more uniform across female ages (in other words, more females in their 30s, Fig. \ref{fig:plot-repro-females}).
```r"}
summarize repro females by year
repro = dplyr::filter(whaleData, pod %in% c("J1","K1","L1"), alive == 1, sexF1M2==1, age >=10, age<=42) repro$pod[which(repro$pod=="J1")]="J" repro$pod[which(repro$pod=="K1")]="K" repro$pod[which(repro$pod=="L1")]="L" repro$year = as.factor(repro$year)
ggplot(repro, aes(x = age, y = year, fill = ..x..)) +
geom_density_ridges_gradient(scale = 1, rel_min_height = 0.01) +
scale_fill_viridis(name = "Female age", option = "C") +
labs(title = 'Ages of reproductive - aged female SRKW') + xlab("Age") + ylab("Year") +
theme(axis.text.y = element_text(angle = 0, hjust = 1)) +
scale_y_discrete(breaks=c("1980","1985","1990","1995","2000","2005","2010","2015","2020"),
labels=c("1980","1985","1990","1995","2000","2005","2010","2015","2020")) +
coord_flip() +
theme_bw()
```r"}
g1 = group_by(repro, year) %>% summarize(m = length(unique(animal)),SRKW="Total") %>%
ggplot(aes(year, m, group=SRKW,color=SRKW)) + geom_line(col="black") + geom_point() +
xlab("Year") + ylab("Reproductive aged females") +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
theme_bw()
g2 = group_by(repro, year, pod) %>% summarize(m = length(unique(animal))) %>%
ggplot(aes(year, m, group=pod, color=pod)) + geom_line() + geom_point() +
xlab("Year") + ylab("Reproductive aged females") +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) + theme_bw() +
scale_color_viridis(discrete = TRUE, end = 0.8) +
scale_fill_viridis(discrete = TRUE, end = 0.8)
gridExtra::grid.arrange(g1, g2, ncol = 1)
For comparison, we can also look at the aggregate number of reproductive females in the NRKW population. This shows a nearly linear growth ovder time (Fig. \ref{fig:ts-reprofemales-nr}).
```r"} repro = dplyr::filter(whaleData, pod %in% c("J1","K1","L1") ==FALSE, alive == 1, sexF1M2==1, age >=10, age<=42, year<=2018)
g1 = group_by(repro, year) %>% summarize(m = length(unique(animal))) %>% ggplot(aes(year, m, group=NA)) + geom_line(col="black") + geom_point() + xlab("Year") + ylab("Reproductive aged females") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + theme_bw() g1
# Changing demographic rates
### Increased evidence of declining fecundity
There are several methodological changes from the projections done previously [@hilborn2012; @ward2013]. First, because we don't indices of salmon abundance available to whales (and none of the existing metrics of salmon abundance have been found to correlate with killer whale demography; PFMC Draft Report 2019), the estimation model was switched from a GLM to a GAM, which allows a smoother over year effects. In terms of interpretation, these smooth terms represent annual variation that may be driven by processes (prey, disease, etc) or changes in data quality or detectability (particularly for birth rates, which are partially a function of things like the fraction of time that particular social groups spend in inland waters where they're more likely to be seen). To evaluate changes in survival or fecundity rates over time, we can examine survival rates estimated using data through 2010, versus estimates using data through `r year.end`. For both instances, we can fit generalized additive models (GAMs) that include the effect of age (seperate splines by sex) and optionally a spline term over time. To not introduce bias, we did not allow animals born < 1970 to be included in the estimation (only the projection). For this sensitivity analysis, the SRKW and NRKW data are combined into a single model, with population as an estimated offset [@ward2013].
Results from these models indicate that of female survival, male survival, and female fecundity, fecundity rates have changed the most since 2010 and have declined (particularly for older ages, Fig. \ref{fig:demo-rates}). Time series of estimated fecundity rates also indicate low rates in recent years (this is not surprising, as there have been very few successful births). Since 1980, there have been several periods associated with high fecundity (late 1980s, mid 2000s, Fig. \ref{fig:ts-demo}). In contrast, estimated survival rates have been relatively flat (Fig. \ref{fig:ts-demo}), with the exception of low survival prior to the population being listed (SRKW reached a peak of 98 animals in 1995, and then dropped to 82 animals in 2003).
```r
# Have survival and fecundity rates changed since workshops?
whaleData = kwdemog::expand(orca)
whaleData = whaleData[whaleData$birth > 1970,]
# A handful of unknown sexes need to be randomly filled in
whaleData$sexF1M2[which(whaleData$sexF1M2==0)] = sample(c(1,2), size=length(which(whaleData$sexF1M2==0)), replace=T)
whaleData$region = ifelse(whaleData$pod %in% c("J1","K1","L1"), "SRKW", "NRKW")
whaleData$alive[which(whaleData$region=="NRKW" & whaleData$time > 2010)] = NA
survival_01 = gam(alive ~ s(age,by=as.factor(sexF1M2)) + s(year), family = "binomial", data = whaleData[which(whaleData$includeSurv == 1 & whaleData$year <= 2010), ])
survival_02 = gam(alive ~ s(age,by=as.factor(sexF1M2)) + s(year) + region, family = "binomial", data = whaleData[which(whaleData$includeSurv == 1), ])
survival_notime_01 = gam(alive ~ s(age,by=as.factor(sexF1M2)) + region, family = "binomial", data = whaleData[which(whaleData$includeSurv == 1 & whaleData$year <= 2010), ])
survival_notime_02 = gam(alive ~ s(age,by=as.factor(sexF1M2)) + region, family = "binomial", data = whaleData[which(whaleData$includeSurv == 1), ])
df = expand.grid("age"=1:40, "year"=2010, "region"="SRKW", "sexF1M2"=1:2, "model"=c("1970-2010",paste0("1970-",year.end)))
df$predict_surv = NA
df$predict_surv[df$model=="1970-2010"] = predict(survival_01, newdata=df[df$model=="1970-2010",], type="response")
df$predict_surv[df$model!="1970-2010"] = predict(survival_02, newdata=df[df$model!="1970-2010",], type="response")
df$predict_surv_notime = NA
df$predict_surv_notime[df$model=="1970-2010"] = predict(survival_notime_01, newdata=df[df$model=="1970-2010",], type="response")
df$predict_surv_notime[df$model!="1970-2010"] = predict(survival_notime_02, newdata=df[df$model!="1970-2010",], type="response")
g1 = ggplot(df[df$sexF1M2==1,], aes(age, predict_surv, group=model, color=model)) +
geom_line() + xlab("Age") + ylab("Female survivial") + ggtitle("Time included as predictor") +
theme_bw()
g1_notime = ggplot(df[df$sexF1M2==1,], aes(age, predict_surv_notime, group=model, color=model)) +
geom_line() + xlab("Age") + ylab("Female survivial") + ggtitle("Time not included as predictor") +
theme_bw()
g1_male = ggplot(df[df$sexF1M2==2 & df$age<=30,], aes(age, predict_surv, group=model, color=model)) +
geom_line() + xlab("Age") + ylab("Male survivial") + theme_bw()
g1_male_notime = ggplot(df[df$sexF1M2==2 & df$age<=30,], aes(age, predict_surv_notime, group=model, color=model)) +
geom_line() + xlab("Age") + ylab("Male survivial") + theme_bw()
# do same with fecundity
#whaleData$Birth[which(whaleData$region == "NRKW" & whaleData$year>2018)] = NA
whales_since76 = as.character(whaleData$animal[whaleData$birth > 1970])
fecundity_01 = gam(gave_birth ~ s(age,k=4) + s(year) + region, family = "binomial", data = whaleData[which(whaleData$animal%in%whales_since76 & whaleData$includeFec==1 & whaleData$sexF1M2=="1" & whaleData$age>=10 & whaleData$age< 43 & whaleData$year <= 2010), ])
fecundity_02 = gam(gave_birth ~ s(age,k=4) + s(year) + region, family = "binomial", data = whaleData[which(whaleData$animal%in%whales_since76 & whaleData$includeFec==1 & whaleData$sexF1M2=="1" & whaleData$age>=10 & whaleData$age< 43), ])
df$predict_fec=NA
df$predict_fec[df$model=="1970-2010"] = predict(fecundity_01, newdata=df[df$model=="1970-2010",], type="response")
df$predict_fec[df$model!="1970-2010"] = predict(fecundity_02, newdata=df[df$model!="1970-2010",], type="response")
g2 = ggplot(df, aes(age, predict_fec, group=model, color=model)) +
geom_line() + ylab("Predicted fecundity") + xlab("Age") + theme_bw()
fecundity_01_notime = gam(gave_birth ~ s(age,k=4) + region, family = "binomial", data = whaleData[which(whaleData$animal%in%whales_since76 & whaleData$includeFec==1 & whaleData$sexF1M2=="1" & whaleData$age>=10 & whaleData$age< 43 & whaleData$year <= 2010), ])
fecundity_02_notime = gam(gave_birth ~ s(age,k=4) + region, family = "binomial", data = whaleData[which(whaleData$animal%in%whales_since76 & whaleData$includeFec==1 & whaleData$sexF1M2=="1" & whaleData$age>=10 & whaleData$age< 43), ])
df$predict_fec=NA
df$predict_fec[df$model=="1970-2010"] = predict(fecundity_01_notime, newdata=df[df$model=="1970-2010",], type="response")
df$predict_fec[df$model!="1970-2010"] = predict(fecundity_02_notime, newdata=df[df$model!="1970-2010",], type="response")
g2_notime = ggplot(df, aes(age, predict_fec, group=model, color=model)) +
geom_line() + ylab("Predicted fecundity") + xlab("Age") + theme_bw()
```r"}
gridExtra::grid.arrange(g1, g1_notime, g1_male, g1_male_notime, g2, g2_notime, ncol=2)
```r",fig.height=8}
print(paste0(report_dir, "srkw_fec_age-year.rds"))
fecundity.model = readRDS(paste0(report_dir, "srkw_fec_age-year.rds"))
survival.model = readRDS(paste0(report_dir, "srkw_surv_age-year.rds"))
df = data.frame("age" = 20, "year"=1981:year.end, region="SRKW", "sexF1M2"=1)
df$pred = predict(fecundity.model, newdata=df, type="response")
df$low = predict(fecundity.model, newdata=df, type="response") - 2*predict(fecundity_02, newdata=df, type="response", se.fit=TRUE)$se.fit
df$hi = predict(fecundity.model, newdata=df, type="response") + 2*predict(fecundity_02, newdata=df, type="response", se.fit=TRUE)$se.fit
g = ggplot(df, aes(year, pred)) +
geom_ribbon(aes(ymin=low, ymax=hi), alpha=0.3, fill="dodgerblue") +
geom_line(color="dodgerblue") + ylab("Fecundity") +
xlab("Year") + theme_bw()
# same for survival
df = data.frame("age" = 20, "year"=1981:year.end, region="SRKW", "sexF1M2"=1,stage="young_female")
df$pred = predict(survival.model, newdata=df, type="response")
df$low = predict(survival.model, newdata=df, type="response") - 2*predict(survival_02, newdata=df, type="response", se.fit=TRUE)$se.fit
df$hi = predict(survival.model, newdata=df, type="response") + 2*predict(survival_02, newdata=df, type="response", se.fit=TRUE)$se.fit
g1 = ggplot(df, aes(year, pred)) +
geom_ribbon(aes(ymin=low, ymax=hi), alpha=0.3, fill="dodgerblue") +
geom_line(color="dodgerblue") + ylab("Female survival") +
xlab("Year") + theme_bw()
df = data.frame("age" = 20, "year"=1981:year.end, region="SRKW", "sexF1M2"=2,stage="young_male")
df$pred = predict(survival.model, newdata=df, type="response")
df$low = predict(survival.model, newdata=df, type="response") - 2*predict(survival.model, newdata=df, type="response", se.fit=TRUE)$se.fit
df$hi = predict(survival.model, newdata=df, type="response") + 2*predict(survival.model, newdata=df, type="response", se.fit=TRUE)$se.fit
g2 = ggplot(df, aes(year, pred)) +
geom_ribbon(aes(ymin=low, ymax=hi), alpha=0.3, fill="dodgerblue") +
geom_line(color="dodgerblue") + ylab("Male survival") +
xlab("Year") + theme_bw()
gridExtra::grid.arrange(g, g1, g2, ncol=1)
\break
Population projections
Given the current population age and sex structure, we performed a series of forecasts or projections, doing 1000 simulations of 25 years for each. Following previous work [@hilborn2012], projections beyond this time frame were not included, as longer term trajectories become negative, resulting in extinction or quasi-extinction. Following previous annual updates, we also used a sex ratio at birth of 55% male / 45% female, because of the historical skew in SRKW sex ratios at birth. Note that this differs from the previous reviews and published projectsions, which assumes a 50:50 sex ratio at birth [@hilborn2012; @ward2013].
The scenarios we considered were:
- Projections using fecundity rates and survival rates estimated over the entire time series (this is done with estimation models ignoring time)
- Projections using fecundity and survival rates estimated for the last 5 years,
r year.end-5tor year.end - Projections using the highest fecundity and survival rates estimated, in the period 1985-1989
Informative priors
Because this analysis was done in a Bayesian framework, we performed the analysis with and without informative prior distributions (scenarios without informative priors were based on Southern Resident data only). Data provided by DFO have updated the NRKW catalog from 2010-2011 to 2018. As such, we can (1) consider informative priors on the effect of age (for survival and fecundity rates), and (2) use informative priors on the effect of year. We included scenarios with just informative priors on age, or informative priors on the age and year terms (but not informative priors on just year alone).
Fecundity was modeled using a Bayesian GAM [@wood2016; @plummer2019], where
$$logit(p_{i,t})=B_{0} + f(year_{i,t}) + f(age_{i,t})$$ where $p_i$ represents the probability of animal $i$ giving birth at time $t$, $B_{0}$ is a population specific intercept, and smooth functions $f()$ are modeled over year and age. As in GAMs fit with maximum likelihood, the smooth functions can be written as $f(x)=\prod { j=1 }^{ J }{ { b }{ j }{ z }{ j }\left( x \right) }$, where $z{j}$ are basis functions for the smoother, $J$ is the dimensionality, and $b_{j}$ are estimated coefficients. Additional details are provided in [@wood2016].
Survival was also modeled using a logistic function, however we implemented a stage- rather than age-based model, following limitations of the data described in [@ward2013] and others. A primary concern for example is that ages of older animals at the start of the killer whale surveys were sometimes guessed based on reproductive histories, and potentially biased. The survival model predicting survival $s_{i,t}$ of animal $i$ at time $t$ can be written as
$$logit(s_{i,t})=B_{0} + B_{stage} + f(year_{i,t})$$ where $B_{0}$ is the population level offset, $B_{stage}$ is a stage-specific effect (stages as described in [@ward2013]) and $f(year_{i,t})$ the smooth term over years.
The fecundity and survival models were first fit to the NRKW data, using weakly informative priors, 50000 Markov chain Monte Carlo (MCMC) iterations, and 3 parallel MCMC chains. Output from the posterior distributions of these models were used to generate multivariate normal prior distributions for the SRKW population models. The effects of age or stage and year were separated, and the posterior mean and variance-covariance matrix of each was summarized.
if(refit_models == TRUE) { # estimate models for NRKW source(paste0(report_dir, "fit_nrkw.r")) # use NRKW prior to do the same for SRKW source(paste0(report_dir, "fit_srkw.r")) # for each of the priors used for SRKW, do projections #whaleData$Birth[which(whaleData$region== "NRKW" & whaleData$year>2018)] = NA # (1) use age-year priors and status-quo projections fecundity.model = readRDS(paste0(report_dir, "srkw_fec_age-year.rds")) survival.model = readRDS(paste0(report_dir, "srkw_surv_age-year.rds")) scenario = "status quo" source(paste0(report_dir,"projections.R")) save.image(paste0(report_dir,"projections_status-quo_age-year.Rdata")) # (2) use age prior and status-quo projections fecundity.model = readRDS(paste0(report_dir, "srkw_fec_age.rds")) survival.model = readRDS(paste0(report_dir, "srkw_surv_age.rds")) scenario = "status quo" source(paste0(report_dir,"projections.R")) save.image(paste0(report_dir,"projections_status-quo_age.Rdata")) # (3) use age-year priors and last-5 projections fecundity.model = readRDS(paste0(report_dir, "srkw_fec_age-year.rds")) survival.model = readRDS(paste0(report_dir, "srkw_surv_age-year.rds")) scenario = "last5" source(paste0(report_dir,"projections.R")) save.image(paste0(report_dir,"projections_last5_age-year.Rdata")) # (4) use age prior and last-5 projections fecundity.model = readRDS(paste0(report_dir, "srkw_fec_age.rds")) survival.model = readRDS(paste0(report_dir, "srkw_surv_age.rds")) scenario = "last5" source(paste0(report_dir,"projections.R")) save.image(paste0(report_dir,"projections_last5_age.Rdata")) # (5) use age-year priors and last-5 projections fecundity.model = readRDS(paste0(report_dir, "srkw_fec_age-year.rds")) survival.model = readRDS(paste0(report_dir, "srkw_surv_age-year.rds")) scenario = "best" source(paste0(report_dir,"projections.R")) save.image(paste0(report_dir,"projections_best_age-year.Rdata")) # (6) use age prior and last-5 projections fecundity.model = readRDS(paste0(report_dir, "srkw_fec_age.rds")) survival.model = readRDS(paste0(report_dir, "srkw_surv_age.rds")) scenario = "best" source(paste0(report_dir,"projections.R")) save.image(paste0(report_dir,"projections_best_age.Rdata")) }
orca$age = year.end - orca$birth orca_current = dplyr::filter(orca, !is.na(age), age < 25, pod %in% c("J001","K001","L001"), is.na(death), sexF1M2 == 1)
Projections done for the 2016 Status Review [@nmfs2016] also showed the population on a downward trajectory, and all three scenarios presented here provide further evidence of a declining population over the next 25 years (Fig. \ref{fig:proj1}). The most optimistic scenario, using demographic rates calculated from the 1985-1989 period, has a trajectory that increases and eventually declines after 2030. Additional runs for this scenario indicated a similar trajectory with a 50:50 sex ratio. Thus, this downward trend is likely driven by the current age and sex structure of young animals in the population, as well as the number of older animals. For example, there's currently r nrow(orca_current) females younger than 25 -- if all survive another 10 years, they will represent the majority of the reproductive females in the population and number fewer than the current population of reproductive females \ref{fig:plot-repro-females}.
A first set of projections represents the SRKW population projected using the NRKW age/stage priors, but not including priors on the year terms (essentially letting each population have different patterns over time). This
```r"}
f1 = paste0(report_dir,"projections_status-quo_age.Rdata") f2 = paste0(report_dir,"projections_last5_age.Rdata") f3 = paste0(report_dir,"projections_best_age.Rdata") f4 = paste0(report_dir,"projections_status-quo_age-year.Rdata") f5 = paste0(report_dir,"projections_last5_age-year.Rdata") f6 = paste0(report_dir,"projections_best_age-year.Rdata") load(f1)
d1 = data.frame("year"=year.end:(year.end+24), "mean"=apply(popSize,2,mean), "low"=apply(popSize,2,quantile,0.05), "hi"=apply(popSize,2,quantile,0.95), "model"="All years")
load(f2) y = paste0((year.end-5),"-",year.end) d2 = data.frame("year"=year.end:(year.end+24), "mean"=apply(popSize,2,mean), "low"=apply(popSize,2,quantile,0.05), "hi"=apply(popSize,2,quantile,0.95), "model"=paste0(year.end-4,"-",year.end))
load(f3) d3 = data.frame("year"=year.end:(year.end+24), "mean"=apply(popSize,2,mean), "low"=apply(popSize,2,quantile,0.05), "hi"=apply(popSize,2,quantile,0.95), "model"="1985-1989")
d = rbind(d1, d2, d3) d = dplyr::rename(d, "Scenario"=model)
ggplot(d, aes(year, mean, group=Scenario, fill=Scenario,col=Scenario)) + geom_ribbon(aes(ymin=low,ymax=hi),col=NA,alpha=0.3) + geom_line() + xlab("Year") + ylab("Population size") + theme_bw() + scale_color_viridis(discrete = TRUE, end = 0.8) + scale_fill_viridis(discrete = TRUE, end = 0.8)
The sensitivity comparing the SRKW projections under different prior distributions illustrated that for most scenarios, there appeared to be little effect of using the prior on year effects, in addition to age. Perhaps the scenario that showed the largest effect was the most optimistic scenario, using survival and birth rates from the 1985-1989 period (Fig. \ref{fig:proj2}).
```r"}
load(f4)
d1 = data.frame("year"=year.end:(year.end+24), "mean"=apply(popSize,2,mean),
"low"=apply(popSize,2,quantile,0.05),
"hi"=apply(popSize,2,quantile,0.95),
"model"="All years")
load(f5)
y = paste0((year.end-5),"-",year.end)
d2 = data.frame("year"=year.end:(year.end+24), "mean"=apply(popSize,2,mean),
"low"=apply(popSize,2,quantile,0.05),
"hi"=apply(popSize,2,quantile,0.95),
"model"=paste0(year.end-4,"-",year.end))
load(f6)
d3 = data.frame("year"=year.end:(year.end+24), "mean"=apply(popSize,2,mean),
"low"=apply(popSize,2,quantile,0.05),
"hi"=apply(popSize,2,quantile,0.95),
"model"="1985-1989")
dz = rbind(d1, d2, d3)
dz = dplyr::rename(dz, "Scenario"=model)
dz$prior = "age-year"
d$prior = "age"
dz = rbind(d, dz)
ggplot(dz, aes(year, mean, group=prior, fill=prior,col=prior)) +
geom_ribbon(aes(ymin=low,ymax=hi),col=NA,alpha=0.3) +
facet_wrap(~Scenario) +
geom_line() + xlab("Year") + ylab("Population size") + theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
scale_color_viridis(discrete = TRUE, end = 0.8) +
scale_fill_viridis(discrete = TRUE, end = 0.8)
It's also important to understand future changes to reproductive SRKW females (age 10-43). These plots use the same data as the projections for the whole population. The uncertainty bands around the projections are tightest for the next 10 years (those are based on animals currently alive in the population) -- but uncertainty increases in the future, as new simulated animals are added (Fig. \ref{fig:proj3}).
```r"}
load(f1)
d1 = data.frame("year"=year.end:(year.end+24), "mean"=apply(nFemales,2,mean), "low"=apply(nFemales,2,quantile,0.05), "hi"=apply(nFemales,2,quantile,0.95), "model"="All years")
load(f2) y = paste0((year.end-5),"-",year.end) d2 = data.frame("year"=year.end:(year.end+24), "mean"=apply(nFemales,2,mean), "low"=apply(nFemales,2,quantile,0.05), "hi"=apply(nFemales,2,quantile,0.95), "model"=paste0(year.end-4,"-",year.end))
load(f3) d3 = data.frame("year"=year.end:(year.end+24), "mean"=apply(nFemales,2,mean), "low"=apply(nFemales,2,quantile,0.05), "hi"=apply(nFemales,2,quantile,0.95), "model"="1985-1989")
d = rbind(d1, d2, d3) d = dplyr::rename(d, "Scenario"=model)
ggplot(d, aes(year, mean, group=Scenario, fill=Scenario,col=Scenario)) + geom_ribbon(aes(ymin=low,ymax=hi),col=NA,alpha=0.3) + geom_line() + xlab("Year") + ylab("Reproductive females") + theme_bw() + scale_color_viridis(discrete = TRUE, end = 0.8) + scale_fill_viridis(discrete = TRUE, end = 0.8)
### Variation in individual reproductive success
One of the factors that may contribute to SRKW declining faster than these projections represent is individual variation in reproductive success. All of the current modeling efforts assume that every female has the same probabilities of giving birth (adjusted for age and year). From an estimation standpoint, it's very difficult to fit more complex models that let individuals have unique rates, when reproductive histories are only partially observed. These generally result in estimates for young animals with wide confidence intervals.
```r
# Have survival and fecundity rates changed since workshops?
whaleData = kwdemog::expand(orca)
whaleData = whaleData[whaleData$birth > 1970,]
# A handful of unknown sexes need to be randomly filled in
whaleData$sexF1M2[which(whaleData$sexF1M2==0)] = sample(c(1,2), size=length(which(whaleData$sexF1M2==0)), replace=T)
whaleData$region = ifelse(whaleData$pod %in% c("J1","K1","L1"), "SRKW", "NRKW")
whales_since76 = as.character(whaleData$animal[whaleData$birth > 1970])
sub = whaleData[which(whaleData$animal%in%whales_since76 & whaleData$includeFec==1 & !is.na(whaleData$gave_birth) & whaleData$sexF1M2=="1" & whaleData$age>=10 & whaleData$age< 43), ]
sub$animal = as.factor(sub$animal)
fecundity_01 = gam(gave_birth ~ s(age,k=4) + s(year) + region + s(animal,bs="re"), family = "binomial", data = sub, method="REML")
extract_ranef(fecundity_01)
data(orca) whaleData = kwdemog::expand(orca) whaleData = whaleData[whaleData$birth > 1970,] # A handful of unknown sexes need to be randomly filled in whaleData$sexF1M2[which(whaleData$sexF1M2==0)] = sample(c(1,2), size=length(which(whaleData$sexF1M2==0)), replace=T) whales_since76 = as.character(orca$animal[orca$birth > 1970]) sub = whaleData[which(whaleData$animal%in%whales_since76 & whaleData$population=="SRKW" & !is.na(whaleData$gave_birth) & whaleData$sexF1M2=="1" & whaleData$age>=10 & whaleData$age< 43), ] sub$animal = as.factor(sub$animal) # filter out females alive in current year currently_alive = dplyr::filter(sub, year == year.end, is.na(death)) # find last birth for each of these last_birth = group_by(dplyr::filter(sub, gave_birth==1), animal) %>% dplyr::summarise(last_birth = max(year[which(gave_birth==1)])) last_birth$last_birth[!is.finite(last_birth$last_birth)] = NA currently_alive = dplyr::left_join(currently_alive, last_birth) %>% dplyr::mutate(unlikely_future_mom="") %>% dplyr::select(animal,age,last_birth,unlikely_future_mom) %>% data.frame() currently_alive$unlikely_future_mom[which(currently_alive$last_birth <= (year.end-10))] = "*" n_duds = length(which(currently_alive$unlikely_future_mom=="*")) n_future = dim(dplyr::filter(currently_alive, age<=33, unlikely_future_mom!="*"))[1] future_moms = dplyr::filter(orca, is.na(death), sexF1M2!=2, population=="SRKW", birth >= (year.end-10))
But we can look at the current potential moms in the population, and identify ones that haven't reproduced in the last decade. While there are r nrow(currently_alive) potential moms, r n_duds haven't given birth in the last decade, because of age or other reasons. Looking forward a decade and filtering these and older animals out, if the remaining all live there will be r n_future. Of the living SRKW animals who aren't yet of reproductive age, there are r length(which(future_moms$sexF1M2==1)) confirmed females and r length(which(future_moms$sexF1M2==0)) of unknown sex -- so if all these animals survive and are female, we'd expect to have r n_future+nrow(future_moms) SRKW moms in 10 years.
knitr::kable(currently_alive)
\break
Acknowledgments
Annual data for SRKW was collected by the Center for Whale Research. Recent updates to the NRKW catalog (initially provided by John Ford in 2011) and help interpreting data were provided by Thomas Doniol-Valcroze and Jared Towers.
References
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.