--- title: 'MigConnectivity package (older functionality)' author: "Jeffrey A. Hostetler, Michael T. Hallworth" date: '`r Sys.Date()`' output: rmarkdown::html_vignette: default rmarkdown::pdf_document: default vignette: | %\VignetteIndexEntry{MigConnectivity package (older functionality)} %\usepackage[utf8]{inputenc} %\VignetteEngine{knitr::rmarkdown} --- ```{r, echo = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", tidy = FALSE) ``` ## Introduction NOTE: this vignette is a bit out of date. Everything in it should still work (let us know if it doesn't), but it doesn't highlight the new functions from versions 0.3 and 0.4, such as `estStrength`, `estTransition`, and `diffMC`. The new [master vignette](https://smbc-nzp.github.io/MigConnectivity/MigConnectivity.html) links to all of those. Technological advancements have spurred rapid growth in the study of migratory connectivity, the spatial and temporal linkage of migratory populations throughout the annual cycle. However, the lack of a quantitative definition for the strength of migratory connectivity (MC) has limited our ability to draw inferences across species, studies, and data types about the seasonal co-occurrence of populations. MC is a standardized metric to quantify migratory connectivity between two phases of the annual cycle. It is independent of data type and accounts for the relative abundance of populations distributed across a seasonal range. Negative values of MC indicate that individuals close to each other in one season are further apart in the other season. If MC = 0, no relationship exists between distances in one season and another; if MC = 1, the relative distances between individuals in one season are the same in the other, although the scale can differ. Three data inputs are needed to calculate MC, and one is optional: 1. The probabilities of movement between regions (i.e., transition probabilities); 1. Distances between regions within the two distinct seasonal ranges (e.g., a matrix of distances between all regions within the breeding range and a matrix of distances between all regions within the non-breeding range); 1. Relative abundance among regions within the seasonal range from which the transition probabilities originate (e.g., relative abundance among breeding regions); 1. (Optional) Total sample size of animals used to estimate the transition probabilities. MC takes the distances between the regions within the seasonal range from which the transition probabilities originate (e.g. breeding) and the distances between the regions within the other seasonal range (e.g. non-breeding) and approximates the correlation between individuals at those distances. Earlier methods (Mantel correlation; rM) used the correlation between two distance matrices of individual animals captured during the breeding phase and recaptured during the non-breeding phase as a measure of the strength of migratory connectivity. Our method builds on this method with a distance-based correlation coefficient approximation, but is not specific to a data type and uses transition probabilities from discreet regions, as opposed to distance matrices between individuals. Further, we include relative abundance within one seasonal range to account for uneven sampling among regions. Cohen, E. B., J.A. Hostetler, M.T. Hallworth, C.S. Rushing, T.S. Sillett, P.P. Marra. 2018. Quantifying the strength of migratory connectivity. Methods in Ecology and Evolution. 9:3. 513-524. ### Overview of Functions 1. `calcMC` calculates the strength of migratory connectivity (MC) with the transitional probabilities between regions, distances between regions within the two seasonal ranges, relative abundance among regions within the seasonal range from which the transition probabilities originate, and (optionally) the total sample size of the transitional probability data. 1. `estMC` estimates the strength of migratory connectivity (MC) with estimates of uncertainty. MC estimated either with transition probability estimates or individual assignments to regions from tracking data. Uses re-sampling to measure MC uncertainty from RMark transition probability matrix estimates and/or JAGS relative abundance MCMC samples OR geolocator and/or GPS data. 1. `simMove` Simulates annual movement of animals between regions across seasons for years and months with individual variability and strength of migratory connectivity (MC). Incorporates movement among regions within a season and movement among regions between seasons. Between seasons, animals either return to the same region each year or return to a different region in the subsequent year (dispersal). Between year dispersal rates occur during the first year (natal) or any subsequent year (breeding or non-breeding). Simulation does not incorporate births or deaths. ## Getting Started ### Installing 'MigConnectivity' from GitHub In order to install the vignette along with the package use the following code with `build_vignettes = TRUE`. **Note** it takes quite a bit longer to download the package when `build_vignettes = TRUE`. ```{r eval = FALSE} install.packages('devtools') devtools::install_github("SMBC-NZP/MigConnectivity", build_vignettes = TRUE) library(MigConnectivity) ``` ```{r,warning=FALSE,message=FALSE,echo=FALSE} library(MigConnectivity) ``` ### Functions `calcMC` Calculate strength of migratory connectivity (MC) `estMC` Estimate the strength of migratory connectivity (MC) while incorporating sampling uncertainty `simMove` Simulates annual and seasonal movement of animals between regions ## Examples ### Utility functions for use in vignette `mlogitMat`: Calculates probability matrix based on exponential decline with distance ```{r} mlogitMat <- function(slope, dist) { preMat <- exp(-slope/mean(dist)*dist) diag(preMat) <- 0 outMat <- t(apply(preMat, 1, function(x) x/(1 + sum(x)))) diag(outMat) <- 1 - rowSums(outMat) return(outMat) } ``` `mlogitMC`: Crude optimization function for developing migratory connectivity pattern based on migratory connectivity strength. Uses same argument (origin.abund) for relative or absolute abundance. ```{r} mlogitMC <- function(slope, MC.in, origin.dist, target.dist, origin.abund, sample.size) { nBreeding <- nrow(origin.dist) nWintering <- nrow(target.dist) psi <- mlogitMat(slope, origin.dist) if (any(psi<0)) return(5*slope^2) MC <- calcMC(origin.dist, target.dist, originRelAbund = origin.abund, psi, sampleSize = sample.size) return((MC.in - MC)^2) } ``` `calcStrengthInd`: rM function for individuals ```{r} calcStrengthInd <- function(originDist, targetDist, locations, resamp=1000, verbose = 0) { nInd <- dim(locations)[1] originDist2 <- targetDist2 <- matrix(0, nInd, nInd) for (i in 1:(nInd-1)) { for (j in (i+1):nInd) { originDist2[i,j] <- originDist[locations[i,1,1,1], locations[j,1,1,1]] targetDist2[i,j] <- targetDist[locations[i,2,1,1], locations[j,2,1,1]] originDist2[j,i] <- originDist[locations[i,1,1,1], locations[j,1,1,1]] targetDist2[j,i] <- targetDist[locations[i,2,1,1], locations[j,2,1,1]] } } return(ncf::mantel.test(originDist2, targetDist2, resamp=resamp, quiet = !verbose)) } ``` `calcPsiMC`: Simple approach to estimate psi matrix (transition probabilities) and MC from simulated or real data (does not incorporate sampling uncertainty) ```{r} calcPsiMC <- function(originDist, targetDist, originRelAbund, locations, years = 1, months = 1, verbose=FALSE) { nOrigin <- nrow(originDist) nTarget <- nrow(targetDist) psiMat <- matrix(0, nOrigin, nTarget) nInd <- dim(locations)[1] nYears <- dim(locations)[3] nMonths <- dim(locations)[4] for (i in 1:nInd) { if (i %% 1000 == 0 && verbose) # cat("Individual", i, "of", nInd, "\n") originMat <- locations[i, 1, years, months] targetMat <- locations[i, 2, years, months] bIndices <- which(!is.na(originMat)) wIndices <- which(!is.na(targetMat)) if (length(bIndices) && length(wIndices)) for (bi in bIndices) for (wi in wIndices) psiMat[originMat[bi], targetMat[wi]] <- psiMat[originMat[bi], targetMat[wi]] + 1 } psiMat <- apply(psiMat, 2, "/", rowSums(psiMat)) MC <- calcMC(originDist, targetDist, psi = psiMat, originRelAbund = originRelAbund, sampleSize = nInd) return(list(psi=psiMat, MC=MC)) } ``` `changeLocations`: transfers the simulated bird locations from the true populations to the researcher defined regions ```{r} changeLocations <- function(animalLoc, breedingSiteTrans, winteringSiteTrans) { animalLoc[,1,,] <- breedingSiteTrans[animalLoc[,1,,]] animalLoc[,2,,] <- winteringSiteTrans[animalLoc[,2,,]] return(animalLoc) } ``` `simLocationError`: Simulates location error with defined bias and variance / covariance matrix ```{r} simLocationError <- function(targetPoints, targetSites, geoBias, geoVCov, projection, verbose = 0, nSim = 1000) { if(!inherits(targetPoints,"sf")){ targetPoints <- sf::st_as_sf(targetPoints)} nAnimals <- nrow(targetPoints) geoBias2 <- matrix(rep(geoBias, nSim), nrow=nSim, byrow=TRUE) target.sample <- rep(NA, nAnimals) target.point.sample <- matrix(NA, nAnimals, 2) for(i in 1:nAnimals){ if (verbose > 0) cat('Animal', i, 'of', nAnimals) draws <- 0 while (is.na(target.sample[i])) { draws <- draws + 1 # Sample random point for each bird # from parametric distribution of NB error ps <- MASS::mvrnorm(n=nSim, mu=cbind(sf::st_coordinates(targetPoints)[i,1], sf::st_coordinates(targetPoints)[i,2]), Sigma=geoVCov)+ geoBias2 colnames(ps) <- c("Longitude","Latitude") point.sample <- sf::st_as_sf( data.frame(ps), coords = c("Longitude","Latitude"), crs = sf::st_crs(targetPoints)) # filtered to stay in NB areas (land) target.sample0 <- unlist(unclass(sf::st_intersects(point.sample, targetSites))) target.sample[i]<-target.sample0[!is.na(target.sample0)][1] } target.point.sample[i, ]<-sf::st_coordinates(point.sample[!is.na(target.sample0),])[1,] if (verbose > 0) cat(' ', draws, 'draws\n') } return(list(targetSample = target.sample, targetPointSample = target.point.sample)) } ``` `toPoly`: Helper function to convert a string of XY coordinates of centroids to polygons ```{r} toPoly <- function(siteCentroids, projection.in = 4326, projection.out = NA, resolution = NA, order.by.input = TRUE) { # This automatically sets the resolution so that all polygons touch and # cover the entire surface. Alternatively the user can supply the # resolution of the cells in the units of the input projection (projection.in) colnames(siteCentroids) <- c("Longitude","Latitude") if(is.na(resolution)){ long <- unique(siteCentroids[,1]) lat <- unique(siteCentroids[,2]) if (length(long)>1) long.res <- abs(long[2]-long[1]) else long.res <- 1 if (length(lat)>1) lat.res <- abs(lat[2]-lat[1]) else lat.res <- 1 resolution <- c(long.res,lat.res) } centers <- sf::st_as_sf(data.frame(siteCentroids), coords = c("Longitude","Latitude"), crs = projection.in) polys <- sf::st_as_sf(sf::st_make_grid(centers, crs = projection.in, cellsize = resolution, offset = apply(siteCentroids, 2, min) - resolution/2)) if(order.by.input){ # Reorders the polygons to match that of the input # orders <- suppressMessages(unlist(unclass(sf::st_intersects(centers, polys, sparse = TRUE)))) polys <- polys[orders,] } if(!is.na(projection.out)){ polys <- sf::st_transform(polys, projection.out) } return(polys) } ``` ### EXAMPLE 1 ### `calcMC` - Calculate strength of migratory connectivity This function calculates the strength of migratory connectivity between populations during two different time periods within the annual cycle. Below, we illustrate how to calculate the strength of MC between three breeding and non-breeding regions. To calculate the strength of MC, you will need: 1. to define the number of breeding and non-breeding regions 1. the distance between the breeding regions and the distance between the non-breeding regions 1. the transition probabilities between each breeding and non-breeding region (`psi`) 1. the relative abundance within each of the regions where individuals originate from (here the breeding regions) 1. (optional) the total sample size of animals used to estimate the transition probabilities. ### Calculate MC between breeding and non-breeding regions Simple example with three breeding and three non-breeding regions Define the number of breeding and non-breeding populations ```{r} nBreeding <- 3 #number of breeding regions nNonBreeding <- 3 #number of non-breeding regions ``` Generate a distance matrix between the different regions ```{r} #distances between centroids of three breeding regions breedDist <- matrix(c(0, 1, 2, 1, 0, 1, 2, 1, 0), nBreeding, nBreeding) #distances between centroids of three non-breeding regions nonBreedDist <- matrix(c(0, 5, 10, 5, 0, 5, 10, 5, 0), nNonBreeding, nNonBreeding) ``` ```{r echo=FALSE} op <- graphics::par(no.readonly = TRUE) on.exit(graphics::par(op)) par(bty="n") plot(c(4,5,6),rep(1,3), ylim=c(-0.05,1.05), xlim=c(-1,11), pch=19, axes=FALSE, yaxt="n", xaxt="n", cex=2, ylab="", xlab="Non-breeding", main="Breeding") points(c(0,5,10),rep(0,3),pch=15,cex=2) ``` Define the transition probabilities between the breeding and non-breeding regions ```{r} #transition probabilities form each breeding to each non-breeding region psi <- matrix(c(0.4, 0.35, 0.25, 0.3, 0.4, 0.3, 0.2, 0.3, 0.5), nBreeding, nNonBreeding, byrow = TRUE) ``` ```{r echo=FALSE} plot(c(4,5,6),rep(1,3), ylim=c(-0.05,1.05), xlim=c(-1,11), pch=19, axes=FALSE, yaxt="n", xaxt="n", cex=2, ylab="", xlab="Non-breeding", main="Breeding") points(c(0,5,10),rep(0,3),pch=15,cex=2) segments(4,1,0,0,lwd=0.4*10) segments(4,1,5,0,lwd=0.35*10) segments(4,1,10,0,lwd=0.25*10) segments(5,1,0,0,lwd=0.3*10, lty = 2) segments(5,1,5,0,lwd=0.4*10, lty = 2) segments(5,1,10,0,lwd=0.3*10, lty = 2) segments(6,1,0,0,lwd=0.2*10, lty = 3) segments(6,1,5,0,lwd=0.3*10, lty = 3) segments(6,1,10,0,lwd=0.5*10, lty = 3) ``` Define the relative abundance within the three breeding regions ```{r} #equal relative abundance among the three breeding regions, must sum to 1 relN <- rep(1/nBreeding, nBreeding) relN # Define total sample size for psi data # for small sample corrected version of MC n <- 250 ``` Calculate the strength of migratory connectivity using the inputs above ```{r} # Calculate the strength of migratory connectivity MC <- calcMC(originDist = breedDist, targetDist = nonBreedDist, originRelAbund = relN, psi = psi) round(MC, 3) MC <- calcMC(originDist = breedDist, targetDist = nonBreedDist, originRelAbund = relN, psi = psi, sampleSize = n) round(MC, 3) ``` ### EXAMPLE 2 ### `estMC` - Estimate strength of migratory connectivity incorporating location and other sampling uncertainty Estimates of MC may also be influenced by error in individual assignment to regions. Data types vary in location accuracy and precision, which are likely to influence the accuracy with which individuals are assigned to regions. To estimate MC and include location uncertainty the following data are needed: 1. A logical vector indicating the type of device the location estimates were derived from (GPS = FALSE, Light-level geolocator = TRUE) 1. Location Bias - a vector that has error estimates for both longitude and latitude 1. Location error - a variance, covariance matrix of longitude and latitude 1. Distance matrix between breeding and non-breeding regions 1. A shapefile of both breeding and non-breeding regions 1. The deployment locations and the 'unknown' locations derived from the tracking devices 1. The relative (or absolute) abundance within each region where the birds originate from (deployment regions) Below is a repeatable example for how to calculate location bias and location error using coordinates derived from light-level geolocators. ```{r, eval = FALSE} # Load in projections data("projections") # Define deployment locations (winter) # captureLocations<-matrix(c(-77.93,18.04, # Jamaica -80.94,25.13, # Florida -66.86,17.97), # Puerto Rico nrow=3, ncol=2, byrow = TRUE) colnames(captureLocations) <- c("Longitude","Latitude") # Convert capture locations into points # CapLocs <- sf::st_as_sf(data.frame(captureLocations), coords = c("Longitude","Latitude"), crs = 4326) # Project Capture locations # CapLocsM<-sf::st_transform(CapLocs, 'ESRI:54027') # Retrieve raw non-breeding locations from github # First grab the identity of the bird so we can loop through the files # For this example we are only interested in the error # around non-breeding locations # here we grab only the birds captured during the non-breeding season # Using paste0 for vignette formatting purposes winterBirds <- dget( "https://raw.githubusercontent.com/SMBC-NZP/MigConnectivity/master/data-raw/GL_NonBreedingFiles/winterBirds.txt" ) # create empty list to store the location data # Non_breeding_files <- vector('list',length(winterBirds)) # Get raw location data from Github # for(i in 1:length(winterBirds)){ Non_breeding_files[[i]] <- dget(paste0("https://raw.githubusercontent.com/", "SMBC-NZP/MigConnectivity/master/data-raw/", "GL_NonBreedingFiles/NonBreeding_", winterBirds[i],".txt")) } # Remove locations around spring Equinox and potential migration points # same NB time frame as Hallworth et al. 2015 # two steps because subset on shapefile doesn't like it in a single step Non_breeding_files <- lapply(Non_breeding_files, FUN = function(x){ month <- as.numeric(format(x$Date,format = "%m")) x[which(month != 3 & month != 4),]}) Jam <- c(1:9) # locations w/in list of winterBirds captured in Jamaica Fla <- c(10:12) # locations w/in list of winterBirds in Florida PR <- c(13:16) # locations w/in list of winterBirds in Puerto Rico # Turn the locations into shapefiles # NB_GL <- lapply(Non_breeding_files, FUN = function(x){ sf::st_as_sf(data.frame(x), coords = c("Longitude", "Latitude"), crs = 4326)}) # Project into UTM projection # NB_GLmeters <- lapply(NB_GL, FUN = function(x){sf::st_transform(x,'ESRI:54027')}) # Process to determine geolocator bias and variance-covariance in meters # # generate empty vector to store data # LongError<-rep(NA,length(winterBirds)) LatError<-rep(NA,length(winterBirds)) # Calculate the error in longitude derived # from geolocators from the true capture location LongError[Jam] <- unlist(lapply(NB_GLmeters[Jam], FUN = function(x){mean(sf::st_coordinates(x)[,1]- sf::st_coordinates(CapLocsM)[1,1])})) LongError[Fla] <- unlist(lapply(NB_GLmeters[Fla], FUN = function(x){mean(sf::st_coordinates(x)[,1]- sf::st_coordinates(CapLocsM)[2,1])})) LongError[PR] <- unlist(lapply(NB_GLmeters[PR], FUN = function(x){mean(sf::st_coordinates(x)[,1]- sf::st_coordinates(CapLocsM)[3,1])})) # Calculate the error in latitude derived from # geolocators from the true capture location LatError[Jam] <- unlist(lapply(NB_GLmeters[Jam], FUN = function(x){mean(sf::st_coordinates(x)[,2]- sf::st_coordinates(CapLocsM)[1,2])})) LatError[Fla] <- unlist(lapply(NB_GLmeters[Fla], FUN = function(x){mean(sf::st_coordinates(x)[,2]- sf::st_coordinates(CapLocsM)[2,2])})) LatError[PR] <- unlist(lapply(NB_GLmeters[PR], FUN = function(x){mean(sf::st_coordinates(x)[,2]- sf::st_coordinates(CapLocsM)[3,2])})) # Get co-variance matrix for error of # known non-breeding deployment sites # lm does multivariate normal models if you give it a matrix dependent variable! geo.error.model <- lm(cbind(LongError,LatError) ~ 1) geo.bias <- coef(geo.error.model) geo.vcov <- vcov(geo.error.model) ``` We measured MC for bootstrapped data of birds tracked from breeding to non-breeding regions using light-level and GPS geolocation when 1) GPS location uncertainty was applied to all individuals, 2) GPS and light-level location uncertainty were applied to individuals with those devices, and 3) light-level location uncertainty was applied to all individuals. Load location data that accompanies the `MigConnectivity` package. The location data are data from breeding Ovenbirds that were fit with Light-level geolocators or PinPoint-10 GPS tags. ```{r} data(OVENdata) # Ovenbird names(OVENdata) ``` The figure below shows the two breeding regions (squares), and the three non-breeding regions (gray scale) used in Cohen et al. (2018) to estimate MC for Ovenbirds tracked with light-level geolocators and PinPoint-10 GPS tags. ```{r, eval = FALSE} install.packages(c('shape', 'terra', 'maps')) ``` ```{r, message = FALSE, warning = FALSE, error=FALSE} library(shape) library(maps) library(terra) # Set the crs to WGS84 originSitesWGS84 <- sf::st_transform(OVENdata$originSites, 4326) targetSitesWGS84 <- sf::st_transform(OVENdata$targetSites, 4326) # Create a simple plot of the origin and and target sites par(mar=c(0,0,0,0)) plot(originSitesWGS84, xlim=c(terra::ext(targetSitesWGS84)[1], terra::ext(targetSitesWGS84)[2]), ylim=c(terra::ext(targetSitesWGS84)[3], terra::ext(originSitesWGS84)[4])) plot(targetSitesWGS84, add = TRUE, col=c("gray70","gray35","gray10")) map("world", add=TRUE) ``` The following code demonstrates how to estimate MC using location data from light-level geolocators and PinPoint-10 GPS tags. ```{r message=FALSE, warning = FALSE, error=FALSE} M<-estMC(isGL=OVENdata$isGL, # Logical vector: light-level geolocator(T)/GPS(F) geoBias = OVENdata$geo.bias, # Light-level geolocator location bias geoVCov = OVENdata$geo.vcov, #Light-level geolocator covariance matrix targetDist = OVENdata$targetDist, # Target location distance matrix originDist = OVENdata$originDist, # Origin location distance matrix targetSites = OVENdata$targetSites, # Non-breeding / target sites originSites = OVENdata$originSites, # Breeding / origin sites targetNames = OVENdata$targetNames, # Names of nonbreeding/target sites originNames = OVENdata$originNames, # Names of breeding/origin sites originPoints = OVENdata$originPoints, # Capture Locations targetPoints = OVENdata$targetPoints, # Target locations from devices originRelAbund = OVENdata$originRelAbund, # Origin relative abundances resampleProjection = sf::st_crs(OVENdata$targetPoints), verbose = 0, # output options - see help ??estMC nSamples = 10) # This is set low for example M ``` Take a closer look at what is included in the output ```{r} str(M, max.level = 2) ``` ```{r echo = FALSE, warning = FALSE, error = FALSE, eval = FALSE} plot(sf::st_geometry(wrld_simple), xlim=sf::st_bbox(OVENdata$targetSites)[c(1,3)], ylim=c(sf::st_bbox(OVENdata$targetSites)[2], sf::st_bbox(OVENdata$originSites)[4])) plot(sf::st_geometry(OVENdata$originSites[1,]),add=TRUE,lwd=1.75) plot(sf::st_geometry(OVENdata$originSites[2,]),add=TRUE,lwd=1.75) plot(sf::st_geometry(OVENdata$targetSites),add=TRUE,lwd=1.5,col=c("gray70","gray35","gray10")) legend(400000,-2105956,legend=paste("MC =",round(M$MC$mean,2), "\u00b1", round(M$MC$se,2)),bty="n",cex=1.8,bg="white",xjust=0) shape::Arrows(x0 = sf::st_coordinates(sf::st_centroid(OVENdata$originSites[1,]))[,1], y0 = sf::st_coordinates(sf::st_centroid(OVENdata$originSites[1,]))[,2], x1 = sf::st_coordinates(sf::st_centroid(OVENdata$targetSites[2,]))[,1]+80000, y1 = sf::st_bbox(OVENdata$targetSites[2,])[4]+150000, arr.length = 0.3, arr.adj = 0.5, arr.lwd = 1, arr.width = 0.4, arr.type = "triangle", lwd = M$psi$mean[1,2]*10, lty = 1) shape::Arrows(sf::st_coordinates(sf::st_centroid(OVENdata$originSites[1,]))[,1], sf::st_coordinates(sf::st_centroid(OVENdata$originSites[1,]))[,2], sf::st_coordinates(sf::st_centroid(OVENdata$targetSites[3,]))[,1], sf::st_bbox(OVENdata$targetSites[3,])[4]+150000, arr.length = 0.3, arr.adj = 0.5, arr.lwd = 1, arr.width = 0.4, arr.type = "triangle", lwd = M$psi$mean[1,3]*10, lty=1) shape::Arrows(sf::st_coordinates(sf::st_centroid(OVENdata$originSites[2,]))[,1], sf::st_coordinates(sf::st_centroid(OVENdata$originSites[2,]))[,2], sf::st_coordinates(sf::st_centroid(OVENdata$targetSites[1,]))[,1], sf::st_bbox(OVENdata$targetSites[1,])[4]+150000, arr.length = 0.3, arr.adj = 0.5, arr.lwd = 1, arr.width = 0.4, arr.type = "triangle", lwd = M$psi$mean[2,1]*10, lty=1) shape::Arrows(sf::st_coordinates(sf::st_centroid(OVENdata$originSites[2,]))[,1], sf::st_coordinates(sf::st_centroid(OVENdata$originSites[2,]))[,2], sf::st_coordinates(sf::st_centroid(OVENdata$targetSites[2,]))[,1]-80000, sf::st_bbox(OVENdata$targetSites[2,])[4]+150000, arr.length = 0.3, arr.adj = 0.5, arr.lwd = 1, arr.width = 0.4, arr.type = "triangle", lwd = M$psi$mean[2,2]*10) box(which="plot") ``` ### EXAMPLE 3 ### `simMove` Simulates position of animals by individual, season, year, and month Calculate the strength of migratory connectivity based on simulation data for 10 years of movement between breeding and non-breeding seasonal ranges. Breeding and non-breeding ranges are equally divided into 100 regions on an ellipsoid globe. See above for details regarding the utility functions. Simulate data to demonstrate the use of `calcMC` ```{r} set.seed(75) # Parameters for simulations nSeasons <- 2 # population with two discrete seasons - breeding / non-breeding nYears <- 10 # Ten years of data nMonths <- 4 # Four months within each season nBreeding <- 100 # Number of populations nWintering <- 100 # Number of populations ``` Generate the spatial arrangement of breeding and non-breeding populations ```{r} breedingPos <- matrix(c(rep(seq(-99,-81,2), each=sqrt(nBreeding)), rep(seq(49,31,-2), sqrt(nBreeding))), nBreeding, 2) winteringPos <- matrix(c(rep(seq(-79,-61,2), each=sqrt(nWintering)), rep(seq(9,-9,-2), sqrt(nWintering))), nWintering, 2) # calculate distance between populations breedDist <- distFromPos(breedingPos, 'ellipsoid') nonbreedDist <- distFromPos(winteringPos, 'ellipsoid') ``` The relative abundance of the study species is needed in each population. Here we generate those data below. ```{r} # Breeding Abundance breedingN <- rep(500, nBreeding) breedingRelN <- breedingN/sum(breedingN) ``` The transition probability (psi) between the breeding and non-breeding populations is calculated below. ```{r eval = FALSE} # Set up psi matrix o <- optimize(mlogitMC, MC.in = 0.25, origin.dist = breedDist, target.dist = nonbreedDist, origin.abund = breedingRelN, sample.size = sum(breedingN), interval = c(1, 3), tol = .Machine$double.eps^0.5)# slope <- o$minimum psi <- mlogitMat(slope, breedDist) ``` ```{r,warning=FALSE,message=FALSE,echo=FALSE} slope <- 1.92093332054914 psi <- mlogitMat(slope, breedDist) ``` Now that we have all the data needed to calculate the strength of migratory connectivity we can use the `calcMC` function in the `MigConnecitivty` package to generate a standardized MC metric. ```{r} # Baseline strength of migratory connectivity MC <- calcMC(originDist = breedDist, targetDist = nonbreedDist, psi = psi, originRelAbund = breedingRelN, sampleSize = sum(breedingN)) # Show Results MC ``` ### `simMove` Simulation of movement between breeding and non-breeding Recycle data generated for the calcMC function - see above ```{r} sims <- simMove(breedingAbund = breedingN, # Breeding relative abundance breedingDist = breedDist, # Breeding distance winteringDist = nonbreedDist, # Non-breeding distance psi = psi, # Transition probabilities nYears = nYears, # Number of years nMonths = nMonths, # Number of months winMoveRate = 0, # winter movement rate sumMoveRate = 0, # breeding movement rate winDispRate = 0, # winter dispersal rate sumDispRate = 0, # summer disperal rate natalDispRate = 0, # natal dispersal rate breedDispRate = 0, # breeding dispersal rate verbose = 0) # verbose output str(sims) ``` ## Additional Examples / Scenarios ### EXAMPLE 4 Calculate migratory connectivity for a range of input values Cohen et al. (*2018*) used simulation to assess the influence of the three data inputs needed to calculate MC, including: 1. matrix of transition probabilities from breeding regions to non-breeding regions 1. matrices of distances between breeding regions and between non-breeding regions 1. relative abundance among breeding regions. Input value simulations included transition probabilities from four breeding to four non-breeding regions and relative abundance measured among breeding regions. There were eight transition probability scenarios: 1. Full Mix 1. Avoid One Site 1. Full Connectivity 1. Half Mix 1. Low 1. Medium 1. One Site Preference 1. Negative There were twelve spatial arrangement / distance scenarios: 1. Base distances, linear/ linear 1. Distance between breeding sites 2 and 3 doubled 1. Distance between breeding sites 2 and 3 halved 1. Distance between breeding sites 3 and 4 doubled 1. Distance between breeding sites 3 and 4 halved 1. Breeding sites on square grid/ winter linear 1. Distance between wintering sites 2 and 3 doubled 1. Distance between wintering sites 2 and 3 halved 1. Distance between wintering sites 3 and 4 doubled 1. Distance between wintering sites 3 and 4 halved 1. Breeding linear, Wintering sites on square grid 1. Wintering and breeding on square grid There were five relative breeding abundance scenarios: 1. Base, all equal 1. Abundance at site B doubled 1. Abundance at site B halved 1. Abundance at site D doubled 1. Abundance at site D halved First we calculate MC for the eight transition probability scenarios (with base spatial arrangement and relative abundance): ```{r} nScenarios1 <- length(samplePsis) # samplePsis - comes with MigConnectivity package # Create vector of length nScenarios1 MC1 <- rep(NA, nScenarios1) # Loop through the different senarios outlined above # for (i in 1:nScenarios1) { MC1[i] <- calcMC(originDist = sampleOriginDist[[1]], targetDist = sampleTargetDist[[1]], psi = samplePsis[[i]], originRelAbund = sampleOriginRelN[[1]]) } # Give meaningful names to the MC1 vector names(MC1) <- names(samplePsis) # Print results round(MC1, 6) ``` Add the scenarios for the spatial arrangements that result in different distances between regions: ```{r} nScenarios2 <- length(sampleOriginPos) MC2 <- matrix(NA, nScenarios1, nScenarios2) rownames(MC2) <- names(samplePsis) colnames(MC2) <- names(sampleOriginPos) for (i in 1:nScenarios1) { for (j in 1:nScenarios2) { MC2[i, j] <- calcMC(originDist = sampleOriginDist[[j]], targetDist = sampleTargetDist[[j]], psi = samplePsis[[i]], originRelAbund = sampleOriginRelN[[1]]) } } # Print results # t(round(MC2, 4)) ``` Another way of comparing results: ```{r} MC.diff2 <- apply(MC2, 2, "-", MC2[ , 1]) t(round(MC.diff2, 4)) ``` ### Add the relative abundance breeding scenarios Calculate MC across abundance and transition probability scenarios for breeding linear / winter linear arrangement: ```{r} nScenarios3 <- length(sampleOriginRelN) MC3 <- matrix(NA, nScenarios1, nScenarios3) rownames(MC3) <- names(samplePsis) colnames(MC3) <- names(sampleOriginRelN) for (i in 1:nScenarios1) { for (j in 1) { for (k in 1:nScenarios3) { MC3[i, k] <- calcMC(originDist = sampleOriginDist[[j]], targetDist = sampleTargetDist[[j]], psi = samplePsis[[i]], originRelAbund = sampleOriginRelN[[k]]) } } } t(round(MC3, 4)) # linear arrangement ``` Calculate MC across abundance and transition probability scenarios for breeding grid/ winter grid arrangement: ```{r} MC4 <- matrix(NA, nScenarios1, nScenarios3) rownames(MC4) <- names(samplePsis) colnames(MC4) <- names(sampleOriginRelN) for (i in 1:nScenarios1) { for (j in nScenarios2) { for (k in 1:nScenarios3) { MC4[i, k] <- calcMC(originDist = sampleOriginDist[[j]], targetDist = sampleTargetDist[[j]], psi = samplePsis[[i]], originRelAbund = sampleOriginRelN[[k]]) } } } t(round(MC4, 4)) # grid arrangement ``` Calculate MC across abundance and transition probability scenarios for for breeding grid, winter linear arrangement ```{r} MC5 <- matrix(NA, nScenarios1, nScenarios3) rownames(MC5) <- names(samplePsis) colnames(MC5) <- names(sampleOriginRelN) for (i in 1:nScenarios1) { for (j in 6) { for (k in 1:nScenarios3) { MC5[i, k] <- calcMC(originDist = sampleOriginDist[[j]], targetDist = sampleTargetDist[[j]], psi = samplePsis[[i]], originRelAbund = sampleOriginRelN[[k]]) } } } t(round(MC5, 4)) # breeding grid, winter linear arrangement ``` ###EXAMPLES 5-12 (from sampling regime simulations in Cohen et al. 2018) Simulation was used to assess the influence of biased or incomplete sampling on the measurement of MC. Below, we measure the influence of several potential sources of sampling error: 1. Incorrect grouping of seasonal ranges into regions 1. Sampling not proportional to abundance 1. Migratory connectivity strength varies across a range 1. Low sample size 1. Uncertainty in assignments to regions from tracking data 1. Uncertainty in estimates of abundance 1. Heterogeneity in detection among regions from mark re-encounter data We also demonstrate how to propagate sampling uncertainty into measurement of MC using a combination of parametric and non-parametric bootstrapping. For most simulations, we compare MC to the results of a distance-based correlation coefficient (Mantel correlation; rM) that does not incorporate relative abundance but has been used as a measure of migratory connectivity. ###EXAMPLE 5 Sampling regime 1 of 3 Researchers divide populations differently than reality Delineation of seasonal ranges into regions I) Breeding range divided along equal longitudinal breaks into ten regions II) Non-breeding range divided along equal longitudinal breaks into ten regions III) Breeding and non-breeding ranges divided along equal longitudinal breaks into ten regions IV) Breeding range divided along the longitudinal and latitudinal midpoint into four regions V) Non-breeding range divided along the longitudinal and latitudinal midpoint into four regions VI) Breeding range divided along the longitudinal and latitudinal midpoint into four regions and non-breeding range divided along equal longitudinal breaks into ten regions VII) Breeding range divided along equal longitudinal breaks into ten regions and non-breeding range divided along the longitudinal and latitudinal midpoint into four regions ```{r, eval=FALSE} #Run functions and parameters above first set.seed(75) scenarios14 <- c("Base", "Breeding10", "Wintering10", "Breeding10Wintering10", "Breeding4", "Wintering4", "Breeding4Wintering10", "Breeding10Wintering4", "Breeding4Wintering4", "CentroidSampleBreeding10", "CentroidSampleBreeding10Wintering10", "CentroidSampleBreeding10Wintering4", "CentroidSampleBreeding4", "CentroidSampleBreeding4Wintering10", "CentroidSampleBreeding4Wintering4") # Each element is for a scenario (see above 1-8), transferring from natural # breeding populations to defined ones breedingSiteTrans14 <- list(1:nBreeding, rep(1:10, each=10), 1:nBreeding, rep(1:10, each=10), c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)), 1:nBreeding, c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)), rep(1:10, each=10), c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)), rep(1:10, each=10), rep(1:10, each=10), rep(1:10, each=10), c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)), c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)), c(rep(1:2, 5, each=5), rep(3:4, 5, each=5))) # Same for non-breeding populations winteringSiteTrans14 <- list(1:nWintering, 1:nWintering, rep(1:10, each=10), rep(1:10, each=10), 1:nWintering, c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)), rep(1:10, each=10), c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)), c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)), 1:nWintering, rep(1:10, each=10), c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)), 1:nWintering, rep(1:10, each=10), c(rep(1:2, 5, each=5), rep(3:4, 5, each=5))) # Examine the transfers in matrix form #lapply(breedingSiteTrans14, matrix, nrow=10, ncol=10) #lapply(winteringSiteTrans14, matrix, nrow=10, ncol=10) #positions of the human defined populations breedingPos14 <- list(breedingPos, rowsum(breedingPos, rep(1:10, each=10))/10, breedingPos, rowsum(breedingPos, rep(1:10, each=10))/10, rowsum(breedingPos, c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)))/25, breedingPos, rowsum(breedingPos, c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)))/25, rowsum(breedingPos, rep(1:10, each=10))/10, rowsum(breedingPos, c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)))/25, rowsum(breedingPos, rep(1:10, each=10))/10, rowsum(breedingPos, rep(1:10, each=10))/10, rowsum(breedingPos, rep(1:10, each=10))/10, rowsum(breedingPos, c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)))/25, rowsum(breedingPos, c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)))/25, rowsum(breedingPos, c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)))/25) winteringPos14 <- list(winteringPos, winteringPos, rowsum(winteringPos, rep(1:10, each=10))/10, rowsum(winteringPos, rep(1:10, each=10))/10, winteringPos, rowsum(winteringPos, c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)))/25, rowsum(winteringPos,rep(1:10, each=10))/10, rowsum(winteringPos, c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)))/25, rowsum(winteringPos, c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)))/25, winteringPos, rowsum(winteringPos, rep(1:10, each=10))/10, rowsum(winteringPos, c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)))/25, winteringPos, rowsum(winteringPos, rep(1:10, each=10))/10, rowsum(winteringPos, c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)))/25) # Calculate distances between defined breeding populations breedDist14 <- lapply(breedingPos14, distFromPos, surface = 'ellipsoid') # Calculate distances between defined non-breeding populations nonbreedDist14 <- lapply(winteringPos14, distFromPos, surface = 'ellipsoid') # Numbers of defined populations nBreeding14 <- c(100, 10, 100, 10, 4, 100, 4, 10, 4, 10, 10, 10, 4, 4, 4) nWintering14 <- c(100, 100, 10, 10, 100, 4, 10, 4, 4, 100, 10, 4, 100, 10, 4) # Relative abundance by scenario and breeding population breedingRelN14 <- lapply(nBreeding14, function(x) rep(1/x, x)) MC14 <- 0.25 # Same for all scenarios in this set nSample14 <- 1000 # Total number sampled per simulation # How many sampled from each natural population (sampling scenarios separate # from definition scenarios) sampleBreeding14 <- list(round(breedingRelN14[[1]]*nSample14), c(rep(0, 22), round(breedingRelN14[[5]][1]*nSample14), rep(0, 4), round(breedingRelN14[[5]][2]*nSample14), rep(0, 44), round(breedingRelN14[[5]][3]*nSample14), rep(0, 4), round(breedingRelN14[[5]][4]*nSample14), rep(0, 22)), rep(c(rep(0, 4), rep(round(breedingRelN14[[2]][1]*nSample14/2), 2), rep(0, 4)), 10)) # for the baseline use the simulation from above, sims animalLoc14base <- sims$animalLoc # Number of scenarios and number of simulations to run nScenarios14 <- length(breedingSiteTrans14) nSims14 <- 100 # Connections between scenarios and sampling regimes scenarioToSampleMap14 <- c(rep(1, 9), rep(3, 3), rep(2, 3)) # Set up data structures for storing results animalLoc14 <- vector("list", nScenarios14) #making an empty list to fill compare14 <- data.frame(Scenario = c("True", scenarios14), MC = c(MC14, rep(NA, nScenarios14)), Mantel = c(MC14, rep(NA, nScenarios14))) compare14.array <- array(NA, c(nSims14, nScenarios14, 2), dimnames = list(1:nSims14, scenarios14, c("MC", "Mantel"))) results14 <- vector("list", nScenarios14) # Run simulations for (sim in 1:nSims14) { cat("Simulation", sim, "of", nSims14, '\n') sim14 <- lapply(sampleBreeding14, simMove, breedingDist = breedDist14[[1]], winteringDist=nonbreedDist14[[1]], psi=psi, nYears=nYears, nMonths=nMonths) for (i in 1:nScenarios14) { cat("\tScenario", i, "\n") animalLoc14[[i]]<-changeLocations( sim14[[scenarioToSampleMap14[i]]]$animalLoc, breedingSiteTrans14[[i]], winteringSiteTrans14[[i]]) results14[[i]] <- calcPsiMC(breedDist14[[i]], nonbreedDist14[[i]], animalLoc14[[i]], originRelAbund = breedingRelN14[[i]], verbose = FALSE) compare14.array[sim, i, 'MC'] <- results14[[i]]$MC compare14.array[sim,i,'Mantel']<-calcStrengthInd(breedDist14[[1]], nonbreedDist14[[1]], sim14[[scenarioToSampleMap14[i]]]$animalLoc, resamp=0)$correlation } } # Compute means for each scenario compare14$MC[1:nScenarios14 + 1] <- apply(compare14.array[,,'MC'], 2, mean) compare14$Mantel[1:nScenarios14 + 1] <- apply(compare14.array[,,'Mantel'], 2, mean) compare14 <- transform(compare14, MC.diff=MC - MC[1], Mantel.diff=Mantel - Mantel[1]) compare14 ``` ###EXAMPLE 6 Sampling regime 2 of 3 Researchers divide populations differently than reality PLUS Different distributions of sampling animals across breeding range PLUS Sample sizes don't always match relative abundances PLUS Compare our approach and simple Mantel approach 1. Base (10 years, uneven abundances but matches sampling) 1. Breeding pops divided into 4 squares, sample across breeding range 1. Breeding pops divided into 4 squares, sample at centroid of each square 1. Sampling high in low abundance populations plus base 1. Scenarios 2 plus 4 1. Scenarios 3 plus 4 ```{r eval = FALSE} set.seed(75) # Transfer between true populations and researcher defined ones (only for # breeding, as not messing with winter populations here) breedingSiteTrans15 <- list(1:100, c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)), c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)), 1:100, c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)), c(rep(1:2, 5, each=5), rep(3:4, 5, each=5))) nScenarios15 <- length(breedingSiteTrans15) nSims15 <- 100 # Basing positions of researcher defined breeding populations on above breedingPos15 <- list(breedingPos, breedingPos14[[5]], breedingPos14[[5]], breedingPos, breedingPos14[[5]], breedingPos14[[5]]) winteringPos15 <- rep(list(winteringPos), nScenarios15) breedDist15 <- lapply(breedingPos15, distFromPos) nonbreedDist15 <- lapply(winteringPos15, distFromPos) nBreeding15 <- rep(c(100, 4, 4), 2) nWintering15 <- rep(100, nScenarios15) # Highest abundance in lower right corner, lowest in top left # Making symmetrical breedingN15base <- rep(NA, 100) for (i in 1:10) #row for (j in 1:10) #column breedingN15base[i+10*(j-1)] <- 500 + 850*i*j sum(breedingN15base) # For researcher defined populations breedingN15 <- lapply(breedingSiteTrans15, rowsum, x=breedingN15base) breedingRelN15 <- lapply(breedingN15, "/", sum(breedingN15base)) nSample15 <- 1000 # Total number sampled per simulation # Number sampled per natural population sampleBreeding15 <- list(round(breedingRelN15[[1]]*nSample15), c(rep(0, 22), round(breedingRelN15[[3]][1]*nSample15), rep(0, 4), round(breedingRelN15[[3]][2]*nSample15), rep(0, 44), round(breedingRelN15[[3]][3]*nSample15), rep(0, 4), round(breedingRelN15[[3]][4]*nSample15), rep(0, 22)), round(breedingRelN15[[1]]*nSample15)[100:1], c(rep(0, 22), round(breedingRelN15[[3]][1]*nSample15), rep(0, 4), round(breedingRelN15[[3]][2]*nSample15), rep(0, 44), round(breedingRelN15[[3]][3]*nSample15), rep(0, 4), round(breedingRelN15[[3]][4]*nSample15), rep(0, 22))[100:1]) # Set up psi matrix o15 <- optimize(mlogitMC, MC.in = 0.25, origin.dist = breedDist15[[1]], target.dist = nonbreedDist15[[1]], origin.abund = breedingN15[[1]]/sum(breedingN15[[1]]), sample.size = sum(breedingN15[[1]]), interval = c(0,10), tol = .Machine$double.eps^0.5) slope15 <- o15$minimum psi15 <- mlogitMat(slope15, breedDist15[[1]]) # Baseline strength of migratory connectivity MC15 <- calcMC(originDist = breedDist15[[1]], targetDist = nonbreedDist15[[1]], psi = psi15, originRelAbund = breedingN15[[1]]/sum(breedingN15[[1]]), sampleSize = sum(breedingN15[[1]])) # Run sampling regimes scenarioToSampleMap15 <- c(1, 1, 2, 3, 3, 4) animalLoc15 <- vector("list", nScenarios15) results15 <- vector("list", nScenarios15) compare15 <- data.frame(Scenario = c("True", "Base", "Breeding4", "CentroidSampleBreeding4", "BiasedSample", "BiasedSampleBreeding4", "BiasedCentroidSampleBreeding4"), MC = c(MC15, rep(NA, nScenarios15)), Mantel = c(MC15, rep(NA, nScenarios15))) compare15.array <- array(NA, c(nSims15, nScenarios15, 2), dimnames = list(1:nSims15, c("Base", "Breeding4", "CentroidSampleBreeding4", "BiasedSample", "BiasedSampleBreeding4", "BiasedCentroidSampleBreeding4"), c("MC", "Mantel"))) for (sim in 1:nSims15) { cat("Simulation", sim, "of", nSims15, '\n') sim15 <- lapply(sampleBreeding15, simMove, breedingDist = breedDist15[[1]], winteringDist=nonbreedDist15[[1]], psi=psi15, nYears=nYears, nMonths=nMonths) for (i in 1:nScenarios15) { cat("\tScenario", i, "\n") animalLoc15[[i]] <- changeLocations( animalLoc=sim15[[scenarioToSampleMap15[i]]]$animalLoc, breedingSiteTrans = breedingSiteTrans15[[i]], winteringSiteTrans = 1:nWintering15[i]) results15[[i]] <- calcPsiMC(originDist = breedDist15[[i]], targetDist = nonbreedDist15[[i]], originRelAbund = breedingRelN15[[i]], locations = animalLoc15[[i]], verbose = FALSE) compare15.array[sim, i, 'MC'] <- results15[[i]]$MC compare15.array[sim, i, 'Mantel'] <- calcStrengthInd(breedDist15[[1]], nonbreedDist15[[1]], sim15[[scenarioToSampleMap15[i]]]$animalLoc, resamp=0)$correlation } } compare15$MC[1:nScenarios15+1] <- apply(compare15.array[,,'MC'], 2, mean, na.rm = TRUE) compare15$Mantel[1:nScenarios15+1] <- apply(compare15.array[,,'Mantel'], 2, mean, na.rm = TRUE) compare15 <- transform(compare15, MC.diff=MC - MC[1], Mantel.diff=Mantel - Mantel[1], MC.prop=MC/MC[1], Mantel.prop=Mantel/Mantel[1]) compare15a <- as.matrix(compare15[1 + 1:nScenarios15, c("MC", "MC.diff", "Mantel", "Mantel.diff")]) rownames(compare15a) <- compare15$Scenario[1 + 1:nScenarios15] ``` ###EXAMPLE 7 Sampling regime 3 of 3 Researchers divide populations differently than reality (simulations) PLUS Different distributions of sampled animals across breeding range PLUS Sample sizes don't always match relative abundances PLUS Compare our approach and simple Mantel approach PLUS MC not same across subsections of range 1. Base (uneven MC (0.15 for NW breeding, 0.3 for SW, 0.45 for NE, and 0.6 for SE), uneven abundances (lowest in NW, highest in SE), sampling proportional to abundance 1. Breeding pops divided into 4 squares, sample across breeding range 1. Breeding pops divided into 4 squares, sample at centroid of each square 1. Sampling high in low abundance populations 1. Scenarios 2 plus 4 1. Scenarios 3 plus 4 ```{r eval = FALSE} set.seed(75) # Transfer between true populations and researcher defined ones # (only for breeding, as not messing with winter populations here) breedingSiteTrans16 <- list(1:100, c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)), c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)), 1:100, c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)), c(rep(1:2, 5, each=5), rep(3:4, 5, each=5))) #lapply(breedingSiteTrans16, matrix, nrow=10, ncol=10) nScenarios16 <- length(breedingSiteTrans16) nSims16 <- 100 # Basing positions of researcher defined breeding populations on above breedingPos16 <- breedingPos15 winteringPos16 <- winteringPos15 breedDist16 <- breedDist15 nonbreedDist16 <- nonbreedDist15 nBreeding16 <- nBreeding15 nWintering16 <- rep(100, nScenarios16) # Highest abundance in lower right corner, lowest in top left # In fact basing on distance from top left population breedingN16base <- breedingN15base breedingN16 <- breedingN15 breedingRelN16 <- lapply(breedingN16, "/", sum(breedingN16base)) # Set up psi matrix # Each quadrant of breeding range has different MC MC.levels16 <- seq(0.15, 0.6, 0.15) nLevels16 <- 4 psi16 <- matrix(NA, nBreeding16[1], nWintering16[1]) for (i in 1:nLevels16) { cat("MC", MC.levels16[i]) # Find a psi matrix that produces the given MC (for whole species) o16a <- optimize(mlogitMC, MC.in = MC.levels16[i], origin.dist = breedDist16[[1]], target.dist = nonbreedDist16[[1]], origin.abund = breedingN16[[1]]/sum(breedingN16[[1]]), sample.size = sum(breedingN16[[1]]), interval=c(0,10), tol=.Machine$double.eps^0.5) slope16a <- o16a$minimum cat(" slope", slope16a, "\n") psi16a <- mlogitMat(slope16a, breedDist16[[1]]) # Then use the rows of that psi matrix only for the one breeding quadrant rows <- 50*(i %/% 3) + rep(1:5, 5) + rep(seq(0, 40, 10), each=5) + ((i-1) %% 2) * 5 psi16[rows, ] <- psi16a[rows, ] } # Baseline strength of migratory connectivity MC16 <- calcMC(originDist = breedDist16[[1]], targetDist = nonbreedDist16[[1]], psi = psi16, originRelAbund = breedingN16[[1]]/sum(breedingN16[[1]]), sampleSize = sum(breedingN16[[1]])) # Set up sampling regimes (different number than number of scenarios) nSample16 <- 1000 sampleBreeding16 <- list(round(breedingRelN16[[1]]*nSample16), c(rep(0, 22), round(breedingRelN16[[3]][1]*nSample16), rep(0, 4), round(breedingRelN16[[3]][2]*nSample16), rep(0, 44), round(breedingRelN16[[3]][3]*nSample16), rep(0, 4), round(breedingRelN16[[3]][4]*nSample16), rep(0, 22)), round(breedingRelN16[[1]]*nSample16)[100:1], c(rep(0, 22), round(breedingRelN16[[3]][1]*nSample16), rep(0, 4), round(breedingRelN16[[3]][2]*nSample16), rep(0, 44), round(breedingRelN16[[3]][3]*nSample16), rep(0, 4), round(breedingRelN16[[3]][4]*nSample16), rep(0, 22))[100:1]) sapply(sampleBreeding16, sum) # Run sampling regimes scenarioToSampleMap16 <- c(1, 1, 2, 3, 3, 4) animalLoc16 <- vector("list", nScenarios16) results16 <- vector("list", nScenarios16) compare16 <- data.frame(Scenario = c("True", "Base", "Breeding4", "CentroidSampleBreeding4", "BiasedSample", "BiasedSampleBreeding4", "BiasedCentroidSampleBreeding4"), MC = c(MC16, rep(NA, nScenarios16)), Mantel = c(MC16, rep(NA, nScenarios16))) compare16.array <- array(NA, c(nSims16, nScenarios16, 2), dimnames = list(1:nSims16, c("Base", "Breeding4", "CentroidSampleBreeding4", "BiasedSample", "BiasedSampleBreeding4", "BiasedCentroidSampleBreeding4"), c("MC", "Mantel"))) set.seed(80) for (sim in 1:nSims16) { cat("Simulation", sim, "of", nSims16, '\n') sim16 <- lapply(sampleBreeding16, simMove, breedingDist = breedDist16[[1]], winteringDist=nonbreedDist16[[1]], psi=psi16, nYears=nYears, nMonths=nMonths) for (i in 1:nScenarios16) { cat("\tScenario", i, "\n") animalLoc16[[i]]<-changeLocations(sim16[[scenarioToSampleMap16[i]]]$animalLoc, breedingSiteTrans16[[i]], 1:nWintering16[i]) results16[[i]] <- calcPsiMC(originDist = breedDist16[[i]], targetDist = nonbreedDist16[[i]], locations = animalLoc16[[i]], originRelAbund = breedingRelN16[[i]], verbose = FALSE) compare16.array[sim, i, 'MC'] <- results16[[i]]$MC compare16.array[sim, i, 'Mantel'] <- calcStrengthInd(breedDist16[[1]], nonbreedDist16[[1]], sim16[[scenarioToSampleMap16[i]]]$animalLoc, resamp=0)$correlation } } compare16$MC[1:nScenarios16+1] <- apply(compare16.array[,,'MC'], 2, mean, na.rm = TRUE) compare16$Mantel[1:nScenarios16+1] <- apply(compare16.array[,,'Mantel'], 2, mean, na.rm = TRUE) compare16 <- transform(compare16, MC.diff=MC - MC[1], Mantel.diff=Mantel - Mantel[1]) compare16a <- as.matrix(compare16[1 + 1:nScenarios16, c('MC','MC.diff', "Mantel", "Mantel.diff")]) rownames(compare16a) <- compare16$Scenario[1 + 1:nScenarios16] ``` ###EXAMPLE 8 `estMC` Estimate strength of migratory connectivity incorporating uncertainty from limited sample size (no other sampling error). Here we simulated data using the 13th scenario from example 5 above ("CentroidSampleBreeding4"; make sure you run example 5 code before running this section). We changed the total sample size from 1000 to 100 animals and the total abundance from 50000 to 2500 animals. ```{r, eval = FALSE} # Load in projections data("projections") set.seed(75) capLocs14 <- lapply(breedingPos14, FUN = function(x){ colnames(x)<-c("Longitude","Latitude") sf::st_as_sf(data.frame(x), coords = c("Longitude","Latitude"), crs = 4326)}) targLocs14 <- lapply(winteringPos14, FUN = function(x){ colnames(x)<-c("Longitude","Latitude") sf::st_as_sf(data.frame(x), coords = c("Longitude","Latitude"), crs = 4326)}) # Relative abundance by scenario and breeding population breedingN14base <- rep(25, nBreeding14[1]) breedingN14 <- lapply(breedingSiteTrans14, rowsum, x=breedingN14base) breedingRelN14 <- lapply(breedingN14, "/", sum(breedingN14base)) MC14 <- 0.25 nSample14.1 <- 100 # Total number sampled per simulation # How many sampled from each natural population (sampling scenarios separate # from definition scenarios) sampleBreeding14.1 <- list(round(breedingRelN14[[1]]*nSample14.1), c(rep(0,22),round(breedingRelN14[[5]][1]*nSample14.1), rep(0,4), round(breedingRelN14[[5]][2]*nSample14.1), rep(0,44),round(breedingRelN14[[5]][3]*nSample14.1), rep(0,4), round(breedingRelN14[[5]][4]*nSample14.1), rep(0, 22)), rep(c(rep(0, 4), rep(round(breedingRelN14[[2]][1]*nSample14.1/2), 2), rep(0, 4)), 10)) lapply(sampleBreeding14.1, matrix, nrow=10, ncol=10) # Number of simulations to run nSims14 <- 100 nSimsLarge14 <- 250 #0riginally 2500 nYears <- 1 nMonths <- 1 # Set up data structures for storing results animalLoc14 <- vector("list", nScenarios14) #making an empty list to fill sim14 <- vector("list", nSimsLarge14) #making an empty list to fill compare14.1.array <- array(NA, c(nSimsLarge14, 2), dimnames = list(1:nSimsLarge14, c("MC", "Mantel"))) results14 <- vector("list", nScenarios14) # Run simulations set.seed(7) system.time(for (sim in 1:nSimsLarge14) { cat("Simulation", sim, "of", nSimsLarge14, '\n') sim14[[sim]] <- lapply(sampleBreeding14.1, simMove, breedingDist = breedDist14[[1]], winteringDist=nonbreedDist14[[1]], psi=psi, nYears=nYears, nMonths=nMonths) for (i in 13) { animalLoc14[[i]] <- changeLocations(sim14[[sim]][[scenarioToSampleMap14[i]]]$animalLoc, breedingSiteTrans14[[i]], winteringSiteTrans14[[i]]) results14[[i]] <- calcPsiMC(breedDist14[[i]], nonbreedDist14[[i]], animalLoc14[[i]], originRelAbund = breedingRelN14[[i]], verbose = FALSE) compare14.1.array[sim,'MC'] <- results14[[i]]$MC compare14.1.array[sim,'Mantel'] <- calcStrengthInd(breedDist14[[1]], nonbreedDist14[[1]], sim14[[sim]][[scenarioToSampleMap14[i]]]$animalLoc, resamp=0)$correlation } }) means14 <- apply(compare14.1.array, 2, mean) vars14 <- apply(compare14.1.array, 2, var) rmse14 <- apply(compare14.1.array, 2, function(x) sqrt(mean((x - MC14)^2))) # Set up data structures for storing estimation results est14.array <- array(NA, c(nSims14, 2), dimnames = list(1:nSims14,c("MC", "Mantel"))) var14.array <- array(NA, c(nSims14, 2), dimnames =list(1:nSims14,c("MC", "Mantel"))) ci14.array <- array(NA, c(nSims14, 2, 2), dimnames = list(1:nSims14, c("MC", "Mantel"), c('lower', 'upper'))) #making an empty list to fill animalLoc14 <- vector("list", nSims14) results14 <- vector("list", nSims14) # Run estimations set.seed(567) sim14.sub <- sim14[sample.int(nSimsLarge14, nSims14, TRUE)] for (sim in 1:nSims14) { cat("Estimation", sim, "of", nSims14, '\n') for (i in 13) {#:nScenarios14) { animalLoc14[[sim]] <- changeLocations(sim14.sub[[sim]][[scenarioToSampleMap14[i]]]$animalLoc, breedingSiteTrans14[[i]], winteringSiteTrans14[[i]]) results14[[sim]] <- estMC(originDist = breedDist14[[i]], targetDist = nonbreedDist14[[i]], originRelAbund = breedingRelN14[[i]], originPoints = capLocs14[[i]][animalLoc14[[sim]][, 1, 1, 1], ], targetPoints = targLocs14[[i]][animalLoc14[[sim]][, 2, 1, 1], ], originAssignment = animalLoc14[[sim]][, 1, 1, 1], targetAssignment = animalLoc14[[sim]][, 2, 1, 1], nSamples = 1000, verbose = 1, calcCorr = TRUE, geoBias = c(0, 0), geoVCov = matrix(0, 2, 2)) est14.array[sim, 'MC'] <- results14[[sim]]$MC$mean est14.array[sim, 'Mantel'] <- results14[[sim]]$corr$mean var14.array[sim, 'MC'] <- results14[[sim]]$MC$se ^ 2 var14.array[sim, 'Mantel'] <- results14[[sim]]$corr$se ^ 2 ci14.array[sim, 'MC', ] <- results14[[sim]]$MC$bcCI ci14.array[sim, 'Mantel', ] <- results14[[sim]]$corr$bcCI } } # Crude point estimates (bootstrap means are better) pointEsts14.1 <- t(sapply(results14, function(x) c(x$MC$point, x$corr$point))) # Summarize results summary(var14.array) vars14 summary(est14.array) summary(pointEsts14.1) means14 colMeans(pointEsts14.1) - MC14 sqrt(colMeans((pointEsts14.1 - MC14)^2)) summary(ci14.array[, "MC", 'lower'] <= MC14 & ci14.array[, "MC", 'upper'] >= MC14) summary(ci14.array[, "Mantel", 'lower'] <= MC14 & ci14.array[, "Mantel", 'upper'] >= MC14) # Plot histograms of bootstrap estimation results est.df <- data.frame(Parameter = rep(c("MC", 'rM'), 2, each = nSims14), Quantity = rep(c('Mean', 'Variance'), each = 2 * nSims14), sim = rep(1:nSims14, 4), Estimate = c(est14.array, var14.array)) trues.df <- data.frame(Parameter = rep(c("MC", 'rM'), 2), Quantity = rep(c('Mean', 'Variance'), each = 2), Value = c(MC14, MC14, vars14)) library(ggplot2) g.est <- ggplot(est.df, aes(Estimate)) + geom_histogram(bins = 15) + facet_grid(Parameter ~ Quantity, scales = 'free_x') + geom_vline(aes(xintercept = Value), data = trues.df) + theme_bw() + scale_x_continuous(breaks = c(0.002, 0.003, 0.1, 0.2, 0.3)) g.est # Summary table of bootstrap estimation results qualities14 <- data.frame(Parameter = rep(c("MC", 'rM'), 5), Quantity = rep(rep(c('Mean', 'Variance'), c(2, 2)), 3, 10), Measure = rep(c('Bias', 'RMSE', "Coverage"), c(4, 4, 2)), Value = c(colMeans(est14.array - MC14), colMeans(var14.array) - vars14, sqrt(colMeans((est14.array - MC14)^2)), sqrt(mean((var14.array[,1]-vars14[1])^2)), sqrt(mean((var14.array[,2]-vars14[2])^2)), mean(ci14.array[, "MC", 'lower'] <= MC14 & ci14.array[, "MC", 'upper'] >= MC14), mean(ci14.array[,"Mantel",'lower']<=MC14 & ci14.array[,"Mantel",'upper']>=MC14))) format(qualities14, digits = 2, scientific = FALSE) ``` ###EXAMPLE 9 `estMC` Estimate strength of migratory connectivity incorporating uncertainty from location error associated with light-level geolocators. Again, we simulated data using the 13th scenario from example 5 above ("CentroidSampleBreeding4"; make sure you run examples 5 and 8 code before running this section). ```{r, eval = FALSE} # Simulation using error associated with light-level geolocators # Assign geolocator bias / variance co-variance matrix geoBias <- c(-10000, 50000) geoVCov <- matrix(c(1.2e+8, 0, 0, 1.2e+9), 2, 2) sqrt(diag(geoVCov)) targLocs14.2 <- lapply(targLocs14, FUN = function(x){sf::st_transform(x, 6933)}) # convert wintering locations to polygons using the helper function WinteringPolys <- lapply(winteringPos14, toPoly, projection.in = 4326, projection.out = 6933) winteringPolys <- lapply(winteringPos14, toPoly, projection.in = 4326, projection.out = 4326) winteringPolys2 <- lapply(winteringPos14, toPoly, projection.in = 4326, projection.out = 'ESRI:54027') # Simulate capture and non-breeding locations of 100 individuals # nSample14.2 <- 100 # Total number sampled per simulation # How many sampled from each natural population (sampling scenarios separate # from definition scenarios) sampleBreeding14.2 <- list(round(breedingRelN14[[1]]*nSample14.2), c(rep(0,22),round(breedingRelN14[[5]][1]*nSample14.2), rep(0, 4),round(breedingRelN14[[5]][2]*nSample14.2), rep(0,44),round(breedingRelN14[[5]][3]*nSample14.2), rep(0, 4),round(breedingRelN14[[5]][4]*nSample14.2), rep(0, 22)), rep(c(rep(0, 4), rep(round(breedingRelN14[[2]][1]*nSample14.2/2), 2), rep(0, 4)), 10)) lapply(sampleBreeding14.2, matrix, nrow=10, ncol=10) # Number of simulations to run nSims14.2 <- 10 nSimsLarge14.2 <- 25 nYears <- 1 nMonths <- 1 # Set up data structures for storing results #making an empty list to fill animalLoc14 <- vector("list", nScenarios14) sim14.2 <- vector("list", nSimsLarge14.2) compare14.2.array <- array(NA, c(nSimsLarge14.2, 2), dimnames = list(1:nSimsLarge14.2, c("MC", "Mantel"))) results14 <- vector("list", nScenarios14) results14.2 <- vector("list", nSimsLarge14.2) # Run simulations set.seed(7) system.time(for (sim in 1:nSimsLarge14.2) { cat("Simulation", sim, "of", nSimsLarge14.2, '\n') sim14.2[[sim]] <- lapply(sampleBreeding14.2, simMove, breedingDist = breedDist14[[1]], winteringDist=nonbreedDist14[[1]], psi=psi, nYears=nYears, nMonths=nMonths) for (i in 13) { # only run one scenario animalLoc14.0 <- sim14.2[[sim]][[scenarioToSampleMap14[i]]]$animalLoc animalLoc14[[i]] <- changeLocations(animalLoc14.0, breedingSiteTrans14[[i]], winteringSiteTrans14[[i]]) breedingTruePoints14 <- capLocs14[[i]][animalLoc14[[i]][,1,1,1], ] winteringTruePoints14 <- targLocs14.2[[i]][animalLoc14.0[,2,1,1], ] results14.2[[sim]] <- simLocationError(winteringTruePoints14, WinteringPolys[[i]], geoBias, geoVCov, projection = sf::st_crs(winteringTruePoints14)) animalLoc14.2 <- animalLoc14[[i]] animalLoc14.2[ , 2, 1, 1] <- results14.2[[sim]]$targetSample results14[[i]] <- calcPsiMC(breedDist14[[i]], nonbreedDist14[[i]], animalLoc14.2, originRelAbund = breedingRelN14[[i]], verbose = FALSE) compare14.2.array[sim, 'MC'] <- results14[[i]]$MC originDist1 <- distFromPos(sf::st_coordinates(breedingTruePoints14)) target.point.sample <- sf::st_as_sf( data.frame(results14.2[[sim]]$targetPointSample), coords = c(1,2), crs = 6933) target.point.sample2 <- sf::st_transform(target.point.sample, 4326) targetDist1 <-distFromPos(sf::st_coordinates(target.point.sample2)) compare14.2.array[sim, 'Mantel'] <- calcMantel(originDist = originDist1, targetDist = targetDist1)$pointCorr } }) means14.2 <- apply(compare14.2.array, 2, mean) vars14.2 <- apply(compare14.2.array, 2, var) rmse14.2 <- apply(compare14.2.array, 2, function(x) sqrt(mean((x - MC14)^2))) # Set up data structures for storing estimation results est14.2.array <- array(NA, c(nSims14.2, 2), dimnames = list(1:nSims14.2, c("MC", "Mantel"))) var14.2.array <- array(NA, c(nSims14.2, 2), dimnames = list(1:nSims14.2, c("MC", "Mantel"))) ci14.2.array <- array(NA, c(nSims14.2, 2, 2), dimnames = list(1:nSims14.2, c("MC", "Mantel"), c('lower', 'upper'))) animalLoc14 <- vector("list", nSims14.2) #making an empty list to fill results14 <- vector("list", nSims14.2) # Run estimations set.seed(567) sim14.2.sub <- sim14.2[sample.int(nSimsLarge14.2, nSims14.2, TRUE)] for (sim in 1:nSims14.2) { cat("Estimation", sim, "of", nSims14.2, '\n') for (i in 13) { points0 <- sf::st_as_sf(data.frame(results14.2[[sim]]$targetPointSample), coords=c("X1","X2"), crs = 6933) points1 <- sf::st_transform(points0, crs = "ESRI:54027") animalLoc14[[sim]] <- changeLocations( sim14.2.sub[[sim]][[scenarioToSampleMap14[i]]]$animalLoc, breedingSiteTrans14[[i]], winteringSiteTrans14[[i]]) results14[[sim]] <- estMC(originDist = breedDist14[[i]], targetDist = nonbreedDist14[[i]], originRelAbund = breedingRelN14[[i]], originPoints = capLocs14[[i]][animalLoc14[[sim]][, 1, 1, 1], ], targetPoints = points1, originAssignment = animalLoc14[[sim]][, 1, 1, 1], targetAssignment = results14.2[[sim]]$targetSample, nSamples = 1000, verbose = 0, calcCorr = TRUE, geoBias = geoBias, geoVCov = geoVCov, isGL = TRUE, targetSites = winteringPolys2[[i]], maintainLegacyOutput = TRUE) est14.2.array[sim, 'MC'] <- results14[[sim]]$meanMC est14.2.array[sim, 'Mantel'] <- results14[[sim]]$meanCorr var14.2.array[sim, 'MC'] <- results14[[sim]]$seMC ^ 2 var14.2.array[sim, 'Mantel'] <- results14[[sim]]$seCorr ^ 2 ci14.2.array[sim, 'MC', ] <- results14[[sim]]$bcCI ci14.2.array[sim, 'Mantel', ] <- results14[[sim]]$bcCICorr } } pointEsts14.2 <- t(sapply(results14, function(x) c(x$pointMC, x$pointCorr))) # Summarize estimation results summary(var14.2.array) vars14.2 summary(est14.2.array) summary(pointEsts14.2) colMeans(pointEsts14.2) - MC14 sqrt(colMeans((pointEsts14.2 - MC14)^2)) means14.2 mean(ci14.2.array[, "MC", 'lower'] <= MC14 & ci14.2.array[, "MC", 'upper'] >= MC14) mean(ci14.2.array[, "Mantel", 'lower'] <= MC14 & ci14.2.array[, "Mantel", 'upper'] >= MC14) sqrt(vars14.2) / MC14 summary(sqrt(var14.2.array)/est14.2.array) # Plot histograms of bootstrap estimation results est14.2.df <- data.frame(Parameter = rep(c("MC", 'rM'), 2, each = nSims14.2), Quantity = rep(c('Mean', 'Variance'), each = 2 * nSims14.2), sim = rep(1:nSims14.2, 4), Estimate = c(est14.2.array, var14.2.array)) trues14.2.df <- data.frame(Parameter = rep(c("MC", 'rM'), 2), Quantity = rep(c('Mean', 'Variance'), each = 2), Value = c(MC14, MC14, vars14.2)) g.est <- ggplot(est14.2.df, aes(Estimate)) + geom_histogram(bins = 15) + facet_grid(Parameter ~ Quantity, scales = 'free_x') + geom_vline(aes(xintercept = Value), data = trues14.2.df) + theme_bw() + scale_x_continuous(breaks = c(0.002, 0.003, 0.1, 0.2, 0.3)) g.est # Summary table of bootstrap estimation results qualities14.2 <- data.frame(Parameter = rep(c("MC", 'rM'), 5), Quantity = rep(rep(c('Mean', 'Variance'), c(2, 2)), 3, 10), Measure = rep(c('Bias', 'RMSE', "Coverage"), c(4, 4, 2)), Value = c(colMeans(est14.2.array - MC14), colMeans(var14.2.array) - vars14.2, sqrt(colMeans((est14.2.array - MC14)^2)), sqrt(mean((var14.2.array[,1] - vars14[1])^2)), sqrt(mean((var14.2.array[,2] - vars14[2])^2)), mean(ci14.2.array[,"MC",'lower']<=MC14 & ci14.2.array[,"MC",'upper']>=MC14), mean(ci14.2.array[, "Mantel", 'lower'] <= MC14 & ci14.2.array[, "Mantel", 'upper'] >= MC14))) format(qualities14.2, digits = 2, scientific = FALSE) # Make summary table and figures of results with (Example 9) and without # (Example 8) location error qualities.all <- rbind(data.frame(DataType = 'GPS', qualities14), data.frame(DataType = 'GL', qualities14.2)) est14.all.df <- rbind(data.frame(DataType = 'GPS', est.df), data.frame(DataType = 'GL', est14.2.df)) trues14.all.df <- data.frame(DataType = rep(c("GPS", "GL"), 2, each = 2), Parameter = rep(c("MC", 'rM'), 4), Quantity = rep(c('Mean', 'Variance'), each = 4), Value = c(rep(MC14, 4), vars14, vars14.2)) g.est <- ggplot(est14.all.df, aes(Estimate)) + geom_histogram(bins = 12) + facet_grid(DataType + Parameter ~ Quantity, scales = 'free_x') + geom_vline(aes(xintercept = Value), data = trues14.all.df) + theme_bw() + scale_x_continuous(breaks = c(0.002, 0.003, 0.1, 0.2, 0.3)) g.est ``` ###EXAMPLE 10 `estMC` Estimate strength of migratory connectivity incorporating uncertainty in estimates of breeding abundance (see estMC help file) ###EXAMPLE 11 `estMC` Estimate strength of migratory connectivity incorporating detection heterogeneity (see estMC help file) ###EXAMPLE 12 `estMC` Estimate strength of migratory connectivity incorporating location uncertainty Estimates of MC may also be influenced by error in individual assignment to regions. Data types vary in location accuracy and precision, which are likely to influence the accuracy with which individuals are assigned to regions. We measured MC for bootstrapped data of birds tracked from breeding to non-breeding regions using light-level and GPS geolocation when 1) GPS location uncertainty was applied to all individuals, 2) GPS and light-level location uncertainty were applied to individuals with those devices, and 3) light-level location uncertainty was applied to all individuals. See Example 2 above and estMC help file. ###EXAMPLE 13 (from process error simulation in Cohen et al. 2018) `simMove` for measuring the influence of dispersal rates on MC Long-distance dispersal occurs when individuals that originate or breed in one population do not return to the same population to breed the next year. To quantify the sensitivity of MC to dispersal, we used simulations in which dispersal probability between breeding regions varied from low to high to measure the sensitivity of MC to dispersal. ```{r, eval = FALSE} ## Simulation ---- nBreeding <- 100 nWintering <- 100 breedingPos <- matrix(c(rep(seq(-99, -81, 2), each=sqrt(nBreeding)), rep(seq(49, 31, -2), sqrt(nBreeding))), nBreeding, 2) winteringPos <- matrix(c(rep(seq(-79, -61, 2), each=sqrt(nWintering)), rep(seq(9, -9, -2), sqrt(nWintering))), nWintering, 2) head(breedingPos) tail(breedingPos) head(winteringPos) tail(winteringPos) breedDist <- distFromPos(breedingPos, 'ellipsoid') nonbreedDist <- distFromPos(winteringPos, 'ellipsoid') # Breeding Abundance breedingN <- rep(5000, nBreeding) breedingRelN <- breedingN/sum(breedingN) # Set up psi matrix o <- optimize(mlogitMC, MC.in = 0.25, origin.dist = breedDist, target.dist = nonbreedDist, origin.abund = breedingRelN, sample.size = sum(breedingN), interval = c(0, 10), tol = .Machine$double.eps^0.5) slope <- o$minimum psi <- mlogitMat(slope, breedDist) # Baseline strength of migratory connectivity MC <- calcMC(breedDist, nonbreedDist, breedingRelN, psi, sum(breedingN)) MC # Other basic simulation parameters ## Dispersal simulations--- set.seed(1516) nYears <- 15 nMonths <- 4 # Each season Drates <- c(0.02, 0.04, 0.08, 0.16, 0.32, 0.64) #rates of dispersal birdLocDisp <- vector('list', length(Drates)) Disp.df <- data.frame(Year=rep(1:nYears, length(Drates)), Rate=rep(Drates, each = nYears), MC = NA) for(i in 1:length(Drates)){ cat('Dispersal Rate', Drates[i], '\n') birdLocDisp[[i]] <- simMove(breedingN, breedDist, nonbreedDist, psi, nYears, nMonths, sumDispRate = Drates[i]) for(j in 1:nYears){ cat('\tYear', j, '\n') temp.results <- calcPsiMC(breedDist, nonbreedDist, breedingRelN, birdLocDisp[[i]]$animalLoc, years=j, verbose = FALSE) Disp.df$MC[j + (i - 1) * nYears] <- temp.results$MC } } # end i loop Disp.df$Year <- Disp.df$Year - 1 #just run once! data.frame(Disp.df, roundMC = round(Disp.df$MC, 2), nearZero = Disp.df$MC < 0.01) # Convert dispersal rates to probabilities of dispersing at least certain distance threshold <- 1000 probFarDisp <- matrix(NA, nBreeding, length(Drates), dimnames = list(NULL, Drates)) for (i in 1:length(Drates)) { for (k in 1:nBreeding) { probFarDisp[k, i] <- sum(birdLocDisp[[i]]$natalDispMat[k, which(breedDist[k, ]>= threshold)]) } } summary(probFarDisp) ```