Introduction

This code works with a dataset on the rare plant Synthyris bullii (formerly Besseya bullii), commonly called kittentails. The data comes from the Chicago Botanic Garden’s Plants of Concern Program, and includes data from 3 separate populations. To protect the state-threatened species, the coordinates of the populations were removed before data was shared. Running the code will briefly explore the data, calculate population growth metrics, test for density dependence, and then run a population viability analysis using a density independent model. Although the population of kittentails does exhibit negative density dependence, for the sake of time and more as an exercise in learning, the density indpendent PVA is sufficient.

Set Up and Exploration

Set working directory

knitr::opts_knit$set(root.dir = 'C:/Users/adavi/Documents/Northwestern/Thesis/Besseya Data')
knitr::opts_chunk$set(echo = TRUE)

Load packages

library(tidyverse)
library(popbio)
library(lubridate)
library(kableExtra)
library(knitr)

#Cite all loaded packages (citations will automatically be added to end of document)
knitr::write_bib(c(.packages()), "packages.bib")

Read in original data

bbull_poc <- read.csv("besseyaBullii_monitoringReports_2021.12.07.csv")

Delete GPS data - this is the shared data.

sbull <- bbull_poc %>% 
  select(-gps_original, 
         -gps_decimal_degrees_wgs84,
         -directions,
         -monitors)

write.csv(sbull, "C:/Users/adavi/Documents/Northwestern/Thesis/Besseya Data/sbull_nocoords", row.names = TRUE)

Explore the data

head(sbull)
tail(sbull)
ncol(sbull)
nrow(sbull)
str(sbull)

Reformat and Calculate Population Growth Metrics

Select single population, add year column

sbull2 <- sbull %>% 
  filter(site_name == "LeRoy Oakes Forest Preserve") ## Should I rename the location?

Convert date monitored to “date” format and add ‘Year’ column with lubridate()

sbull2$date_monitored <- mdy(sbull2$date_monitored)

sbull3 <- sbull2 %>% 
          mutate(year = year(date_monitored))

Simplify dataset - select needed columns only

sbull4 <- sbull3 %>% 
  select(scientific_name, county, year, total_number_counted, count_range, percent_reproductive)

head(sbull4)
tail(sbull4)

Create function to shift cell values

shift <- function(x, n){
  c(x[-(seq(n))], rep(NA, n))
}

Add column with following year count to calculate growth rate

sbull5 <- sbull4 %>% 
          mutate(next_year_count = shift(sbull4$total_number_counted, 1)) %>% 
          relocate(next_year_count, .after = total_number_counted)

Add column for annual growth rate

sbull6 <- sbull5 %>% 
          mutate(growth_rate = (next_year_count - total_number_counted)/(total_number_counted))
sbull6

Visualize Data and Assess Density Dependence

Plot population size over time

sbull_plot1 <- sbull6 %>% 
               ggplot(aes(year, total_number_counted)) +
               geom_point() +
               geom_smooth(method = "loess", se = F) +
               labs(x = "Year",
                    y = "Population Size (N)",
                    caption = "Figure 1. Size of a Synthyris bullii population monitored from 2006 - 2021.") +
               theme_classic()
                   
sbull_plot1

Plot annual growth rate over population sizes - is this population density dependent?

sbull_plot2 <- sbull6 %>% 
               ggplot(aes(total_number_counted, growth_rate)) +
               geom_point() +
               geom_smooth(method = "lm", se = F) +
               labs(x = "Population Size (N)",
                    y = "Annual Growth Rate",
                    caption = " Figure 2. Model of density dependence of a Synthyris bullii population. There is a significant \n negative relationship between size and annual growth rate (slope = -0.008, p <0.05).") +
               geom_hline(yintercept = 0, lty = "dashed") +
               theme_classic()

sbull_plot2

There appears to be a negative relationship between population size and growth rate, is it significant?

sbull_model1 <- lm(growth_rate ~ total_number_counted ,data = sbull6)

sbull_summary1 <- summary(sbull_model1)

sbull_summary2 <- coef(sbull_summary1)

row.names(sbull_summary2) <- c("Intercept", "Population Size (N)")

