knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
install.packages('devtools') devtools::install_github("SMBC-NZP/MigConnectivity", build_vignettes = TRUE) library(MigConnectivity)
library(MigConnectivity) library(sf) set.seed(123)
estStrength
- Estimate the strength of migratory connectivityThe estStrength
function estimates the strength of migratory connectivity (MC) between populations during two different time periods within the annual cycle [@cohen_quantifying_2018; @roberts_migratory_2023]. Below, we illustrate how to estimate MC between three breeding and non-breeding regions.
To estimate the strength of MC, you will need:
psi
) Simple example with four breeding and four nonbreeding regions and simulated data
Define the number of breeding and nonbreeding populations
nBreeding <- 4 #number of breeding regions nNonBreeding <- 4 #number of non-breeding regions
Generate distance matrices between the different regions
originPos <- matrix(c(seq(-99, -81, 6), rep(35, nBreeding)), nBreeding, 2) targetPos <- matrix(c(seq(-93, -87, 2), rep(9, nNonBreeding)), nNonBreeding, 2) #distances between centroids of four breeding regions breedDist <- distFromPos(originPos, surface = "ellipsoid") #distances between centroids of four nonbreeding regions nonBreedDist <- distFromPos(targetPos, surface = "ellipsoid")
op <- graphics::par(no.readonly = TRUE) on.exit(graphics::par(op)) par(bty="n") plot(originPos[,1],originPos[,2], ylim=c(8.5,36), xlim=c(-100, -80), pch=19, axes=FALSE, yaxt="n", xaxt="n", cex=2, ylab="", xlab="Nonbreeding", main="Breeding") points(targetPos[,1],targetPos[,2],pch=15,cex=2)
Define the transition probabilities between the breeding and nonbreeding regions and load previously estimated transition probabilities (using the package RMark)
#transition probabilities form each breeding to each nonbreeding region psiTrue <- samplePsis[["Low"]] psiEst <- getCMRexample(1) # Define total sample size for psi data # for small sample corrected version of MC n <- length(grep('[2-5]', psiEst$data$data$ch))
plot(originPos[,1],originPos[,2], ylim=c(8.5,36), xlim=c(-100, -80), pch=19, axes=FALSE, yaxt="n", xaxt="n", cex=2, ylab="", xlab="Nonbreeding", main="Breeding") points(targetPos[,1],targetPos[,2],pch=15,cex=2) for (b in 1:nBreeding) { for (nb in 1:nNonBreeding) { segments(originPos[b, 1], originPos[b, 2], targetPos[nb, 1], targetPos[nb, 2], lwd = psiTrue[b, nb]*10, lty = b) } }
Define the relative abundance within the four breeding regions, and provide
previously simulated and estimated relative abundances (simple BBS-style analyses
using the functions simCountData
and modelCountDataJAGS
)
#equal relative abundance among the four breeding regions, must sum to 1 relNTrue <- rep(1/nBreeding, nBreeding) relNEst <- abundExamples[[1]]
Calculate the true strength of migratory connectivity using the inputs above
# Calculate the strength of migratory connectivity MCtrue <- calcMC(originDist = breedDist, targetDist = nonBreedDist, originRelAbund = relNTrue, psi = psiTrue)
MCtrue
Estimating the strength of MC by explicitly defining the total sample size of animals to estimate the transition probabilities
MCest <- estStrength(originDist = breedDist, targetDist = nonBreedDist, originRelAbund = relNEst, psi = psiEst, sampleSize = n, originSites = 5:8, targetSites = c(3,2,1,4), nSamples = 1000)
MCest MCest$MC$mean - MCtrue
# Read in the processed Yellow Warbler data # see the Worked Examples Vignette for details on how data structured/created newDir <- tempdir() baseURL <- 'https://github.com/SMBC-NZP/MigConnectivity/blob/devpsi2/data-raw/YEWA/' file.name <- "yewa_estTrans_data.rds" url1 <- paste0(baseURL, file.name, '?raw=true') temp <- paste(newDir, file.name, sep = '/') utils::download.file(url1, temp, mode = 'wb') YEWA_target_sites <- readRDS(temp) YEWAdata <- readRDS(temp) unlink(temp) # take a quick look at the data # str(YEWAdata,1)
First, we estimate (psi) or transition probabilities using the estTransition
function.
Note - in the Yellow Warbler data there are no overlapping data sources for individuals, i.e., we only had one type of data for each bird - so we set the overlap setting to “none”.
## Run analysis for psi psiYEWA <- estTransition(originSites = YEWAdata$originSites, targetSites = YEWAdata$targetSites, originPoints = YEWAdata$originPoints, targetPoints = YEWAdata$targetPoints, originAssignment = YEWAdata$originAssignment, originNames = YEWAdata$originNames, targetNames = YEWAdata$targetNames, nSamples = 10, # Set low for illustration isGL = YEWAdata$isGL, isTelemetry = YEWAdata$isTelemetry, isRaster = YEWAdata$isRaster, isProb = YEWAdata$isProb, captured = YEWAdata$captured, geoBias = YEWAdata$geoBias, geoVCov = YEWAdata$geoVCov, originRaster = YEWAdata$originRaster, verbose = 0, maxTries = 400, resampleProjection = YEWAdata$resampleProjection, nSim = 30, dataOverlapSetting = "none", targetRelAbund = YEWAdata$targetRelAbund, returnAllInput = FALSE)
We then use the resulting object to estimate the strength of migratory connectivity. We first need some additional data, the distances between originSites and targetSites.
# calculate distance between originSites originCentersYEWA <- st_centroid(YEWAdata$originSites) originCentersYEWA <- st_transform(originCentersYEWA, 4326) originDistYEWA <- distFromPos(st_coordinates(originCentersYEWA$geometry)) # calculate distance between targetSites targetCentersYEWA <- st_centroid(YEWAdata$targetSites) targetCentersYEWA <- st_transform(targetCentersYEWA, 4326) targetDistYEWA <- distFromPos(st_coordinates(targetCentersYEWA$geometry))
newDir <- tempdir() baseURL <- 'https://github.com/SMBC-NZP/MigConnectivity/blob/devpsi2/data-raw/YEWA/' file.name <- "psiYEWA.rds" url1 <- paste0(baseURL, file.name, '?raw=true') temp <- paste(newDir, file.name, sep = '/') utils::download.file(url1, temp, mode = 'wb') psiYEWA <- readRDS(temp) unlink(temp)
# Estimate the strength of migratory connectivity YEWAmc <- estStrength(originDist = originDistYEWA, targetDist = targetDistYEWA, originRelAbund = YEWAdata$originRelAbund, psi = psiYEWA, nSamples = 5000, verbose = 0) YEWAmc
#saveRDS(YEWAmc, "YEWAmc.rds")
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.