Rivulus-Genetic-Diversity.Rmd

---
title: "Rivulus-Genetic_Diversity"
author: "Anthony Snead"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Data Import and Formatting

Here we import our genetic data, assign individuals to population groupings, and format the population genetic data for downstream analysis

## Genetic Data Import

Here we import our population genetic data. The file includes columns that are not needed as well as a sister species indicated by the presence of the R22.Dup marker.

```{r, Data Import}
# Read and format the population genetic data for the spatial analysis and the genid object
Population.Data <- read.csv(paste(getwd(), "Data/rivulus_genetic_data.csv", sep = "/")) %>%
  # R22.Dup is a microsatellite marker species to K. hermaphroditus. Therefore, we remove any individuals with that marker as they are not our species of interest.
  dplyr::filter(is.na(R22.Dup) & is.na(R22.Dup.1)) %>%
  #get rid of the columns
  dplyr::select(-R22.Dup, -R22.Dup.1) %>%
  # remove individual with greater than 10% missing data
  dplyr::filter(rowSums((is.na(.[,6:ncol(.)])/32)*100) <= 10) %>%
  # replace the NA values
  replace(is.na(.), "?") %>%
  # gather the loci names into a new column
  tidyr::gather(., Loci, Val, -(1:5), na.rm = FALSE) %>%
  # remove the ".1" from the end of loci names
  dplyr::mutate(Loci = gsub("\\.[0-9]*$","", Loci)) %>%
  # group by all variables except the lovi values
  dplyr::group_by(across(1:6)) %>%
  # join loci values and separate by ":"
  dplyr::summarise(Val = paste(Val[!(is.na(Val)|Val=="")], collapse=":")) %>%
  # make loci names column headings again
  tidyr::spread(Loci, Val)%>%
  # ungroup
  ungroup() %>%
  # rename the ID column
  dplyr::rename(ID = Lines.Loci) %>%
  # get rid of the - male in ID columns
  dplyr::mutate(ID = gsub("- male","", ID)) %>%
  # arrange the data frame
  dplyr::arrange(Long, Lat)
```

## Create Population Groupings

We lump together populations that are within 10 km of each other because the coarsest Ocean General Circulation Model (OGCM) we use is ~9 km resolution. Here we create buffers and merge them together to get our groupsing for the downstream analysis.

```{r, Create Populations}

# make the population groups
Population.Groups <- Population.Data %>%
  # select coordinates
  dplyr::select(Lat, Long) %>%
  # get unique entries for buffers
  unique(.) %>% 
  # convert population data to simple feature.
  sf::st_as_sf(x = .,
               coords = c("Long", "Lat"),
               crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0") %>%
  sf::st_transform(., crs = "+proj=laea +lat_0=37.32 +lon_0=-113.04 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs") %>%
  # make 10km buffer
  sf::st_buffer(., dist = 10000) %>%
  # make valid
  sf::st_make_valid(.) %>%
  # combine areas with overlapping buffers
  sf::st_union(.) %>%
  # make the object valid
  sf::st_make_valid(.) %>%
  # cast to a polygon
  sf::st_cast(., "POLYGON") %>%
  # make a spatial
  sf::as_Spatial(.) 
```

## Assign Population Numbers

We use our population groupings to assign the population number back to the genetic data for downstream analysis.We then rename the populations base on where they are.

```{r, Population Assignment}

# assign the populations to the genetic data and filter out locations with less than 10 samples
Population.Genetic.Data <- Population.Data %>%
  # convert population data to simple feature.
  sf::st_as_sf(x = .,
               coords = c("Long", "Lat"),
               crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0") %>%
  # spatial join to polygon to assign the population ID
  sf::st_join(x = .,
              y = (sf::st_as_sf(Population.Groups,
                               coords = c("Longitude", "Latitude"),
                               crs = "+proj=laea +lat_0=37.32 +lon_0=-113.04 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs") %>%
                     sf::st_transform(., crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")  %>%
                     cbind(., Population_ID = 1:20)),
              join = st_intersects,
              left = FALSE) %>%
  # convert to data frame with geometry column
  as.data.frame(.) %>%
  # split geometry to lat and long
  tidyr::separate(., col = geometry, into = c("Longitude", "Latitude"), sep = ",") %>%
  # extract number 
  dplyr::mutate(Longitude = tidyr::extract_numeric(Longitude),
                Latitude = tidyr::extract_numeric(Latitude)) %>%
  # move the population Id column.
  dplyr::relocate(c(Population_ID), .after = ID) %>%
  # move columns around
  dplyr::relocate(c(Longitude, Latitude, REGION, Year), .before = ID) %>%
  # change numbers to acronyms
  dplyr::mutate(Population_ID = dplyr::case_when(Population_ID == 1 ~ "UH",
                                                 Population_ID == 2 ~ "RH",
                                                 Population_ID == 3 ~ "TC",
                                                 Population_ID == 4 ~ "LC",
                                                 Population_ID == 5 ~ "TA",
                                                 Population_ID == 6 ~ "NC",
                                                 Population_ID == 7 ~ "LK",
                                                 Population_ID == 8 ~ "UK",
                                                 Population_ID == 9 ~ "EI",
                                                 Population_ID == 10 ~ "SS",
                                                 Population_ID == 11 ~ "EG",
                                                 Population_ID == 12 ~ "SL",
                                                 Population_ID == 13 ~ "FL",
                                                 Population_ID == 14 ~ "LB",
                                                 Population_ID == 15 ~ "VR",
                                                 Population_ID == 16 ~ "BB",
                                                 Population_ID == 17 ~ "CC",
                                                 Population_ID == 18 ~ "TB",
                                                 Population_ID == 19 ~ "IR",
                                                 Population_ID == 20 ~ "NS")) %>%
  # group
  dplyr::group_by(Population_ID) %>%
  # remove populations with less than 10 individuals
  dplyr::filter(n() >= 10) %>%
  # ungroup
  dplyr::ungroup(.)

# The acroynmes represent where the samples are from. NS = New Symerna, EG = Everglades, LB = Lower Bogue, EI = Exuma Island, FL = Fort Lauderdale, IR = Indian River, UK = Upper Keys, SL = Saint Lucie, VR = Vero Beach, LC = Long Caye, TC = Twin Cayes, TA = Turneffe Atoll, SS = San Salvador, CC = Charlotte County, BB = Bunche Beach, LK = Lower Keys, UH = Utila Honduras, RH = Roatan Honduras, NC = Northern Caye, TB = Tampa Bay


# the base world data for maps
world <- rnaturalearth::ne_countries(scale = "large", returnclass = "sf")

# map of samples total
Sample.Map <- ggplot2::ggplot(data = world) +
   geom_sf(fill= "antiquewhite", 
            color = "black") +
  geom_point(data = Population.Genetic.Data, aes(x = Longitude, y = Latitude, fill = as.factor(Population_ID)),
             shape=21, size = 2, color = "black") +
  geom_text(data = Population.Genetic.Data, aes(x = Longitude, y = Latitude), label = Population.Genetic.Data$Population_ID, nudge_x = 0.2, nudge_y = 0.2, check_overlap = TRUE) +
    coord_sf(xlim = c(-90,-72),
             ylim = c(15,30))
```

## Convert Genetic Data

Here we convert our genetic data from into a genid object.
```{r, Convert Genetic Data}

# make genind object for population genetic analysis used in the Pairwise Comparison
Population.Genetic.Data.genid <- Population.Genetic.Data %>%
  dplyr::select(-(1:6)) %>% # get the columns we need
  adegenet::df2genind(X = .,
                      sep = ":",
                      ncod = NULL,
                      ind.names = Population.Genetic.Data$ID,
                      loc.names = NULL,
                      pop = Population.Genetic.Data$Population_ID,
                      NA.char= "?",
                      ploidy = 2,
                      type = "codom",
                      strata = Population.Genetic.Data[,3:6],
                      hierarchy = NULL)
```



## Genetic Diversity