knitr::kable(sbull_summary2,
             caption = "Table 1. Synthyris bullii population density dependence. Population size was used to predict annual growth rate. ", 
             col.names = c("Estimate", "Std. Error", "t-value", "p-value")) %>%
             kable_minimal(full_width = F, html_font = "Cambria", position = "left")
Table 1. Synthyris bullii population density dependence. Population size was used to predict annual growth rate.
Estimate Std. Error t-value p-value
Intercept 1.006495 0.4310899 2.334767 0.0362399
Population Size (N) -0.008831 0.0039811 -2.218220 0.0449658

Population Viability Analysis

Plot extinction risk over time across multiple extinction threshold sizes (Ne)

Ne = 5

logN <- log(sbull5$total_number_counted[-1]/sbull5$total_number_counted[-16])

sbull_ext_5 <- countCDFxt(mu=mean(logN), sig2=var(logN), nt=15, tq=15, Nc=53, Ne=5, tmax = 15, Nboot = 100, plot = T)
Figure 3. Extinction risk of Synthyris bullii population with extinction threshold set to 5 individuals.

Figure 3. Extinction risk of Synthyris bullii population with extinction threshold set to 5 individuals.

Ne = 10

sbull_ext_10 <- countCDFxt(mu=mean(logN), sig2=var(logN), nt=15, tq=15, Nc=53, Ne=10, tmax = 15, Nboot = 100, plot = T)
Figure 4. Extinction risk of Synthyris bullii population with extinction threshold set to 10 individuals.

Figure 4. Extinction risk of Synthyris bullii population with extinction threshold set to 10 individuals.

Ne = 20

sbull_ext_20 <- countCDFxt(mu=mean(logN), sig2=var(logN), nt=15, tq=15, Nc=53, Ne=20, tmax = 15, Nboot = 100, plot = T)
Figure 5. Extinction risk of Synthyris bullii population with extinction threshold set to 20 individuals.

Figure 5. Extinction risk of Synthyris bullii population with extinction threshold set to 20 individuals.

Conclusion

This population of S. bullii demonstrates clear negative density dependence. This could be explained by increased competition between individuals or possibly by limited habitat size. The large confidence interval makes interpreting extinction probabilities tricky, and it would be beneficial to incorporate more variables. As a whole, perisistence over the next 15 years looks promising, but continued monitoring will help inform the need to intervene and to understand potential ramifications of changes in climate.

References

Henry, Lionel, and Hadley Wickham. 2020. Purrr: Functional Programming Tools. https://CRAN.R-project.org/package=purrr.
Müller, Kirill, and Hadley Wickham. 2021. Tibble: Simple Data Frames. https://CRAN.R-project.org/package=tibble.
R Core Team. 2021. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.
Stubben, Chris J., and Brook G. Milligan. 2007. “Estimating and Analyzing Demographic Models Using the Popbio Package in r.” Journal of Statistical Software 22 (11).
Stubben, Chris, Brook Milligan, and Patrick Nantel. 2020. Popbio: Construction and Analysis of Matrix Population Models. https://CRAN.R-project.org/package=popbio.
Wickham, Hadley. 2016. Ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org.
———. 2019. Stringr: Simple, Consistent Wrappers for Common String Operations. https://CRAN.R-project.org/package=stringr.
———. 2021a. Forcats: Tools for Working with Categorical Variables (Factors). https://CRAN.R-project.org/package=forcats.
———. 2021b. Tidyr: Tidy Messy Data. https://CRAN.R-project.org/package=tidyr.
———. 2021c. Tidyverse: Easily Install and Load the Tidyverse. https://CRAN.R-project.org/package=tidyverse.
Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.
Wickham, Hadley, Winston Chang, Lionel Henry, Thomas Lin Pedersen, Kohske Takahashi, Claus Wilke, Kara Woo, Hiroaki Yutani, and Dewey Dunnington. 2021. Ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics. https://CRAN.R-project.org/package=ggplot2.
Wickham, Hadley, Romain François, Lionel Henry, and Kirill Müller. 2021. Dplyr: A Grammar of Data Manipulation. https://CRAN.R-project.org/package=dplyr.
Wickham, Hadley, Jim Hester, and Jennifer Bryan. 2021. Readr: Read Rectangular Text Data. https://CRAN.R-project.org/package=readr.