View source: R/estConnectivity.R
estTransition  R Documentation 
Estimation and resampling of uncertainty for psi (transition probabilities between origin sites in one phase of the annual cycle and target sites in another for migratory animals). Data can be from any combination of geolocators (GL), telemetry/GPS, intrinsic markers such as isotopes and genetics, and band/ring reencounter data.
estTransition(
originSites = NULL,
targetSites = NULL,
originPoints = NULL,
targetPoints = NULL,
originAssignment = NULL,
targetAssignment = NULL,
originNames = NULL,
targetNames = NULL,
nSamples = 1000,
isGL = FALSE,
isTelemetry = FALSE,
isRaster = FALSE,
isProb = FALSE,
captured = "origin",
geoBias = NULL,
geoVCov = NULL,
geoBiasOrigin = geoBias,
geoVCovOrigin = geoVCov,
targetRaster = NULL,
originRaster = NULL,
banded = NULL,
reencountered = NULL,
verbose = 0,
alpha = 0.05,
resampleProjection = "ESRI:102010",
nSim = ifelse(any(isRaster & isGL)  any(isRaster & isProb)  any(isGL & isProb),
5000, ifelse(any(isGL), 1000, ifelse(any(isRaster), 10, 1))),
maxTries = 300,
nBurnin = 5000,
nChains = 3,
nThin = 1,
dataOverlapSetting = c("dummy", "none", "named"),
fixedZero = NULL,
targetRelAbund = NULL,
method = c("bootstrap", "MCMC", "moutofnbootstrap"),
m = NULL,
psiPrior = NULL,
returnAllInput = TRUE
)
estPsi(
originSites = NULL,
targetSites = NULL,
originPoints = NULL,
targetPoints = NULL,
originAssignment = NULL,
targetAssignment = NULL,
originNames = NULL,
targetNames = NULL,
nSamples = 1000,
isGL = FALSE,
isTelemetry = FALSE,
isRaster = FALSE,
isProb = FALSE,
captured = "origin",
geoBias = NULL,
geoVCov = NULL,
geoBiasOrigin = geoBias,
geoVCovOrigin = geoVCov,
targetRaster = NULL,
originRaster = NULL,
banded = NULL,
reencountered = NULL,
verbose = 0,
alpha = 0.05,
resampleProjection = "ESRI:102010",
nSim = ifelse(any(isRaster & isGL)  any(isRaster & isProb)  any(isGL & isProb),
5000, ifelse(any(isGL), 1000, ifelse(any(isRaster), 10, 1))),
maxTries = 300,
nBurnin = 5000,
nChains = 3,
nThin = 1,
dataOverlapSetting = c("dummy", "none", "named"),
fixedZero = NULL,
targetRelAbund = NULL,
method = c("bootstrap", "MCMC", "moutofnbootstrap"),
m = NULL,
psiPrior = NULL,
returnAllInput = TRUE
)
originSites 
A polygon spatial layer (sf  MULTIPOLYGON) defining the geographic representation of sites in the origin season. 
targetSites 
A polygon spatial layer (sf  MULTIPOLYGON) defining the geographic representation of sites in the target season. 
originPoints 
A 
targetPoints 
For GL or telemetry data, a 
originAssignment 
Assignment of animals to origin season sites. Either
an integer vector with length number of animals tracked or a matrix of
probabilities with number of animals tracked rows and number of origin sites
columns (and rows summing to 1). The latter only applies to animals released
in the target sites where there is uncertainty about their origin site, for
example from genetic population estimates from the rubias package.
Optional, but some combination of these inputs should be defined. Note that
if 
targetAssignment 
Assignment of animals to target season sites. Either
an integer vector with length number of animals tracked or a matrix of
probabilities with number of animals tracked rows and number of target sites
columns (and rows summing to 1). The latter only applies to animals released
in the origin sites where there is uncertainty about their target site, for
example from genetic population estimates from the rubias package.
Optional, but some combination of these inputs needs to be defined. Note
that if 
originNames 
Optional, but recommended to keep track. Vector of names for the origin sites. If not provided, the function will either try to get these from another input or provide default names (capital letters). 
targetNames 
Optional, but recommended to keep track. Vector of names for the target sites. If not provided, the function will either try to get these from another input or provide default names (numbers). 
nSamples 
Number of postburnin MCMC samples to store ( 
isGL 
Indicates whether or which animals were tracked with geolocators.
Should be either single TRUE or FALSE value, or vector with length of
number of animals tracked, with TRUE or FALSE for each animal in data
(except those in 
isTelemetry 
Indicates whether or which animals were tracked with
telemetry/GPS (no location uncertainty on either end).
Should be either single TRUE or FALSE value, or vector with length of
number of animals tracked, with TRUE or FALSE for each animal in data
(except those in 
isRaster 
Indicates whether or which animals were tracked with
intrinsic markers (e.g., genetics or isotopes), with location uncertainty
expressed as a raster of probabilities by grid cells, either in

isProb 
Indicates whether or which animals were tracked with
intrinsic markers (e.g., genetics or isotopes), with location uncertainty
expressed as a probability table, either in 
captured 
Indicates whether or which animals were captured in the origin sites, the target sites, or neither (another phase of the annual cycle). Location uncertainty will only be applied where the animal was not captured. So this doesn't matter for telemetry data, and is assumed to be "origin" for band return data. Should be either single "origin" (default), "target", or "neither" value, or a character vector with length of number of animals tracked, with "origin", "target", or "neither" for each animal. 
geoBias 
For GL data, vector of length 2 indicating expected bias
in longitude and latitude of 
geoVCov 
For GL data, 2x2 matrix with expected variance/covariance
in longitude and latitude of 
geoBiasOrigin 
For GL data where 
geoVCovOrigin 
For GL data where 
targetRaster 
For intrinsic tracking data, the results of

originRaster 
For intrinsic tracking data, the results of

banded 
For band return data, a vector or matrix of the number of
released animals from each origin site (including those never reencountered
in a target site). If a matrix, the second dimension is taken as the number
of age classes of released animals; the model estimates reencounter
probability by age class but assumes transition probabilities are the same.
Note that this age model is currently implemented only for 
reencountered 
For band return data, either a matrix with B rows and W columns or a B x [number of ages] x W array. Number of animals reencountered on each target site (by age class banded as) by origin site they came from. 
verbose 
0 (default) to 3. 0 prints no output during run (except on
convergence for 
alpha 
Level for confidence/credible intervals provided. Default (0.05) gives 95 percent CI. 
resampleProjection 
Projection when sampling from location uncertainty. Default is Equidistant Conic. The default setting preserves distances around latitude = 0 and longitude = 0. Other projections may work well, depending on the location of sites. Ignored unless data are entered using sites and points and/or rasters. 
nSim 
Tuning parameter for GL or intrinsic data. Affects only the
speed; 1000 seems to work well with our GL data and 10 for our intrinsic
data, but your results may vary. For data combinations, we put the default
higher (5000) to allow for more data conflicts. Should be integer > 0.
Ignored when 
maxTries 
Maximum number of times to run a single GL/intrinsic
bootstrap before exiting with an error. Default is 300; you may want to make
a little higher if your 
nBurnin 
For 
nChains 
For 
nThin 
For 
dataOverlapSetting 
When there is more than one type of data, this
setting allows the user some flexibility for clarifying which type(s) of
data apply to which animals. Setting "dummy" (the default) indicates that
there are dummy values within each dataset for the animals that isGL,
isTelemetry, etc. don't have that data type (FALSE values). If no animals
have a data type, no dummy values are required. If no animals have more than
one type of data, the user can simplify processing their data by choosing
setting "none" here. In this case, there should be no dummy values, and only
the animals with a type of data should be included in that dataset. The
third setting ("named") is not yet implemented, but will eventually allow
another way to allow animals with more than one type of data with named
animals linking records. When there is only one type of data, it is fastest
to leave this on the default. Note that banding data entered through

fixedZero 
When the user has a priori reasons to believe one or more transition probabilities are zero, they can indicate those here, and the model will keep them fixed at zero. This argument should be a matrix with two columns (for row and column of the transition probability matrix) and number of transitions being fixed to zero rows. For MCMC modeling, substantial evidence that a transition fixed to zero isn't zero may cause an error. For bootstrap modeling, a warning will come up if any bootstrap runs generate the transition fixed to zero, and the function will quit with an error if a very large number of runs do (> 10 * nSamples). Fixing transitions to zero may also slow down the bootstrap model somewhat. 
targetRelAbund 
When some/all data have location error at origin sites
(i.e., GL, raster, or probability table data with captured = "target" or
"none"), unless the data were collected in proportion to abundance at target
sites, simulation work indicates substantial bias in transition probability
estimates can result. However, if these data are resampled in proportion to
target site abundance, this bias is removed. This argument allows the user
to provide an estimate of relative abundance at the target sites. Either
a numeric vector of length [number target sites] that sums to 1, or an mcmc
object (such as is produced by 
method 
This important setting lets the user choose the estimation
method used: bootstrap or MCMC (Markov chain Monte Carlo). Bootstrap (the
default) now works with any and all types of data, whereas MCMC currently
only works with banding and telemetry data (enter telemetry data for MCMC
using 
m 
We read that the moutofnbootstrap method may improve the
coverage of confidence intervals for parameters on or near a boundary (0 or
1 in this case). So we're testing that out. This still under development and
not for the end user. In the moutofnbootstrap, m is the number of
samples taken each time (less than the true sample size, n). If the
"moutofnbootstrap" is chosen under 
psiPrior 
matrix with same dimensions as psi. Only relevant when

returnAllInput 
if TRUE (the default) the output includes all of the inputs. If FALSE, only the inputs currently used by another MigConnectivity function are included in the output. Switch this if you're worried about computer memory (and the output will be much slimmer). 
estTransition
returns a list with the elements:
psi
List containing estimates of transition probabilities:
sample
Array of sampled values for psi. nSamples
x
[number of origin sites] x [number of target sites]. Provided to allow
the user to compute own summary statistics.
mean
Main estimate of psi matrix. [number of origin sites]
x [number of target sites].
se
Standard error of psi, estimated from SD of
psi$sample
.
simpleCI
1  alpha
confidence interval for psi,
estimated as alpha/2
and 1  alpha/2
quantiles of
psi$sample
.
bcCI
Biascorrected 1  alpha
confidence interval
for psi. May be preferable to simpleCI
when mean
is the
best estimate of psi. simpleCI
is preferred when
median
is a better estimator. When the mean and median are equal,
these should be identical. Estimated as the
pnorm(2 * z0 + qnorm(alpha / 2))
and
pnorm(2 * z0 + qnorm(1  alpha / 2))
quantiles of sample
,
where z0 is the proportion of sample < mean
.
hpdCI
1  alpha
credible interval for psi,
estimated using the highest posterior density (HPD) method.
median
Median estimate of psi matrix.
point
Simple point estimate of psi matrix, not accounting
for sampling error.
r
List containing estimates of reencounter probabilities at each target site. NULL except when using direct band/ring reencounter data.
input
List containing the inputs to estTransition
.
BUGSoutput
List containing R2jags
output. Only present
when using method
of "MCMC".
estStrength
, plot.estMigConnectivity
,
estMC
, estMantel
##############################################################################
# Examples 1 (banding data: first example is based on common tern banding
# data; the second is made up data to demonstrate data with two ages)
##############################################################################
COTE_banded < c(10360, 1787, 2495, 336)
COTE_reencountered < matrix(c(12, 0, 38, 15,
111, 7, 6, 2,
5, 0, 19, 4,
1123, 40, 41, 7),
4, 4,
dimnames = list(LETTERS[1:4], 1:4))
COTE_psi < estTransition(originNames = LETTERS[1:4],
targetNames = 1:4,
banded = COTE_banded,
reencountered = COTE_reencountered,
verbose = 1,
nSamples = 60000, nBurnin = 20000,
method = "MCMC")
COTE_psi
COTE_banded2 < matrix(rep(COTE_banded, 2), 4, 2)
COTE_reencountered2 < array(c(12, 0, 38, 15, 6, 0, 17, 7,
111, 7, 6, 2, 55, 3, 3, 1,
5, 0, 19, 4, 2, 0, 10, 2,
1123, 40, 41, 7, 660, 20, 20, 3),
c(4, 2, 4),
dimnames = list(LETTERS[1:4], c("J", "A"), 1:4))
COTE_psi2 < estTransition(originNames = LETTERS[1:4],
targetNames = 1:4,
banded = COTE_banded2,
reencountered = COTE_reencountered2,
verbose = 0,
nSamples = 60000, nBurnin = 20000,
method = "MCMC")
COTE_psi2
##############################################################################
# Example 2 (geolocator and telemetry ovenbirds captured on origin sites)
##############################################################################
data(OVENdata) # Ovenbird
nSamplesGLGPS < 100 # Number of bootstrap iterations
# Estimate transition probabilities; treat all data as geolocator
GL_psi < estTransition(isGL=TRUE,
geoBias = OVENdata$geo.bias,
geoVCov = OVENdata$geo.vcov,
targetSites = OVENdata$targetSites,
originSites = OVENdata$originSites,
originPoints = OVENdata$originPoints,
targetPoints = OVENdata$targetPoints,
verbose = 2,
nSamples = nSamplesGLGPS,
resampleProjection=sf::st_crs(OVENdata$targetPoints))
# Treat all data as is
Combined.psi < estTransition(isGL=OVENdata$isGL,
isTelemetry = !OVENdata$isGL,
geoBias = OVENdata$geo.bias, # Lightlevel GL location bias
geoVCov = OVENdata$geo.vcov, # Location covariance matrix
targetSites = OVENdata$targetSites, # Nonbreeding/target sites
originSites = OVENdata$originSites, # Breeding/origin sites
originPoints = OVENdata$originPoints, # Capture Locations
targetPoints = OVENdata$targetPoints, #Device target locations
verbose = 2, # output options
nSamples = nSamplesGLGPS, # This is set low for example
resampleProjection = sf::st_crs(OVENdata$targetPoints))
print(Combined.psi)
# For treating all data as GPS,
# Move the latitude of birds with locations that fall offshore
int < sf::st_intersects(OVENdata$targetPoints, OVENdata$targetSites)
any(lengths(int)<1)
plot(OVENdata$targetPoints)
plot(OVENdata$targetSites,add=TRUE)
tp<sf::st_coordinates(OVENdata$targetPoints)
text(tp[,1], tp[,2], label=c(1:39))
tp[5,2] < 2450000
tp[10,2]< 2240496
tp[1,2]< 2240496
tp[11,2]< 2026511
tp[15,2]< 2031268
tp[16,2]< 2031268
oven_targetPoints<sf::st_as_sf(as.data.frame(tp),
coords = c("X","Y"),
crs = sf::st_crs(OVENdata$targetPoints))
inter < sf::st_intersects(oven_targetPoints, OVENdata$targetSites)
any(lengths(inter)<1)
plot(oven_targetPoints,add=TRUE, col = "green")
plot(oven_targetPoints[lengths(inter)<1,],add=TRUE, col = "darkblue")
# Treat all data as GPS
GPS_psi < estTransition(isTelemetry = TRUE,
targetSites = OVENdata$targetSites, # Nonbreeding/target sites
originSites = OVENdata$originSites, # Breeding/origin sites
originPoints = OVENdata$originPoints, # Capture Locations
targetPoints = oven_targetPoints, # Device target locations
verbose = 2, # output options
nSamples = nSamplesGLGPS) # This is set low for example
##############################################################################
# Example 3 (all released origin; some telemetry, some GL, some probability
# tables, some both GL and probability tables; data modified from ovenbird
# example)
##############################################################################
library(VGAM)
nAnimals < 40
isGL < c(OVENdata$isGL, FALSE)
isTelemetry < c(!OVENdata$isGL, FALSE)
isRaster < rep(FALSE, nAnimals)
isProb < rep(FALSE, nAnimals)
targetPoints < rbind(OVENdata$targetPoints, OVENdata$targetPoints[1,])
targetSites < OVENdata$targetSites
originSites < OVENdata$originSites
resampleProjection < sf::st_crs(OVENdata$targetPoints)
targetNames < OVENdata$targetNames
originNames < OVENdata$originNames
targetAssignment < array(0, dim = c(nAnimals, 3),
dimnames = list(NULL, targetNames))
assignment0 < unclass(sf::st_intersects(x = targetPoints, y = targetSites,
sparse = TRUE))
assignment0[sapply(assignment0, function(x) length(x)==0)] < 0
assignment0 < array(unlist(assignment0), nAnimals)
for (ani in 1:nAnimals) {
if (assignment0[ani]>0)
targetAssignment[ani, assignment0[ani]] < 1
else{
targetAssignment[ani, ] < rdiric(1, c(15, 1, 1))
isProb[ani] < TRUE
}
}
targetAssignment
isProb
nSamplesTry < 100 # Number of bootstrap iterations
originPoints < rbind(OVENdata$originPoints,
OVENdata$originPoints[39,])
system.time(psi3 <
estTransition(isGL = isGL, isRaster = isRaster,
isProb = isProb,
isTelemetry = isTelemetry,
geoBias = OVENdata$geo.bias,
geoVCov = OVENdata$geo.vcov,
targetPoints = targetPoints,
targetAssignment = targetAssignment,
targetSites = targetSites,
resampleProjection = resampleProjection,
nSim = 20000, maxTries = 300,
originSites = originSites,
originPoints = originPoints,
captured = "origin",
originNames = OVENdata$originNames,
targetNames = OVENdata$targetNames,
verbose = 3,
nSamples = nSamplesTry))
psi3
nNonBreeding < nrow(OVENdata$targetSites)
plot(psi3, legend = "top",
main = paste("OVENlike w/", sum(isGL & !isProb), "GL,",
sum(!isGL & isProb), "probs,",
sum(isGL & isProb), "both, and", sum(isTelemetry), "GPS"))
##############################################################################
# Example 4 (add probability animals released on other end)
##############################################################################
nAnimals < 45
captured < rep(c("origin", "target"), c(40, 5))
isGL < c(OVENdata$isGL, rep(FALSE, 6))
isTelemetry < c(!OVENdata$isGL, rep(FALSE, 6))
isRaster < rep(FALSE, nAnimals)
isProb < rep(FALSE, nAnimals)
targetPoints < rbind(OVENdata$targetPoints,
OVENdata$targetPoints[c(1:3,19,23,31),])
targetAssignment < array(0, dim = c(nAnimals, 3),
dimnames = list(NULL, targetNames))
assignment0 < unclass(sf::st_intersects(x = targetPoints, y = targetSites,
sparse = TRUE))
assignment0[sapply(assignment0, function(x) length(x)==0)] < 0
assignment0 < array(unlist(assignment0), nAnimals)
for (ani in 1:nAnimals) {
if (assignment0[ani]>0)
targetAssignment[ani, assignment0[ani]] < 1
else{
targetAssignment[ani, ] < rdiric(1, c(15, 1, 1))
isProb[ani] < TRUE
}
}
targetAssignment
isProb
originPoints < rbind(OVENdata$originPoints,
OVENdata$originPoints[34:39,])
originPoints < sf::st_transform(originPoints, crs = resampleProjection)
originSites < sf::st_transform(OVENdata$originSites,
crs = resampleProjection)
assignment1 < unclass(sf::st_intersects(x = originPoints, y = originSites,
sparse = TRUE))
assignment1[sapply(assignment1, function(x) length(x)==0)] < 0
assignment1 < array(unlist(assignment1), nAnimals)
nOriginSites < nrow(originSites)
originAssignment < array(0, dim = c(nAnimals, nOriginSites),
dimnames = list(NULL, originNames))
for (ani in 1:40) {
originAssignment[ani, assignment1[ani]] < 1
}
for (ani in 41:nAnimals) {
originAssignment[ani, ] < rdiric(1, c(1, 1))
isProb[ani] < TRUE
}
originAssignment
isProb
system.time(psi4 <
estTransition(isGL = isGL, isRaster = isRaster,
isProb = isProb,
isTelemetry = isTelemetry,
geoBias = OVENdata$geo.bias,
geoVCov = OVENdata$geo.vcov,
targetPoints = targetPoints,
targetAssignment = targetAssignment,
targetSites = targetSites,
resampleProjection = resampleProjection,
nSim = 15000, maxTries = 300,
originSites = originSites,
originAssignment = originAssignment,
captured = captured,
originNames = OVENdata$originNames,
targetNames = OVENdata$targetNames,
verbose = 2,
nSamples = nSamplesTry,
targetRelAbund = c(0.1432, 0.3577, 0.4991)))
psi4
plot(psi4, legend = "top",
main = paste(sum(isGL & !isProb), "GL,",
sum(!isGL & isProb & captured == "origin"), "prob.,",
sum(isGL & isProb), "both,",
sum(isTelemetry), "GPS (all\ncaptured origin), and",
sum(isProb & captured == "target"),
"prob. (captured target)"))
MC4 < estStrength(OVENdata$originDist, OVENdata$targetDist,
OVENdata$originRelAbund, psi4,
sampleSize = nAnimals)
MC4
##############################################################################
# Example 5 (all raster, from our OVEN example)
##############################################################################
getCSV < function(filename) {
tmp < tempdir()
url1 < paste0(
'https://github.com/SMBCNZP/MigConnectivity/blob/master/dataraw/',
filename, '?raw=true')
temp < paste(tmp, filename, sep = '/')
utils::download.file(url1, temp, mode = 'wb')
csv < read.csv(temp)
unlink(temp)
return(csv)
}
getRDS < function(speciesDist) {
tmp < tempdir()
extension < '.rds'
filename < paste0(speciesDist, extension)
url1 < paste0(
'https://github.com/SMBCNZP/MigConnectivity/blob/master/dataraw/Spatial_Layers/',
filename, '?raw=true')
temp < paste(tmp, filename, sep = '/')
utils::download.file(url1, temp, mode = 'wb')
shp < readRDS(temp)
unlink(temp)
return(shp)
}
OVENdist < getRDS("OVENdist")
OVENdist < sf::st_as_sf(OVENdist)
OVENdist < sf::st_transform(OVENdist, 4326)
OVENvals < getCSV("deltaDvalues.csv")
OVENvals < OVENvals[grep(x=OVENvals$Sample,"NH", invert = TRUE),]
originSites < getRDS("originSites")
originSites < sf::st_as_sf(originSites)
EVER < length(grep(x=OVENvals$Sample,"EVER"))
JAM < length(grep(x=OVENvals$Sample,"JAM"))
originRelAbund < matrix(c(EVER,JAM),nrow = 1,byrow = TRUE)
originRelAbund < prop.table(originRelAbund,1)
op < sf::st_centroid(originSites)
originPoints < array(NA,c(EVER+JAM,2), list(NULL, c("x","y")))
originPoints[grep(x = OVENvals$Sample,"JAM"),1] < sf::st_coordinates(op)[1,1]
originPoints[grep(x = OVENvals$Sample,"JAM"),2] < sf::st_coordinates(op)[1,2]
originPoints[grep(x = OVENvals$Sample,"EVER"),1]<sf::st_coordinates(op)[2,1]
originPoints[grep(x = OVENvals$Sample,"EVER"),2]<sf::st_coordinates(op)[2,2]
originPoints < sf::st_as_sf(data.frame(originPoints),
coords = c("x", "y"),
crs = sf::st_crs(originSites))
iso < isoAssign(isovalues = OVENvals[,2],
isoSTD = 12, # this value is for demonstration only
intercept = 10, # this value is for demonstration only
slope = 0.8, # this value is for demonstration only
odds = NULL,
restrict2Likely = FALSE,
nSamples = 1000,
sppShapefile = terra::vect(OVENdist),
assignExtent = c(179,60,15,89),
element = "Hydrogen",
period = "GrowingSeason",#this setting for demonstration only
seed = 12345,
verbose=1)
nAnimals < dim(iso$probassign)[3]
isGL <rep(FALSE, nAnimals); isRaster < rep(TRUE, nAnimals)
isProb < rep(FALSE, nAnimals); isTelemetry < rep(FALSE, nAnimals)
targetSites < sf::st_as_sf(iso$targetSites)
targetSites < sf::st_make_valid(targetSites)
targetSites < sf::st_union(targetSites, by_feature = TRUE)
system.time(psi5 <
estTransition(isGL = isGL,
isRaster = isRaster,
isProb = isProb,
isTelemetry = isTelemetry,
targetSites = targetSites,
resampleProjection = resampleProjection,
targetRaster = iso,
originSites = originSites,
originPoints = originPoints,
captured = rep("origin", nAnimals),
verbose = 2,
nSamples = nSamplesTry))
psi5
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.