```{r, unbiased shannon - zahl}

# unbiased shannon estimator from Konopiski 2020
ShannonGen <- function(gInd, estimator = NULL) {
  MLE = function(X) {
    X = X[X > 0]
    n = sum(X)
    - sum(X / n * log(X / n))
  }
  Z = function(X){
    X=X[X>0]
    Y=X[X>1]
    n=sum(X)
    -n*sum(X/n*log(X/n))-(n-1)/n*sum((n-X)*(-X/(n-1)*log(X/(n-1))) )-(n-1)/n*sum(-Y*(Y-1)/(n-1)*log((Y-1)/(n-1)))
  }
  CS = function(X) {
    x = X
    x = x[x > 0]
    n = sum(x)
    f1 = sum(x == 1)
    C_head = 1 - f1 / n
    a = -sum(C_head * (x / n) * log(C_head * (x / n)) / (1 - (1 - C_head *
                                                                (x / n)) ^ n))
    a
  }
  Ch = function(X) {
    x = X
    x = x[x > 0]
    n = sum(x)
    UE <- sum(x / n * (digamma(n) - digamma(x)))
    f1 <- sum(x == 1)
    f2 <- sum(x == 2)
    if (f1 > 0)
    {
      A <-
        1 - ifelse(f2 > 0, (n - 1) * f1 / ((n - 1) * f1 + 2 * f2), (n - 1) * f1 /
                     ((n - 1) * f1 + 2))
      B = sum(x == 1) / n * (1 - A) ^ (-n + 1) * (-log(A) - sum(sapply(1:(n - 1), function(k) {
        1 / k * (1 - A) ^ k})))
    }
    if (f1 == 0) {
      B = 0
    }
    if (f1 == 1 & f2 == 0) {
      B = 0
    }
    UE + B
  }
  inputs <- lapply(adegenet::seppop(gInd),
                   function (i)
                     sapply(adegenet::seploc(i),
                            function(j) colSums(j@tab,na.rm = TRUE)))
  out <- list()
  if (is.null(estimator))
    estimator <- "z"
  
  if ("sh" %in% estimator) {
    out[["Shannon_1949"]] <-
      as.data.frame(lapply(inputs, function (i)
        sapply(i, MLE)))
  }
  if ("z" %in% estimator) {
    out[["Zahl_1977"]] <-
      as.data.frame(lapply(inputs, function (i)
        sapply(i, Z)))
  }
  if ("cs" %in% estimator) {
    out[["Chao_Shen_2003"]] <-
      as.data.frame(lapply(inputs, function (i)
        sapply(i, CS)))
  }
  if ("ch" %in% estimator) {
    out[["Chao_et_al_2013"]] <-
      as.data.frame(lapply(inputs, function (i)
        sapply(i, Ch)))
  }
  return(out)
}

Zahl <- ShannonGen(Population.Genetic.Data.genid)$Zahl_1977 %>%
  # get mean
  dplyr::summarise_all(., mean) %>%
  # pivot
  tidyr::pivot_longer(., cols = 1:ncol(.), names_to = "Pop", values_to = "Z")

```
```{r, Genetic Diversity}
#calculate genotypic diversity metrics
Population.Genetic.Diversity <- poppr::poppr(Population.Genetic.Data.genid, total = TRUE, sample = 999, cutoff = 0.1)

# get rarefied genetic diversity
Population.Genetic.Diversity.Rare <- poppr::diversity_ci(Population.Genetic.Data.genid,
                                                         total = TRUE, rarefy = TRUE,
                                                         n.rare = 10)

# combine genotype rarefied measures with standard genetic diversity measures
Population.Genetic.Diversity.Table <- data.frame(
  Pop = names(apply((t(sapply(adegenet::seppop(Population.Genetic.Data.genid), function(ls) adegenet::summary(ls)$Hobs))),
                    MARGIN = 1, FUN = mean)),
  Hobs = apply((t(sapply(adegenet::seppop(Population.Genetic.Data.genid), function(ls) adegenet::summary(ls)$Hobs))),
               MARGIN = 1, FUN = mean),
  Hexp = apply((t(sapply(adegenet::seppop(Population.Genetic.Data.genid), function(ls) adegenet::summary(ls)$Hexp))),
               MARGIN = 1, FUN = mean),
  Ar = PopGenReport::allel.rich(Population.Genetic.Data.genid, min.alleles = NULL)$mean.richness) %>%
  dplyr::left_join(., (as.data.frame(Population.Genetic.Diversity.Rare$est) %>%
                     tibble::rownames_to_column(., var = "Pop"))) %>%
  dplyr::left_join(., (Population.Genetic.Diversity %>%
                         dplyr::select(Pop, eMLG))) %>%
  dplyr::left_join(., Zahl) %>%
  dplyr::relocate(., "eMLG", .before = "Hobs") %>% 
  dplyr::select(-H, -lambda, -E.5)
```

```{r, sample size versus diversity}
# create table for correlation test
Diversity_N <- Population.Genetic.Diversity.Table %>%
  dplyr::left_join(., (Population.Genetic.Diversity %>%
                         dplyr::select(Pop, N))) %>%
  dplyr::select(N, eMLG, Z, G, Hexp, Hobs, Ar)

#eMLG
cor.test(Diversity_N$N, Diversity_N$eMLG)

#Z
cor.test(Diversity_N$N, Diversity_N$Z)

#G
cor.test(Diversity_N$N, Diversity_N$G)

#Hexp
cor.test(Diversity_N$N, Diversity_N$Hexp)

#Hobs
cor.test(Diversity_N$N, Diversity_N$Hobs)

#Ar
cor.test(Diversity_N$N, Diversity_N$Ar)

```

# Abiotic Variables
## CMS

We use the connectivity-modeling-system (CMS) to estimate ocean current mediated connectivity by simulating an ocean current environment and releasing particles that can represent larva, eggs ,or other propagules. CMS requires ocean current data (nest files), two parameter files, and a release file. It generates either netCDF or ASCII files with the location of each particle at the designated time step. Here we generate the release files and process the output.

### Generate Release Files

The release file is a tab delimited text file with 9 columns (polygon, longitude, latitude, depth, number, year, month, day, and second). Polygon is used with the seascape module if particles are allowed to settle. Even if the particles are not allowed to settle, the polygon column is still required and must be an integer. Longitude, latitude, and depth columns represent the location of the particle release location both geographically and vertically. The number is the number of particles released at that specific combination of location and time. If there is no horizontal diffusivity coefficient, all the particles released at the same location and time will follow the same trajectory. The year, month, day, and second columns dictate the time of release. Particles cannot be released over land, and if they are, CMS will not notify you of the issue. CMS dictates a particle is on land if two of the closest eight points are dictated as land my the nest files.

To get an accurate representation of ocean current connectivity, 1000 locations were selected and 5 particles released form each location at a random point every day for 10 years. Because locations could not be over land, we had to troubleshoot release locations using the example HYCOM files from the specific experiments used in our simulations and additional buffer zones. Therefore, the following code requires importing the HYCOM data, sampling the points, and creating the release files.
#### Release Points

##### HYCOM Data

CMS enables users to use different nest files that are different in resolution. We used two HYCOM experiments (GOM 50.1 & GLBu 0.08). The GOM 50.1 is only the Gulf of Mexico but is available at a 1/25° resolution, while the GLBu 0.08 is global but only at a 1/12° resolution. Here we import sample netCDF files invert them so that the land mask is not equal to one, convert them to polygons, and add a 5 km buffer to them. We then use them later to filter out locations that CMS would dictate as land. 


```{r, HYCOM data}

# get netcdf, convert to polygon, and calculate the buffer
GOM <- raster::raster(paste(getwd(), "Data/Sample_HYCOM/hycom_gomu_501_2012010100_t000.nc", sep = "/"), varname = "water_v") %>%
  # crop the the study extent to make sure that cells outside the domain are not retained inappropriately 
  raster::crop(., raster::extent(c(xmin = -98.1 , xmax = -76.5 , ymin = 18 , ymax = 32))) %>%
  # invert the raster so that the land mask is the only cells with data
  raster::calc(., fun = function(x){
    if(is.na(x[1])){
      return(1)}
    else{ return(NA)}}) %>%
  # convert to polygon
  raster::rasterToPolygons(., dissolve = TRUE) %>%
  # convert to sf object
  sf::st_as_sf(x = .,
               coords = c("Long", "Lat"),
               crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0") %>%
  # project for the buffer
  sf::st_transform(., crs = "+proj=laea +lat_0=37.32 +lon_0=-113.04 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs") %>%
  # calculate the 5 km buffer
  sf::st_buffer(., dist = 5000) %>%
  # unproject because CMS take lat/long
  sf::st_transform(., crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
  
# get netcdf, convert to polygon, and calculate the buffer
GBLu <- raster::raster(paste(getwd(), "Data/Sample_HYCOM/hycom_GLBu0.08_191_1995080100_t000.nc", sep = "/"), varname = "water_v") %>%
  # crop the the study extent to make sure that cells outside the domain are not retained inappropriately 
  raster::crop(., raster::extent(c(xmin = -91 , xmax = -76 , ymin = 18 , ymax = 32))) %>%
  # invert the raster so that the land mask is the only cells with data
  raster::calc(., fun = function(x){
    if(is.na(x[1])){
      return(1)}
    else{ return(NA)}}) %>%
  # convert to polygon
  raster::rasterToPolygons(., dissolve = TRUE) %>%
  # convert to sf object
  sf::st_as_sf(x = .,
               coords = c("Long", "Lat"),
               crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0") %>%
  # project for the buffer
  sf::st_transform(., crs = "+proj=laea +lat_0=37.32 +lon_0=-113.04 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs") %>%
  # calculate the 5 km buffer
  sf::st_buffer(., dist = 5000) %>%
  # unproject because CMS take lat/long
  sf::st_transform(., crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
```

##### Sample Locations

Here we use the population polygons to sample release locations for CMS. This step required troubleshooting to determine the appropriate buffer to add to the population group polygons to ensure that the particles were not considered over land and not overlapping with an adjacent population.

