inst/examples/estNMCExamples.R

In MigConnectivity: Estimate Migratory Connectivity for Migratory Animals

\donttest{
  set.seed(101)
  # Uncertainty in detection (RMark estimates)
  # Number of resampling iterations for generating confidence intervals
  nSamplesCMR <- 100
  nSimulationsCMR <- 10

  # the second intermediate psi scenario, the "low" level
  psiTrue <- samplePsis[["Low"]]
  originRelAbundTrue <- rep(0.25, 4)
  trueNMC <- calcNMC(psiTrue, originRelAbund = originRelAbundTrue)
  trueNMC

  # Storage matrix for samples
  cmrNMCSample <- matrix(NA, nSamplesCMR, nSimulationsCMR)
  summaryCMR <- data.frame(Simulation = 1:nSimulationsCMR, True=trueNMC$NMC,
                           mean=NA, se=NA, lcl=NA, ucl=NA)
  # Get 'RMark' psi estimates and estimate MC from each
  for (r in 1:nSimulationsCMR) {
    cat("Simulation",r,"of",nSimulationsCMR,"\n")
    # Note: getCMRexample() requires a valid internet connection and that GitHub
    # is accessible
    fm <- getCMRexample(r)
    results <- estNMC(psi = fm,
                      originSites = 5:8, targetSites = c(3,2,1,4),
                      nSamples = nSamplesCMR, verbose = 0)
    cmrNMCSample[ , r] <- results$NMC$sample
    summaryCMR$mean[r] <- results$NMC$mean
    summaryCMR$se[r] <- results$NMC$se
    # Calculate confidence intervals using quantiles of sampled MC
    summaryCMR[r, c('lcl', 'ucl')] <- results$NMC$simpleCI
  }

  summaryCMR <- transform(summaryCMR, coverage = (True>=lcl & True<=ucl))
  summaryCMR
  summary(summaryCMR)
  biasCMR <- mean(summaryCMR$mean) - trueNMC$NMC
  biasCMR
  mseCMR <- mean((summaryCMR$mean - trueNMC$NMC)^2)
  mseCMR
  rmseCMR <- sqrt(mseCMR)
  rmseCMR

  # Simulation of BBS data to quantify uncertainty in relative abundance
  nSamplesAbund <- 700 #1700 are stored
  nSimulationsAbund <- 10
  #\dontrun{
  #  nSamplesAbund <- 1700
  #}
  # Storage matrix for samples
  abundNMCaSample <- matrix(NA, nSamplesAbund, nSimulationsAbund)
  summaryAbund <- data.frame(Simulation = 1:nSimulationsAbund,
                             True = trueNMC$NMCa,
                             mean = NA, se = NA, lcl = NA, ucl = NA)
  for (r in 1:nSimulationsAbund) {
    cat("Simulation",r,"of",nSimulationsAbund,"\n")
    row0 <- nrow(abundExamples[[r]]) - nSamplesAbund
    results <- estNMC(originRelAbund = abundExamples[[r]], psi = psiTrue,
                           row0 = row0, nSamples = nSamplesAbund, verbose = 2)
    abundNMCaSample[ , r] <- results$NMCa$sample
    summaryAbund$mean[r] <- results$NMCa$mean
    summaryAbund$se[r] <- results$NMCa$se
    # Calculate confidence intervals using quantiles of sampled MC
    summaryAbund[r, c('lcl', 'ucl')] <- results$NMCa$simpleCI
  }

  summaryAbund <- transform(summaryAbund, coverage = (True >= lcl & True <= ucl))
  summaryAbund
  summary(summaryAbund)
  biasAbund <- mean(summaryAbund$mean) - trueNMC$NMCa
  biasAbund
  mseAbund <- mean((summaryAbund$mean - trueNMC$NMCa)^2)
  mseAbund
  rmseAbund <- sqrt(mseAbund)
  rmseAbund

  # Ovenbird example with GL and GPS data
  data(OVENdata) # Ovenbird

  nSamplesGLGPS <- 100 # Number of bootstrap iterations, set low for example

  # Estimate transition probabilities
  Combined.psi<-estTransition(isGL=OVENdata$isGL, #Light-level geolocator (T/F)
                              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 = 3,   # output options
                  nSamples = nSamplesGLGPS, # This is set low for example
                  resampleProjection = sf::st_crs(OVENdata$targetPoints),
                  nSim = 1000)

  # Can estimate NMC from previous psi estimate
  Combo.NMC1 <- estNMC(psi = Combined.psi,
                       nSamples = nSamplesGLGPS)
  Combo.NMC1

  # Doesn't have to be an estPsi object - can simply be array of psi samples
  Combo.MC2 <- estNMC(psi = Combined.psi$psi$sample, # Array of samples
                      originNames = Combined.psi$input$originNames,
                      nSamples = nSamplesGLGPS)
  Combo.MC2

  # Can estimate NMC from previous psi estimate and abundance estimate
  Combo.NMC3 <- estNMC(psi = Combined.psi,
                       originRelAbund = OVENdata$originRelAbund,
                       nSamples = nSamplesGLGPS)
  Combo.NMC3
}