install.packages('devtools')
devtools::install_github("SMBC-NZP/MigConnectivity", build_vignettes = TRUE)
library(MigConnectivity)
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 et
al. 2018, Roberts et al. 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
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")
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))
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)
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
#> Migratory Connectivity Estimates
#>
#> MC estimate (mean): 0.2276 +/- (SE) 0.01561
#> 95% confidence interval (simple quantile): 0.1979 - 0.2577
#>
#> This is a subset of what's available inside this estMigConnectivity output.
#> For more info, try ?estStrength or ?estMC or str(obj_name, max.level = 2).
MCest$MC$mean - MCtrue
#> [1] 0.03167019
# 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)
#> List of 22
#> $ originSites :Classes 'sf' and 'data.frame': 5 obs. of 3 variables:
#> ..- attr(*, "sf_column")= chr "geometry"
#> ..- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA
#> .. ..- attr(*, "names")= chr [1:2] "population" "originSite"
#> $ targetSites :Classes 'sf' and 'data.frame': 4 obs. of 3 variables:
#> ..- attr(*, "sf_column")= chr "geometry"
#> ..- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA
#> .. ..- attr(*, "names")= chr [1:2] "Region" "targetSite"
#> $ originPoints :Classes 'sf' and 'data.frame': 7 obs. of 3 variables:
#> ..- attr(*, "sf_column")= chr "geometry"
#> ..- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA
#> .. ..- attr(*, "names")= chr [1:2] "bird" "Type"
#> $ targetPoints :Classes 'sf' and 'data.frame': 210 obs. of 2 variables:
#> ..- attr(*, "sf_column")= chr "geometry"
#> ..- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA
#> .. ..- attr(*, "names")= chr "Type"
#> $ originAssignment : num [1:203, 1:5] 6.16e-03 1.91e-17 1.86e-18 1.38e-20 5.45e-11 ...
#> ..- attr(*, "dimnames")=List of 2
#> $ originNames : chr [1:5] "Arctic" "Pacific Northwest" "Southwest" "Central" ...
#> $ targetNames : chr [1:4] "Pacific and Central Mexico" "Atlantic Lowland Mexico" "Central America" "South America"
#> $ nSamples : num 10
#> $ isGL : logi [1:210] TRUE TRUE TRUE TRUE TRUE TRUE ...
#> $ isTelemetry : logi [1:210] FALSE FALSE FALSE FALSE FALSE FALSE ...
#> $ isRaster : logi [1:210] FALSE FALSE FALSE FALSE FALSE FALSE ...
#> $ isProb : logi [1:210] FALSE FALSE FALSE FALSE FALSE FALSE ...
#> $ captured : chr [1:210] "origin" "origin" "origin" "origin" ...
#> $ geoBias : Named num [1:2] -6513 50965
#> ..- attr(*, "names")= chr [1:2] "difflong" "difflat"
#> $ geoVCov : num [1:2, 1:2] 2.59e+08 -8.94e+08 -8.94e+08 4.71e+09
#> ..- attr(*, "dimnames")=List of 2
#> $ verbose : num 2
#> $ maxTries : num 400
#> $ resampleProjection:List of 2
#> ..- attr(*, "class")= chr "crs"
#> $ nSim : num 40
#> $ dataOverlapSetting: chr "none"
#> $ targetRelAbund : Named num [1:4] 0.0881 0.1645 0.4287 0.3187
#> ..- attr(*, "names")= chr [1:4] "Pacific and Central Mexico" "Atlantic Lowland Mexico" "Central America" "South America"
#> $ originRelAbund : Named num [1:5] 0.169 0.125 0.284 0.247 0.175
#> ..- attr(*, "names")= chr [1:5] "Arctic" "Pacific Northwest" "Southwest" "Central" ...
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)
#> Warning: st_centroid assumes attributes are constant over geometries
originCentersYEWA <- st_transform(originCentersYEWA, 4326)
originDistYEWA <- distFromPos(st_coordinates(originCentersYEWA$geometry))
# calculate distance between targetSites
targetCentersYEWA <- st_centroid(YEWAdata$targetSites)
#> Warning: st_centroid assumes attributes are constant over geometries
targetCentersYEWA <- st_transform(targetCentersYEWA, 4326)
targetDistYEWA <- distFromPos(st_coordinates(targetCentersYEWA$geometry))
# Estimate the strength of migratory connectivity
YEWAmc <- estStrength(originDist = originDistYEWA,
targetDist = targetDistYEWA,
originRelAbund = YEWAdata$originRelAbund,
psi = psiYEWA,
nSamples = 5000,
verbose = 0)
YEWAmc
#> Migratory Connectivity Estimates
#>
#> MC estimate (mean): 0.3888 +/- (SE) 0.01555
#> 95% confidence interval (simple quantile): 0.3668 - 0.4146
#>
#> This is a subset of what's available inside this estMigConnectivity output.
#> For more info, try ?estStrength or ?estMC or str(obj_name, max.level = 2).