```{r, Release Locations}

# Get a data frame of points along the edge each polygon group for release points
Population.Points <- Population.Groups %>%
  # convert to sf
  sf::st_as_sf(x = .,
               coords = c("X", "Y"),
               crs = "+proj=laea +lat_0=37.32 +lon_0=-113.04 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs") %>%
  # make 2km buffer to push the points farther from the coast
  sf::st_buffer(., dist = 2000) %>%
  # make valid
  sf::st_make_valid(.) %>%
  # combine areas with overlapping buffers
  sf::st_union(.) %>%
  # make the object valid
  sf::st_make_valid(.) %>%
  # cast to a polygon
  sf::st_cast(., "POLYGON") %>%
  # transform to lat/long cms takes lat/long
  sf::st_transform(., crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0") %>%
  # get points along the edge
  sf::st_segmentize(., dfMaxLength = 0.1) %>%
  # convert to coordinates
  sf::st_coordinates(.) %>%
  # make a data frame
  as.data.frame(.) %>%
  # select important columns
  dplyr::select(X, Y, L2) %>%
  # rename the population ID column
  dplyr::rename(Population_ID = L2) %>%
  # reduce to third decimal
  dplyr::mutate(X = format(round(X, 3), nsmall = 3),
                Y = format(round(Y, 3), nsmall = 3)) %>%
  # get unique points
  dplyr::distinct(.) %>%
  # make sf object
  sf::st_as_sf(.,
               coords = c("X", "Y"),
               crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0") %>%
  # add location columns and columns that indicate if a point is within the land mask for each HYCOM experiment 
  dplyr::mutate(GOM = lengths(sf::st_within(., GOM)),
                GBLu = lengths(sf::st_within(., GBLu)),
                Longitude = sf::st_coordinates(.)[,1],
                Latitude = sf::st_coordinates(.)[,2]) %>%
  # convert to data frame
  as.data.frame(.) %>%
  # filter out locations by the finest resolution nest that over over land
  dplyr::filter((Latitude > 18 & Latitude < 32 & Longitude > -98.1 & Longitude < -76.5 & GOM == 0) |
                  (Latitude < 18 & GBLu == 0) | (Longitude >= -76.5 & GBLu == 0)) %>%
  # select important columns
  dplyr::select(Longitude, Latitude, Population_ID)
```

##### Create Release files

Here we use the locations to create the release files. We first create a master release file data frame that includes all the years.

```{r, Release Data Frame}

# set the seed for reproducability 
set.seed(42)

# make file of release points that are over water
Release.Points <- Population.Points %>%
  # convert the longitude to same format as output for join
  dplyr::mutate(Longitude = 360+Longitude) %>%
  # get rid of population that did not have enough samples
  dplyr::filter(Population_ID !=2 & Population_ID != 15 & Population_ID != 16) %>%
  # group
  dplyr::group_by(Population_ID) %>%
  # sample 1000 locations randomly
  dplyr::sample_n(1000, replace = TRUE) %>%
  # replicate each row by the total number of days from the start of the simulation to the end.
  dplyr::slice(rep(row_number(), (10*12*31))) %>%
  # add columns for CMS and fill Year column
  dplyr::mutate(Depth = 0,
                Number = 10,
                Year = rep(2000:2009, each = (1000*12*31))) %>%
  # group
  dplyr::group_by(Population_ID, Year) %>%
  # fill Month column
  dplyr::mutate(Month = rep(1:12, each = (1000*31))) %>%
  # group
  dplyr::group_by(Population_ID, Year, Month) %>%
  # fill day column
  dplyr::mutate(Day = rep(1:31, each = 1000)) %>%
  # ungroup
  dplyr::ungroup() %>%
  # fill second with a random time
  dplyr::mutate(Second = round(runif(nrow(.), 0, 86400))) %>%
  # filter out all the days that do not exist and Feb 29th because the CMS software ignores it anyways
  dplyr::filter((Month %in% c(1,3,5,7,8,10,12) & Day <= 31) | (Month %in% c(4,6,9,11) & Day <= 30) | (Month == 2 & Day <= 28)) %>%
  # make sure columns are in the correct order for CMS
  dplyr::relocate(Population_ID, .before = Longitude)
```

Even with an HPC, it was more effective to run each year seperatley. Therefore, we write a release file for each year.

```{r, write files}

# write release file for 2000
write.table(Release.Points %>%
               dplyr::filter(Year == 2000),file = paste(getwd(), "CMS/Input/Release_Points_2000.txt", sep = "/"),quote = FALSE, sep = '\t', col.names = FALSE, row.names = FALSE)

# write release file for 2001
write.table(Release.Points %>%
               dplyr::filter(Year == 2001),file = paste(getwd(), "CMS/Input/Release_Points_2001.txt", sep = "/"),quote = FALSE, sep = '\t', col.names = FALSE, row.names = FALSE)

# write release file for 2002
write.table(Release.Points %>%
               dplyr::filter(Year == 2002),file = paste(getwd(), "CMS/Input/Release_Points_2002.txt", sep = "/"),quote = FALSE, sep = '\t', col.names = FALSE, row.names = FALSE)

# write release file for 2003
write.table(Release.Points %>%
               dplyr::filter(Year == 2003),file = paste(getwd(), "CMS/Input/Release_Points_2003.txt", sep = "/"),quote = FALSE, sep = '\t', col.names = FALSE, row.names = FALSE)

# write release file for 2004
write.table(Release.Points %>%
               dplyr::filter(Year == 2004),file = paste(getwd(), "CMS/Input/Release_Points_2004.txt", sep = "/"),quote = FALSE, sep = '\t', col.names = FALSE, row.names = FALSE)

# write release file for 2005
write.table(Release.Points %>%
               dplyr::filter(Year == 2005),file = paste(getwd(), "CMS/Input/Release_Points_2005.txt", sep = "/"),quote = FALSE, sep = '\t', col.names = FALSE, row.names = FALSE)

# write release file for 2006
write.table(Release.Points %>%
               dplyr::filter(Year == 2006),file = paste(getwd(), "CMS/Input/Release_Points_2006.txt", sep = "/"),quote = FALSE, sep = '\t', col.names = FALSE, row.names = FALSE)

# write release file for 2007
write.table(Release.Points %>%
               dplyr::filter(Year == 2007),file = paste(getwd(), "CMS/Input/Release_Points_2007.txt", sep = "/"),quote = FALSE, sep = '\t', col.names = FALSE, row.names = FALSE)

# write release file for 2008
write.table(Release.Points %>%
               dplyr::filter(Year == 2008),file = paste(getwd(), "CMS/Input/Release_Points_2008.txt", sep = "/"),quote = FALSE, sep = '\t', col.names = FALSE, row.names = FALSE)

# write release file for 2009
write.table(Release.Points %>%
               dplyr::filter(Year == 2009),file = paste(getwd(), "CMS/Input/Release_Points_2009.txt", sep = "/"),quote = FALSE, sep = '\t', col.names = FALSE, row.names = FALSE)
```

### Post Processesing 
```{r, CMS Post-Processing Functions}

# the function processes the CMS output files
CMS.Processing.Function <- function(pathlist, polygons, timestep, releasepoints){
  
  library(dplyr)
  library(sf)
  library(sfheaders)
  library(tidyverse)
  
  # read in the population polygon and convert to sf
  Populations <- sf::st_as_sf(polygons,
                               coords = c("X", "Y"),
                               crs = "+proj=laea +lat_0=37.32 +lon_0=-113.04 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs") %>%
    # transform
    sf::st_transform(., crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0") %>%
    # add population ID
    cbind(., Population_ID = (1:20))
  
  releasepoints <- releasepoints %>%
    # get only important columns
    dplyr::select(Population_ID, Longitude, Latitude) %>%
    # get only distinct values
    distinct(.) %>%
    # convert lat and long to character for join
    dplyr::mutate( Longitude = as.character(Longitude),
                   Latitude = as.character(Latitude))
  
  #a function to assign the sources
  Assignment.Fun <- function(path, polygon, points){
    
  # read in the trajectory file
  CMS.Output <- readr::read_table(path, 
                                  col_names = c("Location", "Particle", "Time", "Longitude", "Latitude", "Depth", "Distance", "Exit_Code", "Release_Date"))  %>%
    # remove the columns we do not need
    dplyr::select(-Depth, -Exit_Code, -Release_Date)
  
  # assign the sources
  Assigned.Release <- CMS.Output %>%
    # filter to retain only the first release point
    dplyr::filter(Time == 0) %>%
    # convert the Latitude and longitude to character vectors for the left join
    dplyr::mutate(Longitude = as.character(Longitude),
                  Latitude = as.character(Latitude)) %>%
    # join with release locations
    dplyr::left_join(., points) %>%
    # select the variables we need
    dplyr::select(Location, Population_ID, Particle) %>%
    # relabel the source column
    dplyr::rename(Source = Population_ID) %>%
    # convert location and particle to characters for joining
    dplyr::mutate(Location = as.character(Location),
                  Particle = as.character(Particle))
  
  # Assign the sink and join to the source
  Source.Sink <- CMS.Output %>%
    # get particle that have moved
    dplyr::filter(Time != 0) %>%
    # convert to sf object
    sf::st_as_sf(.,
               coords = c("Longitude", "Latitude"),
               crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0") %>%
    # join it with polygons
    sf::st_join(x = .,
                y = polygon,
              join = st_intersects,
              left = TRUE) %>%
    # convert to data frame
    sfheaders::sf_to_df(., fill = TRUE) %>%
    # select important columns
    dplyr::select(-sfg_id, -point_id) %>%
    # rename columns
    dplyr::rename(Longitude = x,
                 Latitude = y,
                 Sink = Population_ID) %>%
    # convert location and particle to character for left join
    dplyr::mutate(Location = as.character(Location),
                  Particle = as.character(Particle)) %>%
    # join with the sources
    dplyr::left_join(., Assigned.Release) %>%
    # get rid of rows that do not have a sink
    dplyr::filter(!is.na(Sink))
  
  return(Source.Sink)
  }
  
  # calculate connectivity values
  Connectivity.Values <- lapply(pathlist, FUN = Assignment.Fun, polygon = Populations, points = releasepoints) %>% # apply function to all files
    # bind the lists together
    bind_rows(.) %>%
    # group by unique particles for each source sink relationship
    dplyr::group_by(Source, Sink, Location, Particle) %>%
    # arrange by time
    dplyr::arrange(Time, by.group = TRUE) %>%
    # calculate the difference between times
    dplyr::mutate(Time_diff = Time - lag(Time)) %>%
    # mutate to get connectivity values
    dplyr::mutate(Time_diff = replace(Time_diff, Time_diff != timestep, timestep), # retain only time steps that are consecutive and if they are not give it a value of the timestep
                  Time = (sum(Time_diff, na.rm = TRUE)), # get total amount of time in the polygon
                  Distance = min(Distance)) %>% # get min distance traveled
    # get distinction combinations to not count particles twice
    dplyr::distinct(Source, Sink, Location, Particle, Distance, Time, .keep_all = TRUE) %>%
    # group by source sink
    dplyr::group_by(Source, Sink) %>%
    # summarise 
    dplyr::summarise(n = dplyr::n_distinct(Location, Particle), # number of particles
                     Distance = mean(Distance), # average distance
                     Time = mean(Time, na.rm = TRUE)) %>% # mean time within polygon
    # ungroup
    dplyr::ungroup(.)

return(Connectivity.Values)}

```

