estTransition: Estimate psi (transition probabilities between locations in...

View source: R/estConnectivity.R

estTransitionR Documentation

Estimate psi (transition probabilities between locations in two phases of the annual cycle)

Description

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.

Usage

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", "m-out-of-n-bootstrap"),
  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", "m-out-of-n-bootstrap"),
  m = NULL,
  psiPrior = NULL,
  returnAllInput = TRUE
)

Arguments

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 sf or SpatialPoints object, with number of rows or length being the number of animals tracked. Each point indicates the origin location of an animal (or point estimate of same, for GL animals released on target sites). Note that to simplify input of multiple data types both between and for the same animal, if origin points are provided for any animal, they must be provided for all except banding data (can be dummy values), unless dataOverlapSetting is set to "none".

targetPoints

For GL or telemetry data, a sf or SpatialPoints object, with length or number of rows number of animals tracked. Each point indicates the point estimate location of an animal in the target season. Note that to simplify input of multiple data types both between and for the same animal, if target points are provided for any animal, they must be provided for all except banding data (can be dummy values), unless dataOverlapSetting is set to "none".

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 originAssignment is a probability table, animals with known origin sites can have 1 in that column and 0s in all others. Also note that if method is "MCMC", anything in originAssignment and targetAssignment will be assumed to represent animals tracked via telemetry, with known origin and target sites.

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 targetAssignment is a probability table, animals with known target sites can have 1 in that column and 0s in all others.

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 post-burn-in MCMC samples to store (method == "MCMC") OR number of bootstrap runs for method == "bootstrap". In the latter case, animals are sampled with replacement for each. For all, the purpose is to estimate sampling uncertainty.

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 banded, which are handled separately). For TRUE animals, the model applies geoBias and geoVCov to targetPoints where captured == "origin" or "neither" and geoBiasOrigin and geoVCovOrigin to originPoints where captured == "target" or "neither". Geolocator data should be entered as originPoints and targetPoints.

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 banded, which are handled separately). Telemetry data can be entered as points or using the targetAssignment and originAssignment arguments.

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 targetRaster or originRaster. 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 banded, which are handled separately).

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 targetAssignment or originAssignment. 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 banded, which are handled separately).

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 targetPoints, in resampleProjection units (default meters).

geoVCov

For GL data, 2x2 matrix with expected variance/covariance in longitude and latitude of targetPoints, in resampleProjection units (default meters).

geoBiasOrigin

For GL data where captured!="origin", vector of length 2 indicating expected bias in longitude and latitude of originPoints, in resampleProjection units (default meters).

geoVCovOrigin

For GL data where captured!="origin", 2x2 matrix with expected variance/covariance in longitude and latitude of targetPoints, in resampleProjection units (default meters).

targetRaster

For intrinsic tracking data, the results of isoAssign or a similar function of class intrinsicAssign or class RasterBrick/RasterStack, for example from the package assignR. In any case, it expresses location uncertainty on target range, through a raster of probabilities by grid cells.

originRaster

For intrinsic tracking data, the results of isoAssign or a similar function of class intrinsicAssign or class RasterBrick/RasterStack, for example from the package assignR. In any case, it expresses location uncertainty on origin range, through a raster of probabilities by grid cells.

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 method set to "MCMC", and only when banding data is analyzed alone (no telemetry data).

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 method set to "MCMC"). 1 prints an update every 100 samples or bootstraps (or a status bar for "MCMC"). 2 prints an update every sample or bootstrap. 3 also prints the number of draws (for tuning nSim).

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 method is "MCMC".

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 nSim is low and nSamples is high. Set to NULL to never exit. This parameter was added to prevent setups where some sample points never land on target sites from running indefinitely.

nBurnin

For method set to "MCMC", estTransition runs a JAGS multinomial non-Markovian transitions model, for which it needs the number of burn-in samples before beginning to store results. Default 5000.

nChains

For method set to "MCMC", estTransition runs a JAGS multinomial non-Markovian transitions model, for which it needs the number of MCMC chains (to test for convergence). Default 3.

nThin

For method set to "MCMC", estTransition runs a JAGS multinomial non-Markovian transitions model, for which it needs the thinning rate. Default 1.

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 banded and reencountered are assumed to have no overlap with other data types, so none of this applies to those.

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 modelCountDataJAGS) or matrix with at least nSamples rows. If there are more than [number target sites] columns, the relevant columns should be labeled "relN[1]" through "relN[number target sites]".

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 originAssignment and targetAssignment, not originPoints and targetPoints). However, MCMC is usually faster (and may be a bit more accurate). The third option, "m-out-of-n-bootstrap", is still under development and should be left alone.

m

We read that the m-out-of-n-bootstrap 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 m-out-of-n-bootstrap, m is the number of samples taken each time (less than the true sample size, n). If the "m-out-of-n-bootstrap" is chosen under method but this is left blank, currently the default is n/4, rounded up (no idea if that is reasonable).

psiPrior

matrix with same dimensions as psi. Only relevant when method is "MCMC". Each row provides a Dirichlet (https://en.wikipedia.org/wiki/Dirichlet_distribution) prior on the transition probabilities from that origin site. The default (NULL) supplies Dirichlet parameters of all 1s, which is a standard uninformative Dirichlet prior. Setting these to other positive numbers is useful when you think a priori that certain transitions are unlikely, but don't want to rule them out altogether using fixedZero.

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).

Value

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 Bias-corrected 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".

See Also

estStrength, plot.estMigConnectivity, estMC, estMantel

Examples


  ##############################################################################
  # 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, # Light-level 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, # Non-breeding/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/SMBC-NZP/MigConnectivity/blob/master/data-raw/',
                   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/SMBC-NZP/MigConnectivity/blob/master/data-raw/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


SMBC-NZP/MigConnectivity documentation built on March 26, 2024, 4:22 p.m.