```{r, CMS Post-Processing}
# processes the ocean current data
Ocean.Connectivity <- CMS.Processing.Function(pathlist = list.files(paste(getwd(), "CMS/Output/", sep ="/"), recursive  = TRUE, full.names = TRUE), polygons = Population.Groups, releasepoints = Release.Points, timestep = 1800) %>%
  # convert to characters
  dplyr::mutate(across(1:2, ~ as.character(.))) %>% 
  # rename populations 
  dplyr::mutate(across(1:2, ~ dplyr::case_when(. == "1" ~ "UH",
                                               . == "2" ~ "RH",
                                               . == "3" ~ "TC",
                                               . == "4" ~ "LC",
                                               . == "5" ~ "TA",
                                               . == "6" ~ "NC",
                                               . == "7" ~ "LK",
                                               . == "8" ~ "UK",
                                               . == "9" ~ "EI",
                                               . == "10" ~ "SS",
                                               . == "11" ~ "EG",
                                               . == "12" ~ "SL",
                                               . == "13" ~ "FL",
                                               . == "14" ~ "LB",
                                               . == "15" ~ "VR",
                                               . == "16" ~ "BB",
                                               . == "17" ~ "CC",
                                               . == "18" ~ "TB",
                                               . == "19" ~ "IR",
                                               . == "20" ~ "NS"))) %>%
  # filter out populations not analyzed
  dplyr::filter(Source != "RH" & Source != "VR" & Source != "BB" & Sink != "RH" & Sink != "VR" & Sink != "BB")

Ocean.Connectivity.NS <- Ocean.Connectivity %>%
   dplyr::left_join(.,
                    Population.Genetic.Diversity,
                    by = c("Source" = "Pop")) %>%
  dplyr::mutate(Time = (n*Time)) %>%
  dplyr::mutate(Time = Time/sum(Time))

Ocean.Connectivity.Standard <- Ocean.Connectivity %>%
   dplyr::left_join(.,
                    Population.Genetic.Diversity,
                    by = c("Source" = "Pop")) %>%
  dplyr::mutate(Time = (n*Time)*(eMLG/12)) %>%
  dplyr::mutate(Time = Time/sum(Time))
  
```

### Centrality measures

```{r, Centrality Function}
# function makes graphs and estimates centrality measures from the graph
Centrality_Function <- function(Data, Variable, Distance){
  
  # load libraries
  library(tidyverse)
  library(igraph)
  
  # create graph
  Graph.Close <- Data %>%
    # rename to what igraph expects
    dplyr::rename(From = Source,
                  To = Sink,
                  weight = !!Variable) %>%
    # select the important columns
    dplyr::select(From, To, weight) %>%
    # filter out weights less than 0
    dplyr::filter(weight > 0) %>%
    # make directed graph
    igraph::graph_from_data_frame(., directed = TRUE, vertices = NULL)
  
  # create graph
  Graph.Strength <- Data %>%
    # rename to what igraph expects
    dplyr::rename(From = Source,
                  To = Sink,
                  weight = !!Variable) %>%
    # select the important columns
    dplyr::select(From, To, weight) %>%
    # filter out weights less than 0
    dplyr::filter(weight > 0) %>%
    # make directed graph
    igraph::graph_from_data_frame(., directed = TRUE, vertices = NULL)
  
  # if it is not a distance
  if(Distance == FALSE){
    
    # log the weights if not 0 and multiply by negative 1
    E(Graph.Close)$weight <- -1*ifelse(E(Graph.Close)$weight != 0, log(E(Graph.Close)$weight), 0)
  }
  
  # make results data frame fore each centrality measure
  Results <- tibble::enframe((igraph::closeness(Graph.Close, vids = igraph::V(Graph.Close),
                                                 mode = "in", normalized = TRUE)),
                                      name = "Population_ID", value = paste(Variable,"Close.In", sep ="_")) %>%
  dplyr::left_join(., tibble::enframe((igraph::strength(Graph.Strength, vids = igraph::V(Graph.Strength), mode = "in", loops = FALSE)),
                                      name = "Population_ID", value = paste(Variable,"Strength.In", sep ="_")))
  
  #return results
  return(Results)
}
```

```{r, estimate centrality}

# get standardized OC centrality measures
OC.Centrality.NS<- Ocean.Connectivity.NS %>%
  Centrality_Function(.,
                                     Variable = "Time", Distance = FALSE)

# get standardized OC centrality measures
OC.Centrality.Standard <- Ocean.Connectivity.Standard %>%
  Centrality_Function(.,
                                     Variable = "Time", Distance = FALSE)
```
## Distance

Rivulus are unlikely to travel over land. Therefore, we use distance over water as our distance measure.To do so, we first get centroids of all the populations, a land mask, and create a cost transition matrix. We then wrote a function that iterates across all points to get the shortest distance avoiding the land mask.

### Centroids

We get the centroid of each population.

```{r, Centroids}
# get centroids for calculating the distance over water
Centroids <- sf::st_as_sf(Population.Groups,
                                   coords = c("X", "Y"),
                                   crs = "+proj=laea +lat_0=37.32 +lon_0=-113.04 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs") %>%
  # get centroids
  sf::st_centroid(.) %>%
  # transform
  sf::st_transform(., crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0") %>%
  # get coordinates
  sf::st_coordinates(.) 
```

### Land.mask

We get a land mask from the world data we imported for maps.

```{r, Land mask}
# make land mask
Land.mask <- world %>%
  # crop
  sf::st_crop(., xmin = -100, xmax = -65, ymin = 0, ymax = 32) %>%
  # convert to raster
  raster::rasterize(., (raster::raster(nrow = 1000, ncol = 1000, ext=extent(-100, -65, 0, 32))))
```

### Transition object

We create a transition object from the land mask.

```{r, transition object}
# make transition object
Water.Cost <- Land.mask %>%
  # make water 1 and land 999
  raster::calc(., fun = function(x){
    if(is.na(x[1])){
      return(1)}
    else{ return(999)}}) %>%
  # make into transition
  gdistance::transition(., function(x) 1/mean(x), 16) %>%
  # correct
  gdistance::geoCorrection(., type = "c", scl = FALSE)

```

### Ocean distance function

We write a function to iterate across all points and get the distance over water.

```{r, Ocean distance function}
# function to get distances over water
Ocean.Distance.Fun <- function(points, cost_raster, landmask){
  
  library(gdistance)
  library(gissr)
  library(tidyverse)
  library(raster)
  
  # empty list
  Distance <- list()
  
  # loop to get distances
  for(i in 1:nrow(points)){
    # get cost raster for the point
    Cost <- gdistance::accCost(cost_raster, points[i,]) %>%
      # mask land
      raster::mask(., landmask, inverse = TRUE)
    
    # get shortest path and calcualte distance
    Distance[[i]] <- gdistance::shortestPath(cost_raster, points[i,], points[-i,], output = "SpatialLines") %>%
      # get distance
      gissr::sp_length(.) %>%
      # make into data frame
      data.frame(Pop1 = i,
               Pop2 = (1:nrow(points))[-i],
               Distance = .)
  }
  
  return(dplyr::bind_rows(Distance))
}

```

### Ocean Distance

We get the distance over marine habitat for all populations.

```{r, Ocean distances}
Ocean.Distance.df <- Ocean.Distance.Fun(Centroids, Water.Cost, Land.mask) %>%
  # make population acronyms to match the genetic data
  dplyr::mutate(across(1:2, ~ dplyr::case_when(. == 1 ~ "UH",
                                               . == 2 ~ "RH",
                                               . == 3 ~ "TC",
                                               . == 4 ~ "LC",
                                               . == 5 ~ "TA",
                                               . == 6 ~ "NC",
                                               . == 7 ~ "LK",
                                               . == 8 ~ "UK",
                                               . == 9 ~ "EI",
                                               . == 10 ~ "SS",
                                               . == 11 ~ "EG",
                                               . == 12 ~ "SL",
                                               . == 13 ~ "FL",
                                               . == 14 ~ "LB",
                                               . == 15 ~ "VR",
                                               . == 16 ~ "BB",
                                               . == 17 ~ "CC",
                                               . == 18 ~ "TB",
                                               . == 19 ~ "IR",
                                               . == 20 ~ "NS")))

```
## Landscape Metrics

Here we calculate landscape metrics for each population. These metrics are related to landscape configuration including metrics like habitat area and connectedness. We first start by importing and formatting the rasters before calculating the metrics. 

### Raster Formatting

We start by creating a function that will clip, format, and project our rasters.

```{r, Raster Formatting Function}
# function to crop and reproject the land cover data for each population where our habitat is 1 and all other habitat is 2 
Raster.Processing.Fun <- function(path, Clipping.DF, min.code, max.code){
  
  # get libraries we need
  library(raster)
  library(tidyverse)
  library(sf)

  # get list of raster files
  File.List <- list.files(path = path, full.names = TRUE)
  
  # make empty list to hold the rasters
  Cropped.Raster.List <- list()

  # iterate through all the populations
  for(i in 1:nrow(Clipping.DF)) {
    
    # iterate throug all the raster
    for( z in 1:length(File.List)){
      
      # import raster
      Raster <- raster::raster(File.List[[z]])
      
      # test if the extent of the population is within the extent of the raster
      if(sf::st_bbox(Clipping.DF$geometry[i])[[1]] >= sf::st_bbox(Raster)[[1]] &
         sf::st_bbox(Clipping.DF$geometry[i])[[2]] >= sf::st_bbox(Raster)[[2]] &
         sf::st_bbox(Clipping.DF$geometry[i])[[3]] <= sf::st_bbox(Raster)[[3]] &
         sf::st_bbox(Clipping.DF$geometry[i])[[4]] <= sf::st_bbox(Raster)[[4]]) {
        
        # get the cropped raster
        Cropped.Raster.List[i] <- raster::mask(raster::crop(Raster, 
                                                            raster::extent(Clipping.DF[i,])),
                                               Clipping.DF[i,]) %>%
          # project the raster to meters
          raster::projectRaster(., crs = Clipping.DF$Projection[[i]]) %>%
          # recategorize all the cells if they are between two values
          raster::calc(., fun = function(x){
            ifelse(x > min.code & x < max.code, 1, 2)})
      }
    }
    }
  
  # name the rasters by the populations used to crop them
  names(Cropped.Raster.List) <- Clipping.DF$Population_ID
  
  # return the rasters
  return(Cropped.Raster.List)
}
```


We create a data frame that has the population buffers, their names, and the projection that should be used to convert them to UTM.

```{r, Clipping Data frame}
# make sf object for easier graphing and to combine with population labels
Labeled.Pop.Groups <- sf::st_as_sf(Population.Groups, group = TRUE) %>%
  # unproject
  sf::st_transform(., crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0") %>%
  #combine with population labels
  cbind(., Population_ID = c("UH", "RH", "TC", "LC", "TA", "NC", "LK", "UK", "EI", "SS", "EG", "SL", "FL", "LB", "VR", "BB", "CC", "TB", "IR", "NS")) %>%
  # combine with the projection for each population
  cbind(., Projection = c( "+proj=utm +zone=16 +ellps=clrk66 +units=m +no_defs",
                           "+proj=utm +zone=16 +ellps=clrk66 +units=m +no_defs",
                           "+proj=utm +zone=16 +ellps=clrk66 +units=m +no_defs",
                           "+proj=utm +zone=16 +ellps=clrk66 +units=m +no_defs",
                           "+proj=utm +zone=16 +ellps=clrk66 +units=m +no_defs",
                           "+proj=utm +zone=16 +ellps=clrk66 +units=m +no_defs",
                           "+proj=utm +zone=17 +ellps=clrk66 +units=m +no_defs",
                           "+proj=utm +zone=17 +ellps=clrk66 +units=m +no_defs",
                           "+proj=utm +zone=18 +ellps=clrk66 +units=m +no_defs",
                           "+proj=utm +zone=18 +ellps=clrk66 +units=m +no_defs",
                           "+proj=utm +zone=17 +ellps=clrk66 +units=m +no_defs",
                           "+proj=utm +zone=17 +ellps=clrk66 +units=m +no_defs",
                           "+proj=utm +zone=17 +ellps=clrk66 +units=m +no_defs",
                           "+proj=utm +zone=18 +ellps=clrk66 +units=m +no_defs",
                           "+proj=utm +zone=17 +ellps=clrk66 +units=m +no_defs",
                           "+proj=utm +zone=17 +ellps=clrk66 +units=m +no_defs",
                           "+proj=utm +zone=17 +ellps=clrk66 +units=m +no_defs",
                           "+proj=utm +zone=17 +ellps=clrk66 +units=m +no_defs",
                           "+proj=utm +zone=17 +ellps=clrk66 +units=m +no_defs",
                           "+proj=utm +zone=17 +ellps=clrk66 +units=m +no_defs"))

```

We actually format the rasters with the data frame and path to the rasters
```{r, Format Rasters}
# crop and project the rasters
Cropped.Rasters <- Raster.Processing.Fun(path = paste(getwd(),"Data/GIS/", sep = "/"),
                                         Clipping.DF = Labeled.Pop.Groups,
                                         min.code = 9,
                                         max.code = 15)
```

### Calcualte Landscape Metrics

Here we calculate the landscape metrics.

```{r, Landscape Metrics}
# calculate the landscape metrics for each population
Landscape.Metrics <- lapply(Cropped.Rasters, FUN = function(x) {
  x %>%
    # calculate the landscape metrics
    landscapemetrics::calculate_lsm(., level = "class", directions = 8, count_boundary = FALSE,
                                    consider_boundary = FALSE, edge_depth = 1, cell_center = FALSE,
                                    neighbourhood = 8, what = c("lsm_c_ca", 
                                                                "lsm_c_cohesion",
                                                                "lsm_c_ed",
                                                                "lsm_c_np")) %>% 
    # only keep it for our class
    dplyr::filter(class == 1)
}) %>%
  # join all the data frame together and keep the population id. 
  dplyr::bind_rows(., .id = "Population_ID") %>%
  # keep only the columns we want
  dplyr::select(Population_ID, metric, value)

# run correlations between all landscape metrics
Landscape.Corr <- Landscape.Metrics %>%
  # convert from long to wide 
  tidyr::spread(metric, value) %>%
  # select only numerics
  dplyr::select(where(is.numeric)) %>%
  # run correlation
  cor(.)
```

# Analysis

## Data Formating

```{r, Data Formating}
Modeling.Data.NS <- Population.Genetic.Diversity.Table %>%
  dplyr::rename(.,
                "Population_ID" = "Pop") %>%
  # join with diversity metrics
  dplyr::left_join(., (Landscape.Metrics %>%
                         tidyr::spread(metric, value))) %>%
  dplyr::left_join(., OC.Centrality.NS) %>%
  dplyr::mutate(across(8:ncol(.), ~ (.x - mean(.x))/sd(.x)))

# standardized modeing data
Modeling.Data.S <- Population.Genetic.Diversity.Table %>%
  dplyr::rename(.,
                "Population_ID" = "Pop") %>%
  # join with diversity metrics
  dplyr::left_join(., (Landscape.Metrics %>%
                         tidyr::spread(metric, value))) %>%
  dplyr::left_join(., OC.Centrality.Standard) %>%
  dplyr::mutate(across(8:ncol(.), ~ (.x - mean(.x))/sd(.x)))

```

```{r, VIF}
#VIF stepwise with a threshold of 5
usdm::vifstep(as.data.frame(Modeling.Data.NS[, 8:ncol(Modeling.Data.NS)]),  th = 5)
```

```{r, Remove Colinear Variables}
Modeling.Data.Filtered.NS <- Modeling.Data.NS %>%
  dplyr::select(1:7, ca, np, Time_Close.In, Time_Strength.In) %>%
  dplyr::mutate(G = (G)^(3),
                Hexp = (Hexp)^3, # power 3 is best
                Ar = (Ar), # regular is best
                Hobs = Hobs, # regular is best
                Z = (Z)^2)

# standardized
Modeling.Data.Filtered.S <- Modeling.Data.S %>%
  dplyr::select(1:7, ca, np, Time_Close.In, Time_Strength.In) %>%
  dplyr::mutate(G = (G)^(3),
                Hexp = (Hexp)^3, # power 3 is best
                Ar = (Ar), # regular is best
                Hobs = Hobs, # regular is best
                Z = (Z)^2)


```
## Non-Clone Corrected

### Modeling

#### Stoddard and Taylor
```{r, G Formula creation}

Z.Form <- names((Modeling.Data.Filtered.NS %>%
                   dplyr::select(8:ncol(.)))) %>% 
  paste(., collapse = " + ") %>%
  paste("Z ~", ., sep = " ") %>%
  paste(., "ca:Time_Close.In", "ca:Time_Strength.In",
        "np:Time_Close.In", "np:Time_Close.In", sep = " + ")

```

```{r, Z Modeling Not Standardized}

Z.fit.NS <- lm(Z.Form, data = Modeling.Data.Filtered.NS, na.action = "na.fail")

Z.Models.NS <- MuMIn::dredge(Z.fit.NS, beta = "none", m.lim = c(1,NA))

plot(MuMIn::get.models(Z.Models.NS,1)[[1]])

hist(residuals(MuMIn::get.models(Z.Models.NS,1)[[1]]))

summary(MuMIn::get.models(Z.Models.NS,1)[[1]])

confint(MuMIn::get.models(Z.Models.NS,1)[[1]])

summary(MuMIn::get.models(Z.Models.NS,2)[[1]])

confint(MuMIn::get.models(Z.Models.NS,2)[[1]])

```

```{r, Z Modeling Standardized}

Z.fit.S <- lm(Z.Form, data = Modeling.Data.Filtered.S, na.action = "na.fail")

Z.Models.S <- MuMIn::dredge(Z.fit.S, beta = "none", m.lim = c(1,NA))

plot(MuMIn::get.models(Z.Models.S,1)[[1]])

hist(residuals(MuMIn::get.models(Z.Models.S,1)[[1]]))

summary(MuMIn::get.models(Z.Models.S,1)[[1]])

confint(MuMIn::get.models(Z.Models.S,1)[[1]])

summary(MuMIn::get.models(Z.Models.S,2)[[1]])

confint(MuMIn::get.models(Z.Models.S,2)[[1]])

summary(MuMIn::get.models(Z.Models.S,3)[[1]])

confint(MuMIn::get.models(Z.Models.S,3)[[1]])

summary(MuMIn::get.models(Z.Models.S,4)[[1]])

confint(MuMIn::get.models(Z.Models.S,4)[[1]])
```


#### Stoddard and Taylor
```{r, G Formula creation}

G.Form <- names((Modeling.Data.Filtered.NS %>%
                   dplyr::select(7:ncol(.)))) %>% 
  paste(., collapse = " + ") %>%
  paste("G ~", ., sep = " ") %>%
  paste(., "ca:Time_Close.In", "ca:Time_Strength.In",
        "np:Time_Close.In", "np:Time_Close.In", sep = " + ")

```

```{r, G Modeling Not Standardized}

G.fit.NS <- lm(G.Form, data = Modeling.Data.Filtered.NS, na.action = "na.fail")

G.Models.NS <- MuMIn::dredge(G.fit.NS, beta = "none", m.lim = c(1,NA))

plot(MuMIn::get.models(G.Models.NS,1)[[1]])

hist(residuals(MuMIn::get.models(G.Models.NS,1)[[1]]))

summary(MuMIn::get.models(G.Models.NS,1)[[1]])

confint(MuMIn::get.models(G.Models.NS,1)[[1]])

summary(MuMIn::get.models(G.Models.NS,2)[[1]])

confint(MuMIn::get.models(G.Models.NS,2)[[1]])

```

```{r, G Modeling Standardized}

G.fit.S <- lm(G.Form, data = Modeling.Data.Filtered.S, na.action = "na.fail")

G.Models.S <- MuMIn::dredge(G.fit.S, beta = "none", m.lim = c(1,NA))

plot(MuMIn::get.models(G.Models.S,1)[[1]])

hist(residuals(MuMIn::get.models(G.Models.S,1)[[1]]))

summary(MuMIn::get.models(G.Models.S,1)[[1]])

confint(MuMIn::get.models(G.Models.S,1)[[1]])

summary(MuMIn::get.models(G.Models.S,2)[[1]])

confint(MuMIn::get.models(G.Models.S,2)[[1]])

summary(MuMIn::get.models(G.Models.S,3)[[1]])

confint(MuMIn::get.models(G.Models.S,3)[[1]])

summary(MuMIn::get.models(G.Models.S,4)[[1]])

confint(MuMIn::get.models(G.Models.S,4)[[1]])
```

#### Hexp
```{r, Hexp Formula creation}

Hexp.Form <- names((Modeling.Data.Filtered.NS %>%
                      dplyr::select(7:ncol(.)))) %>% 
  paste(., collapse = " + ") %>%
  paste("Hexp ~", ., sep = " ") %>%
  paste(., "ca:Time_Close.In", "ca:Time_Strength.In",
        "np:Time_Close.In", "np:Time_Close.In", sep = " + ")

```

```{r, Hexp Modeling Not Standardized}
Hexp.fit.NS <- lm(Hexp.Form, data = Modeling.Data.Filtered.NS, na.action = "na.fail")

Hexp.Models.NS <- MuMIn::dredge(Hexp.fit.NS, beta = "none", m.lim = c(1,NA))

plot(MuMIn::get.models(Hexp.Models.NS,1)[[1]])

hist(residuals(MuMIn::get.models(Hexp.Models.NS,1)[[1]]))

summary(MuMIn::get.models(Hexp.Models.NS,1)[[1]])

confint(MuMIn::get.models(Hexp.Models.NS,1)[[1]])

summary(MuMIn::get.models(Hexp.Models.NS,2)[[1]])

confint(MuMIn::get.models(Hexp.Models.NS,2)[[1]])
```

```{r, Hexp Modleing Standardized}
Hexp.fit.S <- lm(Hexp.Form, data = Modeling.Data.Filtered.S, na.action = "na.fail")

Hexp.Models.S <- MuMIn::dredge(Hexp.fit.S, beta = "none", m.lim = c(1,NA))

plot(MuMIn::get.models(Hexp.Models.S,1)[[1]])

hist(residuals(MuMIn::get.models(Hexp.Models.S,1)[[1]]))

summary(MuMIn::get.models(Hexp.Models.S,1)[[1]])

confint(MuMIn::get.models(Hexp.Models.S,1)[[1]])

summary(MuMIn::get.models(Hexp.Models.S,2)[[1]])

confint(MuMIn::get.models(Hexp.Models.S,2)[[1]])
```

#### Ar
```{r, Ar Formula creation}

Ar.Form <- names((Modeling.Data.Filtered.NS %>%
                      dplyr::select(7:ncol(.)))) %>% 
  paste(., collapse = " + ") %>%
  paste("Ar ~", ., sep = " ") %>%
  paste(., "ca:Time_Close.In", "ca:Time_Strength.In",
        "np:Time_Close.In", "np:Time_Close.In", sep = " + ")

```

```{r, Ar Modeling Not Standardized}
Ar.fit.NS <- lm(Ar.Form, data = Modeling.Data.Filtered.NS, na.action = "na.fail")

Ar.Models.NS <- MuMIn::dredge(Ar.fit.NS, beta = "none", m.lim = c(1,NA))

plot(MuMIn::get.models(Ar.Models.NS,1)[[1]])

hist(residuals(MuMIn::get.models(Ar.Models.NS,1)[[1]]))

summary(MuMIn::get.models(Ar.Models.NS,1)[[1]])

confint(MuMIn::get.models(Ar.Models.NS,1)[[1]])

summary(MuMIn::get.models(Ar.Models.NS,2)[[1]])

confint(MuMIn::get.models(Ar.Models.NS,2)[[1]])
```

```{r, Ar Modleing Standardized}
Ar.fit.S <- lm(Ar.Form, data = Modeling.Data.Filtered.S, na.action = "na.fail")

Ar.Models.S <- MuMIn::dredge(Ar.fit.S, beta = "none", m.lim = c(1,NA))

plot(MuMIn::get.models(Ar.Models.S,1)[[1]])

hist(residuals(MuMIn::get.models(Ar.Models.S,1)[[1]]))

summary(MuMIn::get.models(Ar.Models.S,1)[[1]])

confint(MuMIn::get.models(Ar.Models.S,1)[[1]])

summary(MuMIn::get.models(Ar.Models.S,2)[[1]])

confint(MuMIn::get.models(Ar.Models.S,2)[[1]])
```

#### Hobs
```{r, Hobs Formula creation}

Hobs.Form <- names((Modeling.Data.Filtered.NS %>%
                      dplyr::select(7:ncol(.)))) %>% 
  paste(., collapse = " + ") %>%
  paste("Hobs ~", ., sep = " ") %>%
  paste(., "ca:Time_Close.In", "ca:Time_Strength.In",
        "np:Time_Close.In", "np:Time_Close.In", sep = " + ")

```

```{r, Hobs Modeling Not Standardized}
Hobs.fit.NS <- lm(Hobs.Form, data = Modeling.Data.Filtered.NS, na.action = "na.fail")

Hobs.Models.NS <- MuMIn::dredge(Hobs.fit.NS, beta = "none", m.lim = c(1,NA))

plot(MuMIn::get.models(Hobs.Models.NS,1)[[1]])

hist(residuals(MuMIn::get.models(Hobs.Models.NS,1)[[1]]))

summary(MuMIn::get.models(Hobs.Models.NS,1)[[1]])

confint(MuMIn::get.models(Hobs.Models.NS,1)[[1]])

summary(MuMIn::get.models(Hobs.Models.NS,2)[[1]])

confint(MuMIn::get.models(Hobs.Models.NS,2)[[1]])
```

```{r, Hexp Modleing Standardized}
Hobs.fit.S <- lm(Hobs.Form, data = Modeling.Data.Filtered.S, na.action = "na.fail")

Hobs.Models.S <- MuMIn::dredge(Hobs.fit.S, beta = "none", m.lim = c(1,NA))

plot(MuMIn::get.models(Hobs.Models.S,1)[[1]])

hist(residuals(MuMIn::get.models(Hobs.Models.S,1)[[1]]))

summary(MuMIn::get.models(Hobs.Models.S,1)[[1]])

confint(MuMIn::get.models(Hobs.Models.S,1)[[1]])

summary(MuMIn::get.models(Hobs.Models.S,2)[[1]])

confint(MuMIn::get.models(Hobs.Models.S,2)[[1]])
```

# Figures

```{r, Map}

# the base world data for maps
world <- rnaturalearth::ne_countries(scale = "large", returnclass = "sf")

Map_df <- Population.Genetic.Data %>%
  dplyr::select(Population_ID, Longitude, Latitude) %>%
  distinct(.) %>%
  dplyr::mutate(Population_Name = dplyr::case_when(Population_ID == "NS" ~ "New Smyrna (NS)",
                                                   Population_ID == "EG" ~ "Everglades (EG)",
                                                   Population_ID == "LB" ~ "Lower Bogue (LB)",
                                                   Population_ID == "EI" ~ "Exuma Island (EI)",
                                                   Population_ID == "FL" ~ "Fort Lauderdale (FL)",
                                                   Population_ID == "IR" ~ "Indian River (IR)",
                                                   Population_ID == "UK" ~ "Upper Keys (UK)",
                                                   Population_ID == "SL" ~ "Saint Lucie (SL)",
                                                   Population_ID == "LC" ~ "Long Caye (LC)",
                                                   Population_ID == "TC" ~ "Twin Cayes (TC)",
                                                   Population_ID == "TA" ~ "Turneffe Atoll (TA)",
                                                   Population_ID == "SS" ~ "San Salvador (SS)",
                                                   Population_ID == "CC" ~ "Charlotte County (CC)",
                                                   Population_ID == "LK" ~ "Lower Keys (LK)",
                                                   Population_ID == "UH" ~ " Utila, Honduras (UH)",
                                                   Population_ID == "NC" ~ "Northern Caye (NC)",
                                                   Population_ID == "TB" ~ "Tampa Bay (TB)"))

# make map
Sample_Map <- ggplot2::ggplot(data = world) +
  # add world polygons
  geom_sf(fill= "antiquewhite", 
            color = "black") +
  coord_sf(xlim = c(-90,-70.5),
             ylim = c(15.6,31)) +
  # all of the next sections are just adding the circles
  ggforce::geom_mark_ellipse(data = Map_df, aes(x = Longitude, y = Latitude,
                                                                   group = Population_ID, label = Population_Name,
                                                                   x0 = -80,
                                                                   y0 = 31,
                                                                   filter = Population_ID == "NS"),
                             label.fontface = "plain",
                             con.cap = 0, expand = 0.012, 
                             label.fontsize = 8, label.family = "serif", con.border = "one", 
                             label.buffer = unit(0.000001, "mm"), 
                             label.margin = margin(0.5, 0.5, 0.5, 0.5, "mm")) +
  ggforce::geom_mark_ellipse(data = Map_df, aes(x = Longitude, y = Latitude,
                                                                   group = Population_ID, label = Population_Name,
                                                                   x0 = -75.5,
                                                                   y0 = 29.5,
                                                                   filter = Population_ID == "IR"), 
                             label.fontface = "plain",
                             con.cap = 0, expand = 0.012, 
                             label.fontsize = 8, label.family = "serif", con.border = "one", 
                             label.buffer = unit(0.000001, "mm"), 
                             label.margin = margin(0.5, 0.5, 0.5, 0.5, "mm")) +
  ggforce::geom_mark_ellipse(data = Map_df, aes(x = Longitude, y = Latitude,
                                                                   group = Population_ID, label = Population_Name,
                                                                   x0 = -76,
                                                                   y0 = 28.5,
                                                                   filter = Population_ID == "SL"),
                             label.fontface = "plain",
                             con.cap = 0, expand = 0.012, 
                             label.fontsize = 8, label.family = "serif", con.border = "one", 
                             label.buffer = unit(0.000001, "mm"), 
                             label.margin = margin(0.5, 0.5, 0.5, 0.5, "mm")) +
  ggforce::geom_mark_ellipse(data = Map_df, aes(x = Longitude, y = Latitude,
                                                                   group = Population_ID, label = Population_Name,
                                                                   x0 = -76,
                                                                   y0 = 28,
                                                                   filter = Population_ID == "FL"), 
                             label.fontface = "plain",
                             con.cap = 0, expand = 0.012, 
                             label.fontsize = 8, label.family = "serif", con.border = "one", 
                             label.buffer = unit(0.000001, "mm"), 
                             label.margin = margin(0.5, 0.5, 0.5, 0.5, "mm")) +
  ggforce::geom_mark_ellipse(data = Map_df, aes(x = Longitude, y = Latitude,
                                                                   group = Population_ID, label = Population_Name,
                                                                   x0 = -73,
                                                                   y0 = 27,
                                                                   filter = Population_ID == "LB"),
                             label.fontface = "plain",
                             con.cap = 0, expand = 0.012, 
                             label.fontsize = 8, label.family = "serif", con.border = "one", 
                             label.buffer = unit(0.000001, "mm"), 
                             label.margin = margin(0.5, 0.5, 0.5, 0.5, "mm")) +
  ggforce::geom_mark_ellipse(data = Map_df, aes(x = Longitude, y = Latitude,
                                                                   group = Population_ID, label = Population_Name,
                                                                   x0 = -78.5,
                                                                   y0 = 25,
                                                                   filter = Population_ID == "EI"),
                             label.fontface = "plain",
                             con.cap = 0, expand = 0.012, 
                             label.fontsize = 8, label.family = "serif", con.border = "one", 
                             label.buffer = unit(0.000001, "mm"), 
                             label.margin = margin(0.5, 0.5, 0.5, 0.5, "mm")) +
  ggforce::geom_mark_ellipse(data = Map_df, aes(x = Longitude, y = Latitude,
                                                                   group = Population_ID, label = Population_Name,
                                                                   x0 = -73,
                                                                   y0 = 25,
                                                                   filter = Population_ID == "SS"),
                             label.fontface = "plain",
                             con.cap = 0, expand = 0.012, 
                             label.fontsize = 8, label.family = "serif", con.border = "one", 
                             label.buffer = unit(0.000001, "mm"), 
                             label.margin = margin(0.5, 0.5, 0.5, 0.5, "mm")) +
  ggforce::geom_mark_ellipse(data = Map_df, aes(x = Longitude, y = Latitude,
                                                                   group = Population_ID, label = Population_Name,
                                                                   x0 = -84.5,
                                                                   y0 = 30.8,
                                                                   filter = Population_ID == "TB"),
                             label.fontface = "plain",
                             con.cap = 0, expand = 0.012, 
                             label.fontsize = 8, label.family = "serif", con.border = "one", 
                             label.buffer = unit(0.000001, "mm"), 
                             label.margin = margin(0.5, 0.5, 0.5, 0.5, "mm")) +
  ggforce::geom_mark_ellipse(data = Map_df, aes(x = Longitude, y = Latitude,
                                                                   group = Population_ID, label = Population_Name,
                                                                   x0 = -87,
                                                                   y0 = 29.5,
                                                                   filter = Population_ID == "CC"),
                             label.fontface = "plain",
                             con.cap = 0, expand = 0.012, 
                             label.fontsize = 8, label.family = "serif", con.border = "one", 
                             label.buffer = unit(0.000001, "mm"), 
                             label.margin = margin(0.5, 0.5, 0.5, 0.5, "mm")) +
  scale_colour_identity() +
  scale_fill_brewer(palette = "Paired") +
  # add scale bar and place in bottom right
  ggspatial::annotation_scale(location = "bl",
                   width_hint = 0.2,
                   text_family = "serif",
                   height = unit(1.5, "mm")) +
  # place north arrow in top right
  ggspatial::annotation_north_arrow(location = "tr", 
                         which_north = "true", 
        style = ggspatial::north_arrow_fancy_orienteering) +
  #make dashed grid lines and fill the ocean with blue
  theme(panel.grid.major = element_line(color = gray(.5), linetype = "dashed", size = 0.5), 
        panel.background = element_rect(fill = "aliceblue"),
        panel.border = element_rect(colour = "black", linewidth = 1, fill = NA),
        text= element_text(family="serif", size = 8),
        axis.title = element_text(size = 12),
        legend.position="none") +
  # add location points
  geom_point(data = Map_df, aes(x = Longitude, y = Latitude),
             shape = 21, fill = "darksalmon")  +
  # Add red box around Florida keys
  geom_rect(aes(xmin= -81.8,
                xmax= -80.2, 
                ymin= 24.4, 
                ymax= 25.5), 
            color="black", 
            linewidth = 0.2, 
            fill = NA) +
  geom_segment(aes(x = -81.8, y = 25.5, xend = -84.15, yend = 28.3), 
            color="black", 
            linewidth = 0.2) +
  geom_segment(aes(x = -81.8, y = 24.4, xend = -84.2, yend = 23.8), 
            color="black", 
            linewidth = 0.2) + 
  # Add red box around Central America
  geom_rect(aes(xmin= -89.7,
                xmax= -86.7, 
                ymin= 15.95, 
                ymax= 18), 
            color="black", 
            linewidth = 0.2, 
            fill = NA) +
  geom_segment(aes(x = -89.7, y = 15.95, xend = -90.9, yend = 18.1), 
            color="black", 
            linewidth = 0.2) +
  geom_segment(aes(x = -86.7, y = 15.95, xend = -83, yend = 18.1), 
            color="black", 
            linewidth = 0.2)


Map_SF_Insert <- ggplot2::ggplot(data = world) +
  # add world polygons
  geom_sf(fill= "antiquewhite", 
            color = "black") +
  coord_sf(xlim = c(-81.57, 
                    -80.45),
           ylim = c(24.6,
                    25.35)) +
  # adding the circles just repeated code
 ggforce::geom_mark_ellipse(data = Map_df, aes(x = Longitude, y = Latitude,
                                                                   group = Population_ID, label = Population_Name,
                                                                   x0 = -81.3,
                                                                   y0 = 24.8,
                                                                   filter = Population_ID == "LK"), 
                            label.fontface = "plain",
                             con.cap = 0, expand = 0.015, 
                             label.fontsize = 8, label.family = "serif", con.border = "one", 
                             label.buffer = unit(0.000001, "mm"), 
                             label.margin = margin(0.5, 0.5, 0.5, 0.5, "mm")) +
  ggforce::geom_mark_ellipse(data = Map_df, aes(x = Longitude, y = Latitude,
                                                                   group = Population_ID, label = Population_Name,
                                                                   x0 = -80.9,
                                                                   y0 = 25,
                                                                   filter = Population_ID == "UK"), 
                             label.fontface = "plain",
                             con.cap = 0, expand = 0.015, 
                             label.fontsize = 8, label.family = "serif", con.border = "one", 
                             label.buffer = unit(0.000001, "mm"), 
                             label.margin = margin(0.5, 0.5, 0.5, 0.5, "mm")) +
  ggforce::geom_mark_ellipse(data = Map_df, aes(x = Longitude, y = Latitude,
                                                                   group = Population_ID, label = Population_Name,
                                                                   x0 = -81.2,
                                                                   y0 = 25.3,
                                                                   filter = Population_ID == "EG"), 
                             label.fontface = "plain",
                             con.cap = 0, expand = 0.015, 
                             label.fontsize = 8, label.family = "serif", con.border = "one", 
                             label.buffer = unit(0.000001, "mm"), 
                             label.margin = margin(0.5, 0.5, 0.5, 0.5, "mm")) +
  # add location points
  geom_point(data = Map_df, aes(x = Longitude, y = Latitude),
             shape = 21, fill = "darksalmon")  + 
  scale_fill_brewer(palette = "Paired") +
  # add scale bar and place in bottom right
  ggspatial::annotation_scale(location = "br",
                   width_hint = 0.4,
                   height = unit(0.15, "cm"),
                   text_family = "serif") +
  # place north arrow in top right
  ggspatial::annotation_north_arrow(location = "tr", 
                         which_north = "true", 
        style = ggspatial::north_arrow_fancy_orienteering,
        height = unit(.7, "cm"),
        width = unit(.7, "cm"),) +
  theme_void() +
  #make dashed grid lines and fill the ocean with blue
  theme(panel.background = element_rect(fill = "aliceblue"),
        legend.position="none",
        #make border
        panel.border = element_rect(colour = "black",
                                    fill=NA),
        #plot margins
        plot.margin = margin(-5,-10,0,-10),
        #remove axis text
        axis.text = element_blank(),
        #remove title
        axis.title = element_blank(),
        #remove ticks
        axis.ticks=element_blank(),
        #change plot title to serif
        plot.title = element_text(hjust = 0.5,
                                  family = "serif",
                                  size = 10,
                                  face = "bold"),
        plot.background = element_rect(fill = "aliceblue",
                                       color = "transparent")) + 
  ggtitle(label = "South Florida") 


# make central american insert
Map_CA_Insert <- ggplot2::ggplot(data = world) +
  # add world polygons
  geom_sf(fill= "antiquewhite", 
            color = "black") +
  coord_sf(xlim = c(-89.7, 
                    -86.7),
           ylim = c(16,
                    18)) + 
  # add the circles 
 ggforce::geom_mark_ellipse(data = Map_df, aes(x = Longitude, y = Latitude,
                                                                   group = Population_ID, label = Population_Name,
                                                                   x0 = -89,
                                                                   y0 = 16.3,
                                                                   filter = Population_ID == "UH"), 
                            label.fontface = "plain",
                             con.cap = 0, expand = 0.03, 
                             label.fontsize = 8, label.family = "serif", con.border = "one", 
                             label.buffer = unit(0.00000001, "mm"), 
                             label.margin = margin(0.5, 0.5, 0.5, 0.5, "mm")) +
  ggforce::geom_mark_ellipse(data = Map_df, aes(x = Longitude, y = Latitude,
                                                                   group = Population_ID, label = Population_Name,
                                                                   x0 = -89.6,
                                                                   y0 = 16.8,
                                                                   filter = Population_ID == "TC"),
                             label.fontface = "plain",
                             con.cap = 0, expand = 0.03, 
                             label.fontsize = 8, label.family = "serif", con.border = "one", 
                             label.buffer = unit(0.00000001, "mm"), 
                             label.margin = margin(0.5, 0.5, 0.5, 0.5, "mm")) +
  ggforce::geom_mark_ellipse(data = Map_df, aes(x = Longitude, y = Latitude,
                                                                   group = Population_ID, label = Population_Name,
                                                                   x0 = -89.5,
                                                                   y0 = 17.2,
                                                                   filter = Population_ID == "TA"), 
                             label.fontface = "plain",
                             con.cap = 0, expand = 0.03, 
                             label.fontsize = 8, label.family = "serif", con.border = "one", 
                             label.buffer = unit(0.0000001, "mm"), 
                             label.margin = margin(0.5, 0.5, 0.5, 0.5, "mm")) +
  ggforce::geom_mark_ellipse(data = Map_df, aes(x = Longitude, y = Latitude,
                                                                   group = Population_ID, label = Population_Name,
                                                                   x0 = -89.6,
                                                                   y0 = 17.6,
                                                                   filter = Population_ID == "NC"),
                             label.fontface = "plain",
                             con.cap = 0, expand = 0.03, 
                             label.fontsize = 8, label.family = "serif", con.border = "one", 
                             label.buffer = unit(0.000001, "mm"), 
                             label.margin = margin(0.5, 0.5, 0.5, 0.5, "mm")) +
  ggforce::geom_mark_ellipse(data = Map_df, aes(x = Longitude, y = Latitude,
                                                                   group = Population_ID, label = Population_Name,
                                                                   x0 = -86.7,
                                                                   y0 = 16.8,
                                                                   filter = Population_ID == "LC"),
                             label.fontface = "plain",
                             con.cap = 0, expand = 0.03, 
                             label.fontsize = 8, label.family = "serif", con.border = "one", 
                             label.buffer = unit(0.00000001, "mm"), 
                             label.margin = margin(0.5, 0.5, 0.5, 0.5, "mm")) +
  # add location points
  geom_point(data = Map_df, aes(x = Longitude, y = Latitude),
             shape = 21, fill = "darksalmon")  +
  # add scale bar and place in bottom right
  ggspatial::annotation_scale(location = "bl",
                   width_hint = 0.3,
                   height = unit(0.12, "cm"),
                   text_family = "serif") +
  # place north arrow in top right
  ggspatial::annotation_north_arrow(location = "tr", 
                         which_north = "true", 
        style = ggspatial::north_arrow_fancy_orienteering,
        height = unit(.7, "cm"),
        width = unit(.7, "cm"),) +
  theme_void() +
  scale_x_continuous(breaks= -88, minor_breaks=NULL) +
  scale_y_continuous(breaks= 18, minor_breaks=NULL) +
  scale_fill_brewer(palette = "Paired") +
  #make dashed grid lines and fill the ocean with blue
  theme(panel.background = element_rect(fill = "aliceblue"),
        legend.position="none",
        #make border
        panel.border = element_rect(colour = "black",
                                    fill=NA),
        #plot margins
        plot.margin = margin(-5,-10,0,-10),
        #remove axis text
        axis.text = element_blank(),
        #remove title
        axis.title = element_blank(),
        #remove ticks
        axis.ticks=element_blank(),
        #change plot title to serif
        plot.title = element_text(hjust = 0.5,
                                  family = "serif",
                                  size = 10,
                                  face = "bold"),
        plot.background = element_rect(fill = "aliceblue",
                                       color = "transparent")) + 
  ggtitle(label = "Central America") 

# make one map with cowplot
Sample_Map_Final <- cowplot::plot_grid(cowplot::plot_grid((cowplot::ggdraw() +
             cowplot::draw_plot(Sample_Map) +
               cowplot::draw_plot(Map_SF_Insert, 
                         x = .122,
                         y = .542, 
                         width = .235, 
                         height = .295) +
               cowplot::draw_plot(Map_CA_Insert, 
                         x = .119,
                         y = .245, 
                         width = .29, 
                         height = .3)))) 


```