inst/doc/olderFunctionality.R

In MigConnectivity: Estimate Migratory Connectivity for Migratory Animals

## ----echo = FALSE-------------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", tidy = FALSE)

## ----eval = FALSE-------------------------------------------------------------
# install.packages('devtools')
# 
# devtools::install_github("SMBC-NZP/MigConnectivity", build_vignettes = TRUE)
# 
# library(MigConnectivity)

## ----warning=FALSE,message=FALSE,echo=FALSE-----------------------------------
library(MigConnectivity)

## -----------------------------------------------------------------------------
mlogitMat <- function(slope, dist) {
   preMat <- exp(-slope/mean(dist)*dist)
   diag(preMat) <- 0
   outMat <- t(apply(preMat, 1, function(x) x/(1 + sum(x))))
   diag(outMat) <- 1 - rowSums(outMat)
   return(outMat)
}


## -----------------------------------------------------------------------------
mlogitMC <- function(slope, 
                     MC.in, 
                     origin.dist,
                     target.dist,
                     origin.abund,
                     sample.size) {
  
    nBreeding <- nrow(origin.dist)
    nWintering <- nrow(target.dist)
    psi <- mlogitMat(slope, origin.dist)
    
  if (any(psi<0))
    
    return(5*slope^2)
    
  MC <- calcMC(origin.dist, 
               target.dist, 
               originRelAbund = origin.abund,
               psi, 
               sampleSize = sample.size)
  
  return((MC.in - MC)^2)
}

## -----------------------------------------------------------------------------
calcStrengthInd <- function(originDist, 
                            targetDist, 
                            locations, 
                            resamp=1000, 
                            verbose = 0) {
nInd <- dim(locations)[1]
originDist2 <- targetDist2 <- matrix(0, nInd, nInd)
for (i in 1:(nInd-1)) {
 for (j in (i+1):nInd) {
      originDist2[i,j] <- originDist[locations[i,1,1,1], locations[j,1,1,1]]
      targetDist2[i,j] <- targetDist[locations[i,2,1,1], locations[j,2,1,1]]
      originDist2[j,i] <- originDist[locations[i,1,1,1], locations[j,1,1,1]]
      targetDist2[j,i] <- targetDist[locations[i,2,1,1], locations[j,2,1,1]]
    }
  }
return(ncf::mantel.test(originDist2, targetDist2, resamp=resamp, quiet = !verbose))
} 

## -----------------------------------------------------------------------------
calcPsiMC <- function(originDist, 
                      targetDist, 
                      originRelAbund, 
                      locations, 
                      years = 1, 
                      months = 1, 
                      verbose=FALSE) {
  
  nOrigin <- nrow(originDist)
  nTarget <- nrow(targetDist)
  psiMat <- matrix(0, nOrigin, nTarget)
  nInd <- dim(locations)[1]
  nYears <- dim(locations)[3]
  nMonths <- dim(locations)[4]
for (i in 1:nInd) {
  if (i %% 1000 == 0 && verbose) #
      cat("Individual", i, "of", nInd, "\n")
    originMat <- locations[i, 1, years, months]
    targetMat <- locations[i, 2, years, months]
    bIndices <- which(!is.na(originMat))
    wIndices <- which(!is.na(targetMat))
  if (length(bIndices) && length(wIndices))
    for (bi in bIndices)
      for (wi in wIndices)
       psiMat[originMat[bi], targetMat[wi]] <- psiMat[originMat[bi], 
                                                      targetMat[wi]] + 1
  }
  psiMat <- apply(psiMat, 2, "/", rowSums(psiMat))
  MC <- calcMC(originDist, targetDist, psi = psiMat,
               originRelAbund = originRelAbund, sampleSize = nInd)
  return(list(psi=psiMat, MC=MC))
}

## -----------------------------------------------------------------------------
changeLocations <- function(animalLoc, 
                            breedingSiteTrans, 
                            winteringSiteTrans) {
  animalLoc[,1,,] <- breedingSiteTrans[animalLoc[,1,,]]
  animalLoc[,2,,] <- winteringSiteTrans[animalLoc[,2,,]]
  return(animalLoc)
}

## -----------------------------------------------------------------------------
simLocationError <- function(targetPoints, 
                             targetSites, 
              							 geoBias, 
              							 geoVCov, 
                             projection, 
              							 verbose = 0, 
              							 nSim = 1000) {
  
if(!inherits(targetPoints,"sf")){
  targetPoints <- sf::st_as_sf(targetPoints)}

nAnimals <- nrow(targetPoints)

geoBias2 <- matrix(rep(geoBias, nSim), nrow=nSim, byrow=TRUE)
target.sample <- rep(NA, nAnimals)
target.point.sample <- matrix(NA, nAnimals, 2)

for(i in 1:nAnimals){
  if (verbose > 0)
    cat('Animal', i, 'of', nAnimals)
  draws <- 0
  while (is.na(target.sample[i])) {
    draws <- draws + 1
    # Sample random point for each bird 
    # from parametric distribution of NB error
    ps <- MASS::mvrnorm(n=nSim, 
                  mu=cbind(sf::st_coordinates(targetPoints)[i,1],
                           sf::st_coordinates(targetPoints)[i,2]),
                  Sigma=geoVCov)+
                  geoBias2
    
    colnames(ps) <- c("Longitude","Latitude")
    
    point.sample <- sf::st_as_sf(
                  data.frame(ps),
                  coords = c("Longitude","Latitude"),
                  crs = sf::st_crs(targetPoints))
	  
    # filtered to stay in NB areas (land)
    
    target.sample0 <- unlist(unclass(sf::st_intersects(point.sample,
                                                       targetSites)))
    target.sample[i]<-target.sample0[!is.na(target.sample0)][1]
  }
	
  target.point.sample[i, ]<-sf::st_coordinates(point.sample[!is.na(target.sample0),])[1,]
  if (verbose > 0)
    cat(' ', draws, 'draws\n')
}
  return(list(targetSample = target.sample, 
              targetPointSample = target.point.sample))
}

## -----------------------------------------------------------------------------
toPoly <- function(siteCentroids,
                   projection.in = 4326,
                   projection.out = NA,
                   resolution = NA,
                   order.by.input = TRUE)
{
# This automatically sets the resolution so that all polygons touch and
# cover the entire surface. Alternatively the user can supply the 
# resolution of the cells in the units of the input projection (projection.in) 
  colnames(siteCentroids) <- c("Longitude","Latitude")
  
  if(is.na(resolution)){
     long <- unique(siteCentroids[,1])
     lat <- unique(siteCentroids[,2])
     if (length(long)>1)
       long.res <- abs(long[2]-long[1])
     else
       long.res <- 1
     if (length(lat)>1)
       lat.res <- abs(lat[2]-lat[1])
     else
       lat.res <- 1
     resolution <- c(long.res,lat.res)
  }
  centers <- sf::st_as_sf(data.frame(siteCentroids),
                          coords = c("Longitude","Latitude"),
                          crs = projection.in)
  polys <- sf::st_as_sf(sf::st_make_grid(centers,
                                 crs = projection.in, 
                                 cellsize = resolution, 
                                 offset = apply(siteCentroids, 2, min) - 
                                   resolution/2))
  

  if(order.by.input){
    # Reorders the polygons to match that of the input #
    orders <- suppressMessages(unlist(unclass(sf::st_intersects(centers,
                                                         polys,
                                                         sparse = TRUE))))
    polys <- polys[orders,] 
  }
  
  if(!is.na(projection.out)){
    polys <- sf::st_transform(polys, projection.out)
  }
  
  return(polys)
}

## -----------------------------------------------------------------------------
nBreeding <- 3 #number of breeding regions
nNonBreeding <- 3 #number of non-breeding regions

## -----------------------------------------------------------------------------
#distances between centroids of three breeding regions
breedDist <- matrix(c(0, 1, 2,
                      1, 0, 1,
                      2, 1, 0), nBreeding, nBreeding) 

#distances between centroids of three non-breeding regions
nonBreedDist <- matrix(c(0, 5, 10,
                         5, 0, 5,
                        10, 5, 0), nNonBreeding, nNonBreeding)

## ----echo=FALSE---------------------------------------------------------------
op <- graphics::par(no.readonly = TRUE)
on.exit(graphics::par(op))
par(bty="n")
plot(c(4,5,6),rep(1,3),
     ylim=c(-0.05,1.05),
     xlim=c(-1,11),
     pch=19,
     axes=FALSE,
     yaxt="n",
     xaxt="n",
     cex=2,
     ylab="",
     xlab="Non-breeding",
     main="Breeding")
points(c(0,5,10),rep(0,3),pch=15,cex=2)

## -----------------------------------------------------------------------------
#transition probabilities form each breeding to each non-breeding region
psi <- matrix(c(0.4, 0.35, 0.25,
                0.3, 0.4, 0.3,
                0.2, 0.3, 0.5), nBreeding, nNonBreeding, byrow = TRUE) 


## ----echo=FALSE---------------------------------------------------------------
plot(c(4,5,6),rep(1,3),
     ylim=c(-0.05,1.05),
     xlim=c(-1,11),
     pch=19,
     axes=FALSE,
     yaxt="n",
     xaxt="n",
     cex=2,
     ylab="",
     xlab="Non-breeding",
     main="Breeding")
points(c(0,5,10),rep(0,3),pch=15,cex=2)

segments(4,1,0,0,lwd=0.4*10)
segments(4,1,5,0,lwd=0.35*10)
segments(4,1,10,0,lwd=0.25*10)

segments(5,1,0,0,lwd=0.3*10, lty = 2)
segments(5,1,5,0,lwd=0.4*10, lty = 2)
segments(5,1,10,0,lwd=0.3*10, lty = 2)

segments(6,1,0,0,lwd=0.2*10, lty = 3)
segments(6,1,5,0,lwd=0.3*10, lty = 3)
segments(6,1,10,0,lwd=0.5*10, lty = 3)

## -----------------------------------------------------------------------------
#equal relative abundance among the three breeding regions, must sum to 1                      
relN <- rep(1/nBreeding, nBreeding) 

relN

# Define total sample size for psi data 
# for small sample corrected version of MC 

n <- 250


## -----------------------------------------------------------------------------
# Calculate the strength of migratory connectivity 
MC <- calcMC(originDist = breedDist,
             targetDist = nonBreedDist,
             originRelAbund = relN,
             psi = psi)

round(MC, 3)

MC <- calcMC(originDist = breedDist,
             targetDist = nonBreedDist,
             originRelAbund = relN, 
             psi = psi,
             sampleSize = n)

round(MC, 3)

## ----eval = FALSE-------------------------------------------------------------
# # Load in projections
# data("projections")
# 
# # Define deployment locations (winter) #
# captureLocations<-matrix(c(-77.93,18.04,   # Jamaica
#                            -80.94,25.13,   # Florida
#                            -66.86,17.97),  # Puerto Rico
#                             nrow=3, ncol=2, byrow = TRUE)
# 
# colnames(captureLocations) <- c("Longitude","Latitude")
# 
# # Convert capture locations into points #
# 
# CapLocs <- sf::st_as_sf(data.frame(captureLocations),
#                       coords = c("Longitude","Latitude"),
#                       crs = 4326)
# 
# # Project Capture locations #
# 
# CapLocsM<-sf::st_transform(CapLocs, 'ESRI:54027')
# 
# # Retrieve raw non-breeding locations from github
# # First grab the identity of the bird so we can loop through the files
# # For this example we are only interested in the error
# # around non-breeding locations
# # here we grab only the birds captured during the non-breeding season
# # Using paste0 for vignette formatting purposes
# 
# winterBirds <- dget(
# "https://raw.githubusercontent.com/SMBC-NZP/MigConnectivity/master/data-raw/GL_NonBreedingFiles/winterBirds.txt"
#                     )
# 
# # create empty list to store the location data #
# Non_breeding_files <- vector('list',length(winterBirds))
# 
# # Get raw location data from Github #
# for(i in 1:length(winterBirds)){
# Non_breeding_files[[i]] <- dget(paste0("https://raw.githubusercontent.com/",
#                                         "SMBC-NZP/MigConnectivity/master/data-raw/",
#                                         "GL_NonBreedingFiles/NonBreeding_",
#                                         winterBirds[i],".txt"))
# }
# 
# # Remove locations around spring Equinox and potential migration points
# # same NB time frame as Hallworth et al. 2015
# 
# # two steps because subset on shapefile doesn't like it in a single step
# 
# Non_breeding_files <- lapply(Non_breeding_files,
#                       FUN = function(x){
#                       month <- as.numeric(format(x$Date,format = "%m"))
#                                x[which(month != 3 & month != 4),]})
# 
# 
# Jam <- c(1:9)   # locations w/in list of winterBirds captured in Jamaica
# Fla <- c(10:12) # locations w/in list of winterBirds in Florida
# PR <- c(13:16)  # locations w/in list of winterBirds in Puerto Rico
# 
# # Turn the locations into shapefiles #
# 
# NB_GL <- lapply(Non_breeding_files,
#                 FUN = function(x){
#                   sf::st_as_sf(data.frame(x),
#                                coords = c("Longitude",
#                                           "Latitude"),
#                                crs = 4326)})
# 
# # Project into UTM projection #
# 
# NB_GLmeters <- lapply(NB_GL,
#                       FUN = function(x){sf::st_transform(x,'ESRI:54027')})
# 
# # Process to determine geolocator bias and variance-covariance in meters #
# 
# # generate empty vector to store data #
# LongError<-rep(NA,length(winterBirds))
# LatError<-rep(NA,length(winterBirds))
# 
# # Calculate the error in longitude derived
# # from geolocators from the true capture location
# 
# LongError[Jam] <- unlist(lapply(NB_GLmeters[Jam],
#                          FUN = function(x){mean(sf::st_coordinates(x)[,1]-
#                                                   sf::st_coordinates(CapLocsM)[1,1])}))
# 
# LongError[Fla] <- unlist(lapply(NB_GLmeters[Fla],
#                          FUN = function(x){mean(sf::st_coordinates(x)[,1]-
#                                                   sf::st_coordinates(CapLocsM)[2,1])}))
# 
# LongError[PR] <- unlist(lapply(NB_GLmeters[PR],
#                         FUN = function(x){mean(sf::st_coordinates(x)[,1]-
#                                                  sf::st_coordinates(CapLocsM)[3,1])}))
# 
# # Calculate the error in latitude derived from
# # geolocators from the true capture location
# 
# LatError[Jam] <- unlist(lapply(NB_GLmeters[Jam],
#                         FUN = function(x){mean(sf::st_coordinates(x)[,2]-
#                                                   sf::st_coordinates(CapLocsM)[1,2])}))
# 
# LatError[Fla] <- unlist(lapply(NB_GLmeters[Fla],
#                         FUN = function(x){mean(sf::st_coordinates(x)[,2]-
#                                                  sf::st_coordinates(CapLocsM)[2,2])}))
# 
# LatError[PR] <- unlist(lapply(NB_GLmeters[PR],
#                         FUN = function(x){mean(sf::st_coordinates(x)[,2]-
#                                                   sf::st_coordinates(CapLocsM)[3,2])}))
# 
# # Get co-variance matrix for error of
# # known non-breeding deployment sites
# 
# # lm does multivariate normal models if you give it a matrix dependent variable!
# 
# geo.error.model <- lm(cbind(LongError,LatError) ~ 1)
# 
# geo.bias <- coef(geo.error.model)
# geo.vcov <- vcov(geo.error.model)

## -----------------------------------------------------------------------------
data(OVENdata) # Ovenbird 

names(OVENdata)

## ----eval = FALSE-------------------------------------------------------------
# install.packages(c('shape', 'terra', 'maps'))

## ----message = FALSE, warning = FALSE, error=FALSE----------------------------
library(shape)
library(maps)
library(terra)

# Set the crs to WGS84
originSitesWGS84 <- sf::st_transform(OVENdata$originSites, 4326)
targetSitesWGS84 <- sf::st_transform(OVENdata$targetSites, 4326)

# Create a simple plot of the origin and and target sites 
par(mar=c(0,0,0,0))
plot(originSitesWGS84,
     xlim=c(terra::ext(targetSitesWGS84)[1],
            terra::ext(targetSitesWGS84)[2]),

     ylim=c(terra::ext(targetSitesWGS84)[3],
            terra::ext(originSitesWGS84)[4]))

plot(targetSitesWGS84,
     add = TRUE,
     col=c("gray70","gray35","gray10"))

map("world", add=TRUE)

## ----message=FALSE, warning = FALSE, error=FALSE------------------------------
M<-estMC(isGL=OVENdata$isGL, # Logical vector: light-level geolocator(T)/GPS(F)
         geoBias = OVENdata$geo.bias, # Light-level geolocator location bias
         geoVCov = OVENdata$geo.vcov, #Light-level geolocator covariance matrix
         targetDist = OVENdata$targetDist, # Target location distance matrix
         originDist = OVENdata$originDist, # Origin location distance matrix
         targetSites = OVENdata$targetSites, # Non-breeding / target sites
         originSites = OVENdata$originSites, # Breeding / origin sites 
         targetNames = OVENdata$targetNames, # Names of nonbreeding/target sites
         originNames = OVENdata$originNames, # Names of breeding/origin sites
         originPoints = OVENdata$originPoints, # Capture Locations 
         targetPoints = OVENdata$targetPoints, # Target locations from devices
         originRelAbund = OVENdata$originRelAbund, # Origin relative abundances
         resampleProjection = sf::st_crs(OVENdata$targetPoints), 
         verbose = 0,   # output options - see help ??estMC
         nSamples = 10) # This is set low for example 

M

## -----------------------------------------------------------------------------
str(M, max.level = 2)

## ----echo = FALSE, warning = FALSE, error = FALSE, eval = FALSE---------------
# plot(sf::st_geometry(wrld_simple),
#      xlim=sf::st_bbox(OVENdata$targetSites)[c(1,3)],
#      ylim=c(sf::st_bbox(OVENdata$targetSites)[2],
#             sf::st_bbox(OVENdata$originSites)[4]))
# plot(sf::st_geometry(OVENdata$originSites[1,]),add=TRUE,lwd=1.75)
# plot(sf::st_geometry(OVENdata$originSites[2,]),add=TRUE,lwd=1.75)
# plot(sf::st_geometry(OVENdata$targetSites),add=TRUE,lwd=1.5,col=c("gray70","gray35","gray10"))
# 
# legend(400000,-2105956,legend=paste("MC =",round(M$MC$mean,2), "\u00b1", round(M$MC$se,2)),bty="n",cex=1.8,bg="white",xjust=0)
# 
# shape::Arrows(x0 = sf::st_coordinates(sf::st_centroid(OVENdata$originSites[1,]))[,1],
#               y0 = sf::st_coordinates(sf::st_centroid(OVENdata$originSites[1,]))[,2],
#               x1 = sf::st_coordinates(sf::st_centroid(OVENdata$targetSites[2,]))[,1]+80000,
#               y1 = sf::st_bbox(OVENdata$targetSites[2,])[4]+150000,
#               arr.length = 0.3,
#               arr.adj = 0.5,
#               arr.lwd = 1,
#               arr.width = 0.4,
#               arr.type = "triangle",
#               lwd = M$psi$mean[1,2]*10,
#               lty = 1)
# 
# shape::Arrows(sf::st_coordinates(sf::st_centroid(OVENdata$originSites[1,]))[,1],
#              sf::st_coordinates(sf::st_centroid(OVENdata$originSites[1,]))[,2],
#               sf::st_coordinates(sf::st_centroid(OVENdata$targetSites[3,]))[,1],
#               sf::st_bbox(OVENdata$targetSites[3,])[4]+150000,
#               arr.length = 0.3,
#               arr.adj = 0.5,
#               arr.lwd = 1,
#               arr.width = 0.4,
#               arr.type = "triangle",
#               lwd = M$psi$mean[1,3]*10,
#               lty=1)
# 
# shape::Arrows(sf::st_coordinates(sf::st_centroid(OVENdata$originSites[2,]))[,1],
#               sf::st_coordinates(sf::st_centroid(OVENdata$originSites[2,]))[,2],
#               sf::st_coordinates(sf::st_centroid(OVENdata$targetSites[1,]))[,1],
#               sf::st_bbox(OVENdata$targetSites[1,])[4]+150000,
#               arr.length = 0.3,
#               arr.adj = 0.5,
#               arr.lwd = 1,
#               arr.width = 0.4,
#               arr.type = "triangle",
#               lwd = M$psi$mean[2,1]*10,
#               lty=1)
# 
# shape::Arrows(sf::st_coordinates(sf::st_centroid(OVENdata$originSites[2,]))[,1],
#               sf::st_coordinates(sf::st_centroid(OVENdata$originSites[2,]))[,2],
#               sf::st_coordinates(sf::st_centroid(OVENdata$targetSites[2,]))[,1]-80000,
#               sf::st_bbox(OVENdata$targetSites[2,])[4]+150000,
#               arr.length = 0.3,
#               arr.adj = 0.5,
#               arr.lwd = 1,
#               arr.width = 0.4,
#               arr.type = "triangle",
#               lwd = M$psi$mean[2,2]*10)
# 
# box(which="plot")

## -----------------------------------------------------------------------------
set.seed(75)

# Parameters for simulations

nSeasons <- 2 # population with two discrete seasons - breeding / non-breeding
nYears <- 10 # Ten years of data
nMonths <- 4 # Four months within each season

nBreeding <- 100 # Number of populations
nWintering <- 100 # Number of populations 

## -----------------------------------------------------------------------------
breedingPos <- matrix(c(rep(seq(-99,-81,2), each=sqrt(nBreeding)),
                        rep(seq(49,31,-2), sqrt(nBreeding))), nBreeding, 2)
winteringPos <- matrix(c(rep(seq(-79,-61,2), each=sqrt(nWintering)),
                         rep(seq(9,-9,-2), sqrt(nWintering))), nWintering, 2)

# calculate distance between populations
breedDist <- distFromPos(breedingPos, 'ellipsoid') 

nonbreedDist <- distFromPos(winteringPos, 'ellipsoid') 

## -----------------------------------------------------------------------------
# Breeding Abundance
breedingN <- rep(500, nBreeding)
breedingRelN <- breedingN/sum(breedingN)

## ----eval = FALSE-------------------------------------------------------------
# # Set up psi matrix
# o <- optimize(mlogitMC, MC.in = 0.25, origin.dist = breedDist,
#               target.dist = nonbreedDist, origin.abund = breedingRelN,
#               sample.size = sum(breedingN),
#               interval = c(1, 3), tol = .Machine$double.eps^0.5)#
# 
# slope <- o$minimum
# 
# psi <- mlogitMat(slope, breedDist)

## ----warning=FALSE,message=FALSE,echo=FALSE-----------------------------------
slope <- 1.92093332054914
psi <- mlogitMat(slope, breedDist)

## -----------------------------------------------------------------------------
# Baseline strength of migratory connectivity
MC <- calcMC(originDist = breedDist, 
             targetDist = nonbreedDist,
             psi = psi, 
             originRelAbund = breedingRelN,
             sampleSize = sum(breedingN))

# Show Results
MC

## -----------------------------------------------------------------------------
sims <- simMove(breedingAbund = breedingN, # Breeding relative abundance
               breedingDist = breedDist,  # Breeding distance
               winteringDist = nonbreedDist, # Non-breeding distance
               psi = psi, # Transition probabilities
               nYears = nYears, # Number of years
               nMonths = nMonths, # Number of months
               winMoveRate = 0, # winter movement rate
               sumMoveRate = 0, # breeding movement rate
               winDispRate = 0, # winter dispersal rate
               sumDispRate = 0, # summer disperal rate
               natalDispRate = 0, # natal dispersal rate
               breedDispRate = 0, # breeding dispersal rate
               verbose = 0) # verbose output


str(sims)

## -----------------------------------------------------------------------------
nScenarios1 <- length(samplePsis) # samplePsis - comes with MigConnectivity package

# Create vector of length nScenarios1
MC1 <- rep(NA, nScenarios1)

# Loop through the different senarios outlined above #
for (i in 1:nScenarios1) {
  MC1[i] <- calcMC(originDist = sampleOriginDist[[1]],
                 targetDist = sampleTargetDist[[1]],
                 psi = samplePsis[[i]],
                 originRelAbund = sampleOriginRelN[[1]])
}

# Give meaningful names to the MC1 vector
names(MC1) <- names(samplePsis)

# Print results 
round(MC1, 6)

## -----------------------------------------------------------------------------
nScenarios2 <- length(sampleOriginPos)

MC2 <- matrix(NA, nScenarios1, nScenarios2)

rownames(MC2) <- names(samplePsis)

colnames(MC2) <- names(sampleOriginPos)

for (i in 1:nScenarios1) {
  for (j in 1:nScenarios2) {
    MC2[i, j] <- calcMC(originDist = sampleOriginDist[[j]],
                      targetDist = sampleTargetDist[[j]],
                      psi = samplePsis[[i]],
                      originRelAbund = sampleOriginRelN[[1]])
  }
}

# Print results #
t(round(MC2, 4))

## -----------------------------------------------------------------------------

MC.diff2 <- apply(MC2, 2, "-", MC2[ , 1])

t(round(MC.diff2, 4))

## -----------------------------------------------------------------------------
nScenarios3 <- length(sampleOriginRelN)

MC3 <- matrix(NA, nScenarios1, nScenarios3)

rownames(MC3) <- names(samplePsis)
colnames(MC3) <- names(sampleOriginRelN)

for (i in 1:nScenarios1) {
for (j in 1) {
for (k in 1:nScenarios3) {
MC3[i, k] <- calcMC(originDist = sampleOriginDist[[j]],
                    targetDist = sampleTargetDist[[j]], 
                    psi = samplePsis[[i]],
                    originRelAbund = sampleOriginRelN[[k]])
}
}
}

t(round(MC3, 4)) # linear arrangement

## -----------------------------------------------------------------------------
MC4 <- matrix(NA, nScenarios1, nScenarios3)

rownames(MC4) <- names(samplePsis)
colnames(MC4) <- names(sampleOriginRelN)

for (i in 1:nScenarios1) {
for (j in nScenarios2) {
for (k in 1:nScenarios3) {
MC4[i, k] <- calcMC(originDist = sampleOriginDist[[j]],
                    targetDist = sampleTargetDist[[j]],
                    psi = samplePsis[[i]],
                    originRelAbund = sampleOriginRelN[[k]])
}
}
}

t(round(MC4, 4)) # grid arrangement

## -----------------------------------------------------------------------------
MC5 <- matrix(NA, nScenarios1, nScenarios3)

rownames(MC5) <- names(samplePsis)
colnames(MC5) <- names(sampleOriginRelN)

for (i in 1:nScenarios1) {
for (j in 6) {
for (k in 1:nScenarios3) {
MC5[i, k] <- calcMC(originDist = sampleOriginDist[[j]],
                    targetDist = sampleTargetDist[[j]], 
                    psi = samplePsis[[i]], 
                    originRelAbund = sampleOriginRelN[[k]])
}
}
}

t(round(MC5, 4)) # breeding grid, winter linear arrangement

## ----eval=FALSE---------------------------------------------------------------
# #Run functions and parameters above first
# set.seed(75)
# 
# scenarios14 <- c("Base",
#                  "Breeding10",
#                  "Wintering10",
#                  "Breeding10Wintering10",
#                  "Breeding4",
#                  "Wintering4",
#                  "Breeding4Wintering10",
#                  "Breeding10Wintering4",
#                  "Breeding4Wintering4",
#                  "CentroidSampleBreeding10",
#                  "CentroidSampleBreeding10Wintering10",
#                  "CentroidSampleBreeding10Wintering4",
#                  "CentroidSampleBreeding4",
#                  "CentroidSampleBreeding4Wintering10",
#                  "CentroidSampleBreeding4Wintering4")
# 
# # Each element is for a scenario (see above 1-8), transferring from natural
# # breeding populations to defined ones
# breedingSiteTrans14 <- list(1:nBreeding,
#                            rep(1:10, each=10),
#                            1:nBreeding,
#                            rep(1:10, each=10),
#                            c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                            1:nBreeding,
#                            c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                            rep(1:10, each=10),
#                            c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                            rep(1:10, each=10),
#                            rep(1:10, each=10),
#                            rep(1:10, each=10),
#                            c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                            c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                            c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)))
# 
# # Same for non-breeding populations
# winteringSiteTrans14 <- list(1:nWintering,
#                              1:nWintering,
#                              rep(1:10, each=10),
#                              rep(1:10, each=10),
#                              1:nWintering,
#                              c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                              rep(1:10, each=10),
#                              c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                              c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                              1:nWintering,
#                              rep(1:10, each=10),
#                              c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                              1:nWintering,
#                              rep(1:10, each=10),
#                              c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)))
# 
# # Examine the transfers in matrix form
# #lapply(breedingSiteTrans14, matrix, nrow=10, ncol=10)
# #lapply(winteringSiteTrans14, matrix, nrow=10, ncol=10)
# 
# #positions of the human defined populations
# breedingPos14 <- list(breedingPos,
#                       rowsum(breedingPos, rep(1:10, each=10))/10,
#                       breedingPos,
#                       rowsum(breedingPos, rep(1:10, each=10))/10,
#                       rowsum(breedingPos, c(rep(1:2, 5, each=5),
#                                             rep(3:4, 5, each=5)))/25,
#                       breedingPos,
#                       rowsum(breedingPos, c(rep(1:2, 5, each=5),
#                                             rep(3:4, 5, each=5)))/25,
#                       rowsum(breedingPos, rep(1:10, each=10))/10,
#                       rowsum(breedingPos, c(rep(1:2, 5, each=5),
#                                             rep(3:4, 5, each=5)))/25,
#                       rowsum(breedingPos, rep(1:10, each=10))/10,
#                       rowsum(breedingPos, rep(1:10, each=10))/10,
#                       rowsum(breedingPos, rep(1:10, each=10))/10,
#                       rowsum(breedingPos, c(rep(1:2, 5, each=5),
#                                             rep(3:4, 5, each=5)))/25,
#                       rowsum(breedingPos, c(rep(1:2, 5, each=5),
#                                             rep(3:4, 5, each=5)))/25,
#                       rowsum(breedingPos, c(rep(1:2, 5, each=5),
#                                             rep(3:4, 5, each=5)))/25)
# 
# winteringPos14 <- list(winteringPos,
#                        winteringPos,
#                        rowsum(winteringPos, rep(1:10, each=10))/10,
#                        rowsum(winteringPos, rep(1:10, each=10))/10,
#                        winteringPos,
#                        rowsum(winteringPos, c(rep(1:2, 5, each=5),
#                                               rep(3:4, 5, each=5)))/25,
#                        rowsum(winteringPos,rep(1:10, each=10))/10,
#                        rowsum(winteringPos, c(rep(1:2, 5, each=5),
#                                               rep(3:4, 5, each=5)))/25,
#                        rowsum(winteringPos, c(rep(1:2, 5, each=5),
#                                               rep(3:4, 5, each=5)))/25,
#                        winteringPos,
#                        rowsum(winteringPos, rep(1:10, each=10))/10,
#                        rowsum(winteringPos, c(rep(1:2, 5, each=5),
#                                               rep(3:4, 5, each=5)))/25,
#                        winteringPos,
#                        rowsum(winteringPos, rep(1:10, each=10))/10,
#                        rowsum(winteringPos, c(rep(1:2, 5, each=5),
#                                               rep(3:4, 5, each=5)))/25)
# 
# 
# # Calculate distances between defined breeding populations
# breedDist14 <- lapply(breedingPos14, distFromPos, surface = 'ellipsoid')
# 
# # Calculate distances between defined non-breeding populations
# nonbreedDist14 <- lapply(winteringPos14, distFromPos, surface = 'ellipsoid')
# 
# # Numbers of defined populations
# nBreeding14 <- c(100, 10, 100, 10, 4, 100, 4, 10, 4, 10, 10, 10, 4, 4, 4)
# 
# nWintering14 <- c(100, 100, 10, 10, 100, 4, 10, 4, 4, 100, 10, 4, 100, 10, 4)
# 
# # Relative abundance by scenario and breeding population
# breedingRelN14 <- lapply(nBreeding14, function(x) rep(1/x, x))
# 
# MC14 <- 0.25 # Same for all scenarios in this set
# 
# nSample14 <- 1000 # Total number sampled per simulation
# 
# # How many sampled from each natural population (sampling scenarios separate
# # from definition scenarios)
# 
# sampleBreeding14 <- list(round(breedingRelN14[[1]]*nSample14),
#                          c(rep(0, 22), round(breedingRelN14[[5]][1]*nSample14),
#                            rep(0, 4), round(breedingRelN14[[5]][2]*nSample14),
#                            rep(0, 44), round(breedingRelN14[[5]][3]*nSample14),
#                            rep(0, 4), round(breedingRelN14[[5]][4]*nSample14),
#                            rep(0, 22)),
#                          rep(c(rep(0, 4),
#                              rep(round(breedingRelN14[[2]][1]*nSample14/2), 2),
#                              rep(0, 4)), 10))
# 
# # for the baseline use the simulation from above, sims
# 
# animalLoc14base <- sims$animalLoc
# 
# 
# # Number of scenarios and number of simulations to run
# nScenarios14 <- length(breedingSiteTrans14)
# nSims14 <- 100
# 
# # Connections between scenarios and sampling regimes
# scenarioToSampleMap14 <- c(rep(1, 9), rep(3, 3), rep(2, 3))
# 
# # Set up data structures for storing results
# animalLoc14 <- vector("list", nScenarios14) #making an empty list to fill
# 
# compare14 <- data.frame(Scenario = c("True", scenarios14),
#                         MC = c(MC14, rep(NA, nScenarios14)),
#                         Mantel = c(MC14, rep(NA, nScenarios14)))
# 
# compare14.array <- array(NA, c(nSims14, nScenarios14, 2), dimnames =
#                            list(1:nSims14,
#                            scenarios14,
#                            c("MC", "Mantel")))
# 
# results14 <- vector("list", nScenarios14)
# 
# # Run simulations
# for (sim in 1:nSims14) {
#   cat("Simulation", sim, "of", nSims14, '\n')
#   sim14 <- lapply(sampleBreeding14,
#                   simMove,
#                   breedingDist = breedDist14[[1]],
#                   winteringDist=nonbreedDist14[[1]],
#                   psi=psi,
#                   nYears=nYears,
#                   nMonths=nMonths)
# 
#   for (i in 1:nScenarios14) {
#    cat("\tScenario", i, "\n")
#    animalLoc14[[i]]<-changeLocations(
#      sim14[[scenarioToSampleMap14[i]]]$animalLoc,
#                                      breedingSiteTrans14[[i]],
#                                      winteringSiteTrans14[[i]])
# 
#    results14[[i]] <- calcPsiMC(breedDist14[[i]], nonbreedDist14[[i]],
#                                animalLoc14[[i]],
#                                originRelAbund = breedingRelN14[[i]],
#                                verbose = FALSE)
# 
#    compare14.array[sim, i, 'MC'] <- results14[[i]]$MC
# 
#    compare14.array[sim,i,'Mantel']<-calcStrengthInd(breedDist14[[1]],
#                                                     nonbreedDist14[[1]],
#                                    sim14[[scenarioToSampleMap14[i]]]$animalLoc,
#                                                     resamp=0)$correlation
#   }
# }
# 
# # Compute means for each scenario
# compare14$MC[1:nScenarios14 + 1] <- apply(compare14.array[,,'MC'], 2, mean)
# compare14$Mantel[1:nScenarios14 + 1] <- apply(compare14.array[,,'Mantel'], 2,
#                                               mean)
# 
# compare14 <- transform(compare14, MC.diff=MC - MC[1],
#                        Mantel.diff=Mantel - Mantel[1])
# compare14
# 

## ----eval = FALSE-------------------------------------------------------------
# set.seed(75)
# 
# # Transfer between true populations and researcher defined ones (only for
# # breeding, as not messing with winter populations here)
# 
# breedingSiteTrans15 <- list(1:100,
#                             c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                             c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                             1:100,
#                             c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                             c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)))
# 
# nScenarios15 <- length(breedingSiteTrans15)
# 
# nSims15 <- 100
# 
# # Basing positions of researcher defined breeding populations on above
# breedingPos15 <- list(breedingPos,
#                       breedingPos14[[5]],
#                       breedingPos14[[5]],
#                       breedingPos,
#                       breedingPos14[[5]],
#                       breedingPos14[[5]])
# 
# winteringPos15 <- rep(list(winteringPos), nScenarios15)
# 
# breedDist15 <- lapply(breedingPos15, distFromPos)
# 
# nonbreedDist15 <- lapply(winteringPos15, distFromPos)
# 
# nBreeding15 <- rep(c(100, 4, 4), 2)
# 
# nWintering15 <- rep(100, nScenarios15)
# 
# # Highest abundance in lower right corner, lowest in top left
# # Making symmetrical
# 
# breedingN15base <- rep(NA, 100)
# for (i in 1:10) #row
#   for (j in 1:10)  #column
#     breedingN15base[i+10*(j-1)] <- 500 + 850*i*j
# 
# sum(breedingN15base)
# 
# # For researcher defined populations
# breedingN15 <- lapply(breedingSiteTrans15, rowsum, x=breedingN15base)
# 
# breedingRelN15 <- lapply(breedingN15, "/", sum(breedingN15base))
# 
# nSample15 <- 1000 # Total number sampled per simulation
# 
# # Number sampled per natural population
# sampleBreeding15 <- list(round(breedingRelN15[[1]]*nSample15),
#                          c(rep(0, 22), round(breedingRelN15[[3]][1]*nSample15),
#                            rep(0, 4), round(breedingRelN15[[3]][2]*nSample15),
#                            rep(0, 44), round(breedingRelN15[[3]][3]*nSample15),
#                            rep(0, 4), round(breedingRelN15[[3]][4]*nSample15),
#                            rep(0, 22)),
#                          round(breedingRelN15[[1]]*nSample15)[100:1],
#                          c(rep(0, 22), round(breedingRelN15[[3]][1]*nSample15),
#                            rep(0, 4), round(breedingRelN15[[3]][2]*nSample15),
#                            rep(0, 44), round(breedingRelN15[[3]][3]*nSample15),
#                            rep(0, 4), round(breedingRelN15[[3]][4]*nSample15),
#                            rep(0, 22))[100:1])
# 
# # Set up psi matrix
# o15 <- optimize(mlogitMC, MC.in = 0.25, origin.dist = breedDist15[[1]],
#                 target.dist = nonbreedDist15[[1]],
#                 origin.abund = breedingN15[[1]]/sum(breedingN15[[1]]),
#                 sample.size = sum(breedingN15[[1]]),
#                 interval = c(0,10),
#                 tol = .Machine$double.eps^0.5)
# 
# slope15 <- o15$minimum
# 
# psi15 <- mlogitMat(slope15, breedDist15[[1]])
# 
# # Baseline strength of migratory connectivity
# MC15 <- calcMC(originDist = breedDist15[[1]],
#                targetDist = nonbreedDist15[[1]],
#                psi = psi15,
#                originRelAbund = breedingN15[[1]]/sum(breedingN15[[1]]),
#                sampleSize = sum(breedingN15[[1]]))
# 
# # Run sampling regimes
# scenarioToSampleMap15 <- c(1, 1, 2, 3, 3, 4)
# 
# animalLoc15 <- vector("list", nScenarios15)
# 
# results15 <- vector("list", nScenarios15)
# 
# compare15 <- data.frame(Scenario = c("True",
#                                      "Base",
#                                      "Breeding4",
#                                      "CentroidSampleBreeding4",
#                                      "BiasedSample",
#                                      "BiasedSampleBreeding4",
#                                      "BiasedCentroidSampleBreeding4"),
#                         MC = c(MC15, rep(NA, nScenarios15)),
#                         Mantel = c(MC15, rep(NA, nScenarios15)))
# 
# compare15.array <- array(NA, c(nSims15, nScenarios15, 2),
#                          dimnames = list(1:nSims15,
#                                     c("Base", "Breeding4",
#                                       "CentroidSampleBreeding4",
#                                       "BiasedSample", "BiasedSampleBreeding4",
#                                       "BiasedCentroidSampleBreeding4"),
#                                     c("MC", "Mantel")))
# 
# 
# for (sim in 1:nSims15) {
#   cat("Simulation", sim, "of", nSims15, '\n')
# 
#   sim15 <- lapply(sampleBreeding15,
#                   simMove,
#                   breedingDist = breedDist15[[1]],
#                   winteringDist=nonbreedDist15[[1]],
#                   psi=psi15,
#                   nYears=nYears,
#                   nMonths=nMonths)
# 
# for (i in 1:nScenarios15) {
#   cat("\tScenario", i, "\n")
#     animalLoc15[[i]] <- changeLocations(
#                         animalLoc=sim15[[scenarioToSampleMap15[i]]]$animalLoc,
#                         breedingSiteTrans = breedingSiteTrans15[[i]],
#                         winteringSiteTrans = 1:nWintering15[i])
# 
#     results15[[i]] <- calcPsiMC(originDist = breedDist15[[i]],
#                                 targetDist = nonbreedDist15[[i]],
#                                 originRelAbund = breedingRelN15[[i]],
#                                 locations = animalLoc15[[i]],
#                                 verbose = FALSE)
# 
# compare15.array[sim, i, 'MC'] <- results15[[i]]$MC
# 
# compare15.array[sim, i, 'Mantel'] <- calcStrengthInd(breedDist15[[1]],
#                                                      nonbreedDist15[[1]],
#                                    sim15[[scenarioToSampleMap15[i]]]$animalLoc,
#                                                     resamp=0)$correlation
#   }
# }
# 
# compare15$MC[1:nScenarios15+1] <- apply(compare15.array[,,'MC'], 2, mean,
#                                         na.rm = TRUE)
# 
# compare15$Mantel[1:nScenarios15+1] <- apply(compare15.array[,,'Mantel'], 2,
#                                             mean, na.rm = TRUE)
# compare15 <- transform(compare15, MC.diff=MC - MC[1],
#                        Mantel.diff=Mantel - Mantel[1],
#                        MC.prop=MC/MC[1],
#                        Mantel.prop=Mantel/Mantel[1])
# 
# 
# compare15a <- as.matrix(compare15[1 + 1:nScenarios15,
#                                   c("MC", "MC.diff", "Mantel", "Mantel.diff")])
# 
# rownames(compare15a) <- compare15$Scenario[1 + 1:nScenarios15]

## ----eval = FALSE-------------------------------------------------------------
# set.seed(75)
# 
# # Transfer between true populations and researcher defined ones
# # (only for breeding, as not messing with winter populations here)
# breedingSiteTrans16 <- list(1:100, c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                             c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)), 1:100,
#                             c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)),
#                             c(rep(1:2, 5, each=5), rep(3:4, 5, each=5)))
# 
# #lapply(breedingSiteTrans16, matrix, nrow=10, ncol=10)
# 
# nScenarios16 <- length(breedingSiteTrans16)
# 
# nSims16 <- 100
# 
# # Basing positions of researcher defined breeding populations on above
# breedingPos16 <- breedingPos15
# winteringPos16 <- winteringPos15
# 
# breedDist16 <- breedDist15
# nonbreedDist16 <- nonbreedDist15
# nBreeding16 <- nBreeding15
# nWintering16 <- rep(100, nScenarios16)
# 
# # Highest abundance in lower right corner, lowest in top left
# # In fact basing on distance from top left population
# 
# breedingN16base <- breedingN15base
# breedingN16 <- breedingN15
# breedingRelN16 <- lapply(breedingN16, "/", sum(breedingN16base))
# 
# # Set up psi matrix
# # Each quadrant of breeding range has different MC
# 
# MC.levels16 <- seq(0.15, 0.6, 0.15)
# 
# nLevels16 <- 4
# 
# psi16 <- matrix(NA, nBreeding16[1], nWintering16[1])
# 
# 
# for (i in 1:nLevels16) {
#   cat("MC", MC.levels16[i])
#   # Find a psi matrix that produces the given MC (for whole species)
#   o16a <- optimize(mlogitMC, MC.in = MC.levels16[i],
#                    origin.dist = breedDist16[[1]],
#                    target.dist = nonbreedDist16[[1]],
#                    origin.abund = breedingN16[[1]]/sum(breedingN16[[1]]),
#                    sample.size = sum(breedingN16[[1]]),
#                    interval=c(0,10), tol=.Machine$double.eps^0.5)
# 
#   slope16a <- o16a$minimum
# 
#   cat(" slope", slope16a, "\n")
# 
#   psi16a <- mlogitMat(slope16a, breedDist16[[1]])
# 
#   # Then use the rows of that psi matrix only for the one breeding quadrant
#   rows <- 50*(i %/% 3) + rep(1:5, 5) + rep(seq(0, 40, 10), each=5) +
#     ((i-1) %% 2) * 5
#   psi16[rows, ] <- psi16a[rows, ]
# }
# 
# 
# # Baseline strength of migratory connectivity
# MC16 <- calcMC(originDist = breedDist16[[1]],
#                targetDist = nonbreedDist16[[1]],
#                psi = psi16,
#                originRelAbund = breedingN16[[1]]/sum(breedingN16[[1]]),
#                sampleSize = sum(breedingN16[[1]]))
# 
# 
# # Set up sampling regimes (different number than number of scenarios)
# nSample16 <- 1000
# sampleBreeding16 <- list(round(breedingRelN16[[1]]*nSample16),
#                          c(rep(0, 22), round(breedingRelN16[[3]][1]*nSample16),
#                            rep(0, 4), round(breedingRelN16[[3]][2]*nSample16),
#                            rep(0, 44), round(breedingRelN16[[3]][3]*nSample16),
#                            rep(0, 4), round(breedingRelN16[[3]][4]*nSample16),
#                            rep(0, 22)),
#                          round(breedingRelN16[[1]]*nSample16)[100:1],
#                          c(rep(0, 22), round(breedingRelN16[[3]][1]*nSample16),
#                            rep(0, 4), round(breedingRelN16[[3]][2]*nSample16),
#                            rep(0, 44), round(breedingRelN16[[3]][3]*nSample16),
#                            rep(0, 4), round(breedingRelN16[[3]][4]*nSample16),
#                            rep(0, 22))[100:1])
# 
# sapply(sampleBreeding16, sum)
# 
# # Run sampling regimes
# scenarioToSampleMap16 <- c(1, 1, 2, 3, 3, 4)
# animalLoc16 <- vector("list", nScenarios16)
# results16 <- vector("list", nScenarios16)
# compare16 <- data.frame(Scenario = c("True",
#                                      "Base",
#                                      "Breeding4",
#                                      "CentroidSampleBreeding4",
#                                      "BiasedSample",
#                                      "BiasedSampleBreeding4",
#                                      "BiasedCentroidSampleBreeding4"),
#                         MC = c(MC16, rep(NA, nScenarios16)),
#                         Mantel = c(MC16, rep(NA, nScenarios16)))
# 
# compare16.array <- array(NA, c(nSims16, nScenarios16, 2),
#                          dimnames = list(1:nSims16,
#                                          c("Base", "Breeding4",
#                                            "CentroidSampleBreeding4",
#                                            "BiasedSample",
#                                            "BiasedSampleBreeding4",
#                                            "BiasedCentroidSampleBreeding4"),
#                                          c("MC", "Mantel")))
# 
# set.seed(80)
# for (sim in 1:nSims16) {
#   cat("Simulation", sim, "of", nSims16, '\n')
# 
#   sim16 <- lapply(sampleBreeding16, simMove, breedingDist = breedDist16[[1]],
#                   winteringDist=nonbreedDist16[[1]], psi=psi16, nYears=nYears,
#                   nMonths=nMonths)
# 
# for (i in 1:nScenarios16) {
#   cat("\tScenario", i, "\n")
#   animalLoc16[[i]]<-changeLocations(sim16[[scenarioToSampleMap16[i]]]$animalLoc,
#                                       breedingSiteTrans16[[i]],
#                                       1:nWintering16[i])
# 
#   results16[[i]] <- calcPsiMC(originDist = breedDist16[[i]],
#                               targetDist = nonbreedDist16[[i]],
#                               locations = animalLoc16[[i]],
#                               originRelAbund = breedingRelN16[[i]],
#                               verbose = FALSE)
# 
# compare16.array[sim, i, 'MC'] <- results16[[i]]$MC
# 
# compare16.array[sim, i, 'Mantel'] <- calcStrengthInd(breedDist16[[1]],
#                                                          nonbreedDist16[[1]],
#                                    sim16[[scenarioToSampleMap16[i]]]$animalLoc,
#                                                          resamp=0)$correlation
#   }
# }
# 
# compare16$MC[1:nScenarios16+1] <- apply(compare16.array[,,'MC'], 2, mean, na.rm = TRUE)
# 
# compare16$Mantel[1:nScenarios16+1] <- apply(compare16.array[,,'Mantel'], 2,
#                                             mean, na.rm = TRUE)
# 
# compare16 <- transform(compare16, MC.diff=MC - MC[1],
#                        Mantel.diff=Mantel - Mantel[1])
# 
# compare16a <- as.matrix(compare16[1 + 1:nScenarios16, c('MC','MC.diff',
#                                                         "Mantel",
#                                                         "Mantel.diff")])
# 
# rownames(compare16a) <- compare16$Scenario[1 + 1:nScenarios16]

## ----eval = FALSE-------------------------------------------------------------
# # Load in projections
# data("projections")
# 
# set.seed(75)
# capLocs14 <- lapply(breedingPos14,
#                     FUN = function(x){
#                       colnames(x)<-c("Longitude","Latitude")
#                       sf::st_as_sf(data.frame(x),
#                                    coords = c("Longitude","Latitude"),
#                                    crs = 4326)})
# 
# targLocs14 <- lapply(winteringPos14,
#                      FUN = function(x){
#                       colnames(x)<-c("Longitude","Latitude")
#                       sf::st_as_sf(data.frame(x),
#                                    coords = c("Longitude","Latitude"),
#                                    crs = 4326)})
# 
# 
# # Relative abundance by scenario and breeding population
# breedingN14base <- rep(25, nBreeding14[1])
# breedingN14 <- lapply(breedingSiteTrans14, rowsum, x=breedingN14base)
# breedingRelN14 <- lapply(breedingN14, "/", sum(breedingN14base))
# 
# 
# MC14 <- 0.25
# 
# nSample14.1 <- 100 # Total number sampled per simulation
# 
# # How many sampled from each natural population (sampling scenarios separate
# # from definition scenarios)
# sampleBreeding14.1 <- list(round(breedingRelN14[[1]]*nSample14.1),
#                          c(rep(0,22),round(breedingRelN14[[5]][1]*nSample14.1),
#                            rep(0,4), round(breedingRelN14[[5]][2]*nSample14.1),
#                            rep(0,44),round(breedingRelN14[[5]][3]*nSample14.1),
#                            rep(0,4), round(breedingRelN14[[5]][4]*nSample14.1),
#                            rep(0, 22)),
#                            rep(c(rep(0, 4),
#                            rep(round(breedingRelN14[[2]][1]*nSample14.1/2), 2),
#                            rep(0, 4)), 10))
# 
# lapply(sampleBreeding14.1, matrix, nrow=10, ncol=10)
# 
# # Number of simulations to run
# nSims14 <- 100
# nSimsLarge14 <- 250 #0riginally 2500
# nYears <- 1
# nMonths <- 1
# 
# # Set up data structures for storing results
# animalLoc14 <- vector("list", nScenarios14) #making an empty list to fill
# sim14 <- vector("list", nSimsLarge14) #making an empty list to fill
# 
# compare14.1.array <- array(NA, c(nSimsLarge14, 2),
#                            dimnames = list(1:nSimsLarge14,
#                                            c("MC", "Mantel")))
# 
# results14 <- vector("list", nScenarios14)
# 
# # Run simulations
# set.seed(7)
# 
# system.time(for (sim in 1:nSimsLarge14) {
# 
#   cat("Simulation", sim, "of", nSimsLarge14, '\n')
# 
#   sim14[[sim]] <- lapply(sampleBreeding14.1,
#                          simMove,
#                          breedingDist = breedDist14[[1]],
#                          winteringDist=nonbreedDist14[[1]],
#                          psi=psi,
#                          nYears=nYears,
#                          nMonths=nMonths)
# 
# for (i in 13) {
#  animalLoc14[[i]] <- changeLocations(sim14[[sim]][[scenarioToSampleMap14[i]]]$animalLoc,
#                                      breedingSiteTrans14[[i]],
#                                      winteringSiteTrans14[[i]])
# 
#  results14[[i]] <- calcPsiMC(breedDist14[[i]],
#                              nonbreedDist14[[i]],
#                              animalLoc14[[i]],
#                              originRelAbund = breedingRelN14[[i]],
#                              verbose = FALSE)
# 
#   compare14.1.array[sim,'MC'] <- results14[[i]]$MC
# 
#   compare14.1.array[sim,'Mantel'] <- calcStrengthInd(breedDist14[[1]],
#                                                   nonbreedDist14[[1]],
#                             sim14[[sim]][[scenarioToSampleMap14[i]]]$animalLoc,
#                                                   resamp=0)$correlation
#   }
# })
# 
# means14 <- apply(compare14.1.array, 2, mean)
# 
# vars14 <- apply(compare14.1.array, 2, var)
# 
# rmse14 <- apply(compare14.1.array, 2, function(x) sqrt(mean((x - MC14)^2)))
# 
# # Set up data structures for storing estimation results
# est14.array <- array(NA, c(nSims14, 2),
#                      dimnames = list(1:nSims14,c("MC", "Mantel")))
# 
# var14.array <- array(NA, c(nSims14, 2),
#                      dimnames =list(1:nSims14,c("MC", "Mantel")))
# 
# ci14.array <- array(NA, c(nSims14, 2, 2),
#                     dimnames = list(1:nSims14,
#                                     c("MC", "Mantel"),
#                                     c('lower', 'upper')))
# 
# #making an empty list to fill
# 
# animalLoc14 <- vector("list", nSims14)
# 
# results14 <- vector("list", nSims14)
# 
# # Run estimations
# set.seed(567)
# 
# sim14.sub <- sim14[sample.int(nSimsLarge14, nSims14, TRUE)]
# 
# for (sim in 1:nSims14) {
#   cat("Estimation", sim, "of", nSims14, '\n')
# 
#  for (i in 13) {#:nScenarios14) {
#   animalLoc14[[sim]] <-
#     changeLocations(sim14.sub[[sim]][[scenarioToSampleMap14[i]]]$animalLoc,
#                                         breedingSiteTrans14[[i]],
#                                         winteringSiteTrans14[[i]])
# 
#   results14[[sim]] <- estMC(originDist = breedDist14[[i]],
#                             targetDist = nonbreedDist14[[i]],
#                             originRelAbund = breedingRelN14[[i]],
#               originPoints = capLocs14[[i]][animalLoc14[[sim]][, 1, 1, 1], ],
#               targetPoints = targLocs14[[i]][animalLoc14[[sim]][, 2, 1, 1], ],
#                             originAssignment = animalLoc14[[sim]][, 1, 1, 1],
#                             targetAssignment = animalLoc14[[sim]][, 2, 1, 1],
#                             nSamples = 1000,
#                             verbose = 1,
#                             calcCorr = TRUE,
#                             geoBias = c(0, 0),
#                             geoVCov = matrix(0, 2, 2))
# 
#     est14.array[sim, 'MC'] <- results14[[sim]]$MC$mean
#     est14.array[sim, 'Mantel'] <- results14[[sim]]$corr$mean
#     var14.array[sim, 'MC'] <- results14[[sim]]$MC$se ^ 2
#     var14.array[sim, 'Mantel'] <- results14[[sim]]$corr$se ^ 2
#     ci14.array[sim, 'MC', ] <- results14[[sim]]$MC$bcCI
#     ci14.array[sim, 'Mantel', ] <- results14[[sim]]$corr$bcCI
#   }
# }
# 
# # Crude point estimates (bootstrap means are better)
# pointEsts14.1 <- t(sapply(results14, function(x) c(x$MC$point, x$corr$point)))
# 
# # Summarize results
# summary(var14.array)
# vars14
# summary(est14.array)
# summary(pointEsts14.1)
# means14
# colMeans(pointEsts14.1) - MC14
# sqrt(colMeans((pointEsts14.1 - MC14)^2))
# summary(ci14.array[, "MC", 'lower'] <= MC14 &
#           ci14.array[, "MC", 'upper'] >= MC14)
# summary(ci14.array[, "Mantel", 'lower'] <= MC14 &
#           ci14.array[, "Mantel", 'upper'] >= MC14)
# 
# # Plot histograms of bootstrap estimation results
# est.df <- data.frame(Parameter = rep(c("MC", 'rM'), 2, each = nSims14),
#                      Quantity = rep(c('Mean', 'Variance'), each = 2 * nSims14),
#                      sim = rep(1:nSims14, 4),
#                      Estimate = c(est14.array, var14.array))
# 
# trues.df <- data.frame(Parameter = rep(c("MC", 'rM'), 2),
#                        Quantity = rep(c('Mean', 'Variance'), each = 2),
#                        Value = c(MC14, MC14, vars14))
# 
# library(ggplot2)
# g.est <- ggplot(est.df, aes(Estimate)) + geom_histogram(bins = 15) +
#   facet_grid(Parameter ~ Quantity, scales = 'free_x') +
#   geom_vline(aes(xintercept = Value), data = trues.df) +
#   theme_bw() + scale_x_continuous(breaks = c(0.002, 0.003, 0.1, 0.2, 0.3))
# 
# g.est
# 
# 
# # Summary table of bootstrap estimation results
# qualities14 <- data.frame(Parameter = rep(c("MC", 'rM'), 5),
#                           Quantity = rep(rep(c('Mean', 'Variance'), c(2, 2)),
#                                          3, 10),
#                           Measure = rep(c('Bias', 'RMSE', "Coverage"),
#                                         c(4, 4, 2)),
#                           Value = c(colMeans(est14.array - MC14),
#                                     colMeans(var14.array) - vars14,
#                                     sqrt(colMeans((est14.array - MC14)^2)),
#                                     sqrt(mean((var14.array[,1]-vars14[1])^2)),
#                                     sqrt(mean((var14.array[,2]-vars14[2])^2)),
#                                     mean(ci14.array[, "MC", 'lower'] <= MC14 &
#                                          ci14.array[, "MC", 'upper'] >= MC14),
#                                     mean(ci14.array[,"Mantel",'lower']<=MC14 &
#                                          ci14.array[,"Mantel",'upper']>=MC14)))
# 
# format(qualities14, digits = 2, scientific = FALSE)
# 

## ----eval = FALSE-------------------------------------------------------------
# # Simulation using error associated with light-level geolocators
# # Assign geolocator bias / variance co-variance matrix
# 
# geoBias <- c(-10000, 50000)
# geoVCov <- matrix(c(1.2e+8, 0, 0, 1.2e+9), 2, 2)
# sqrt(diag(geoVCov))
# 
# targLocs14.2 <- lapply(targLocs14,
#                        FUN = function(x){sf::st_transform(x, 6933)})
# 
# # convert wintering locations to polygons using the helper function
# WinteringPolys <- lapply(winteringPos14,
#                          toPoly,
#                          projection.in = 4326,
#                          projection.out = 6933)
# 
# winteringPolys <- lapply(winteringPos14,
#                          toPoly,
#                          projection.in = 4326,
#                          projection.out = 4326)
# 
# winteringPolys2 <- lapply(winteringPos14,
#                           toPoly,
#                           projection.in = 4326,
#                           projection.out = 'ESRI:54027')
# 
# 
# # Simulate capture and non-breeding locations of 100 individuals #
# 
# nSample14.2 <- 100 # Total number sampled per simulation
# 
# # How many sampled from each natural population (sampling scenarios separate
# # from definition scenarios)
# sampleBreeding14.2 <- list(round(breedingRelN14[[1]]*nSample14.2),
#                          c(rep(0,22),round(breedingRelN14[[5]][1]*nSample14.2),
#                            rep(0, 4),round(breedingRelN14[[5]][2]*nSample14.2),
#                            rep(0,44),round(breedingRelN14[[5]][3]*nSample14.2),
#                            rep(0, 4),round(breedingRelN14[[5]][4]*nSample14.2),
#                            rep(0, 22)),
#                          rep(c(rep(0, 4),
#                                rep(round(breedingRelN14[[2]][1]*nSample14.2/2),
#                                    2),
#                                rep(0, 4)), 10))
# 
# lapply(sampleBreeding14.2, matrix, nrow=10, ncol=10)
# 
# # Number of simulations to run
# nSims14.2 <- 10
# nSimsLarge14.2 <- 25
# nYears <- 1
# nMonths <- 1
# 
# # Set up data structures for storing results
# 
# #making an empty list to fill
# animalLoc14 <- vector("list", nScenarios14)
# sim14.2 <- vector("list", nSimsLarge14.2)
# 
# compare14.2.array <- array(NA, c(nSimsLarge14.2, 2),
#                            dimnames = list(1:nSimsLarge14.2,
#                                            c("MC", "Mantel")))
# 
# results14 <- vector("list", nScenarios14)
# results14.2 <- vector("list", nSimsLarge14.2)
# 
# # Run simulations
# set.seed(7)
# 
# system.time(for (sim in 1:nSimsLarge14.2) {
#   cat("Simulation", sim, "of", nSimsLarge14.2, '\n')
#   sim14.2[[sim]] <- lapply(sampleBreeding14.2,
#                            simMove,
#                            breedingDist = breedDist14[[1]],
#                            winteringDist=nonbreedDist14[[1]],
#                            psi=psi,
#                            nYears=nYears,
#                            nMonths=nMonths)
# 
# for (i in 13) { # only run one scenario
#   animalLoc14.0 <- sim14.2[[sim]][[scenarioToSampleMap14[i]]]$animalLoc
#   animalLoc14[[i]] <- changeLocations(animalLoc14.0,
#                                       breedingSiteTrans14[[i]],
#                                       winteringSiteTrans14[[i]])
# 
#     breedingTruePoints14 <- capLocs14[[i]][animalLoc14[[i]][,1,1,1], ]
# 
#     winteringTruePoints14 <- targLocs14.2[[i]][animalLoc14.0[,2,1,1], ]
# 
#     results14.2[[sim]] <- simLocationError(winteringTruePoints14,
#                                            WinteringPolys[[i]],
#                                            geoBias,
#                                            geoVCov,
#                                            projection = sf::st_crs(winteringTruePoints14))
# 
#     animalLoc14.2 <- animalLoc14[[i]]
# 
#     animalLoc14.2[ , 2, 1, 1] <- results14.2[[sim]]$targetSample
# 
#     results14[[i]] <- calcPsiMC(breedDist14[[i]],
#                                 nonbreedDist14[[i]],
#                                 animalLoc14.2,
#                                 originRelAbund = breedingRelN14[[i]],
#                                 verbose = FALSE)
# 
#     compare14.2.array[sim, 'MC'] <- results14[[i]]$MC
# 
#     originDist1 <- distFromPos(sf::st_coordinates(breedingTruePoints14))
# 
#     target.point.sample <- sf::st_as_sf(
#       data.frame(results14.2[[sim]]$targetPointSample),
#       coords = c(1,2),
#       crs = 6933)
# 
#     target.point.sample2 <- sf::st_transform(target.point.sample, 4326)
# 
#     targetDist1 <-distFromPos(sf::st_coordinates(target.point.sample2))
# 
#     compare14.2.array[sim, 'Mantel'] <- calcMantel(originDist = originDist1,
#                                                    targetDist = targetDist1)$pointCorr
#   }
# })
# 
# means14.2 <- apply(compare14.2.array, 2, mean)
# vars14.2 <- apply(compare14.2.array, 2, var)
# rmse14.2 <- apply(compare14.2.array, 2, function(x) sqrt(mean((x - MC14)^2)))
# 
# # Set up data structures for storing estimation results
# est14.2.array <- array(NA, c(nSims14.2, 2),
#                        dimnames = list(1:nSims14.2,
#                                        c("MC", "Mantel")))
# 
# var14.2.array <- array(NA, c(nSims14.2, 2),
#                        dimnames = list(1:nSims14.2,
#                                        c("MC", "Mantel")))
# 
# ci14.2.array <- array(NA, c(nSims14.2, 2, 2),
#                       dimnames = list(1:nSims14.2,
#                                       c("MC", "Mantel"),
#                                       c('lower', 'upper')))
# 
# animalLoc14 <- vector("list", nSims14.2) #making an empty list to fill
# results14 <- vector("list", nSims14.2)
# 
# # Run estimations
# 
# set.seed(567)
# sim14.2.sub <- sim14.2[sample.int(nSimsLarge14.2, nSims14.2, TRUE)]
# 
# for (sim in 1:nSims14.2) {
#   cat("Estimation", sim, "of", nSims14.2, '\n')
#  for (i in 13) {
#   points0 <- sf::st_as_sf(data.frame(results14.2[[sim]]$targetPointSample),
#                           coords=c("X1","X2"),
#                           crs = 6933)
# 
#   points1 <- sf::st_transform(points0, crs = "ESRI:54027")
# 
#   animalLoc14[[sim]] <- changeLocations(
#     sim14.2.sub[[sim]][[scenarioToSampleMap14[i]]]$animalLoc,
#                                         breedingSiteTrans14[[i]],
#                                         winteringSiteTrans14[[i]])
# 
#   results14[[sim]] <- estMC(originDist = breedDist14[[i]],
#                             targetDist = nonbreedDist14[[i]],
#                             originRelAbund = breedingRelN14[[i]],
#                 originPoints = capLocs14[[i]][animalLoc14[[sim]][, 1, 1, 1], ],
#                             targetPoints = points1,
#                             originAssignment = animalLoc14[[sim]][, 1, 1, 1],
#                             targetAssignment = results14.2[[sim]]$targetSample,
#                             nSamples = 1000,
#                             verbose = 0,
#                             calcCorr = TRUE,
#                             geoBias = geoBias,
#                             geoVCov = geoVCov,
#                             isGL = TRUE,
#                             targetSites = winteringPolys2[[i]],
#                             maintainLegacyOutput = TRUE)
# 
#     est14.2.array[sim, 'MC'] <- results14[[sim]]$meanMC
# 
#     est14.2.array[sim, 'Mantel'] <- results14[[sim]]$meanCorr
# 
#     var14.2.array[sim, 'MC'] <- results14[[sim]]$seMC ^ 2
# 
#     var14.2.array[sim, 'Mantel'] <- results14[[sim]]$seCorr ^ 2
# 
#     ci14.2.array[sim, 'MC', ] <- results14[[sim]]$bcCI
# 
#     ci14.2.array[sim, 'Mantel', ] <- results14[[sim]]$bcCICorr
# 
#   }
# }
# 
# pointEsts14.2 <- t(sapply(results14, function(x) c(x$pointMC, x$pointCorr)))
# 
# # Summarize estimation results
# summary(var14.2.array)
# vars14.2
# summary(est14.2.array)
# summary(pointEsts14.2)
# colMeans(pointEsts14.2) - MC14
# sqrt(colMeans((pointEsts14.2 - MC14)^2))
# means14.2
# mean(ci14.2.array[, "MC", 'lower'] <= MC14 &
#        ci14.2.array[, "MC", 'upper'] >= MC14)
# mean(ci14.2.array[, "Mantel", 'lower'] <= MC14 &
#        ci14.2.array[, "Mantel", 'upper'] >= MC14)
# sqrt(vars14.2) / MC14
# summary(sqrt(var14.2.array)/est14.2.array)
# 
# # Plot histograms of bootstrap estimation results
# est14.2.df <- data.frame(Parameter = rep(c("MC", 'rM'), 2, each = nSims14.2),
#                      Quantity = rep(c('Mean', 'Variance'),
#                                     each = 2 * nSims14.2),
#                      sim = rep(1:nSims14.2, 4),
#                      Estimate = c(est14.2.array, var14.2.array))
# 
# trues14.2.df <- data.frame(Parameter = rep(c("MC", 'rM'), 2),
#                        Quantity = rep(c('Mean', 'Variance'), each = 2),
#                        Value = c(MC14, MC14, vars14.2))
# 
# g.est <- ggplot(est14.2.df, aes(Estimate)) + geom_histogram(bins = 15) +
#   facet_grid(Parameter ~ Quantity, scales = 'free_x') +
#   geom_vline(aes(xintercept = Value), data = trues14.2.df) +
#   theme_bw() + scale_x_continuous(breaks = c(0.002, 0.003, 0.1, 0.2, 0.3))
# 
# g.est
# 
# 
# # Summary table of bootstrap estimation results
# qualities14.2 <- data.frame(Parameter = rep(c("MC", 'rM'), 5),
#                             Quantity = rep(rep(c('Mean', 'Variance'),
#                                                c(2, 2)), 3, 10),
#                             Measure = rep(c('Bias', 'RMSE', "Coverage"),
#                                           c(4, 4, 2)),
#                             Value = c(colMeans(est14.2.array - MC14),
#                                       colMeans(var14.2.array) - vars14.2,
#                                       sqrt(colMeans((est14.2.array - MC14)^2)),
#                                       sqrt(mean((var14.2.array[,1] -
#                                                    vars14[1])^2)),
#                                       sqrt(mean((var14.2.array[,2] -
#                                                    vars14[2])^2)),
#                                       mean(ci14.2.array[,"MC",'lower']<=MC14 &
#                                            ci14.2.array[,"MC",'upper']>=MC14),
#                               mean(ci14.2.array[, "Mantel", 'lower'] <= MC14 &
#                                    ci14.2.array[, "Mantel", 'upper'] >= MC14)))
# 
# format(qualities14.2, digits = 2, scientific = FALSE)
# 
# 
# 
# # Make summary table and figures of results with (Example 9) and without
# # (Example 8) location error
# 
# qualities.all <- rbind(data.frame(DataType = 'GPS', qualities14),
#                        data.frame(DataType = 'GL', qualities14.2))
# 
# 
# est14.all.df <- rbind(data.frame(DataType = 'GPS', est.df),
#                        data.frame(DataType = 'GL', est14.2.df))
# 
# trues14.all.df <- data.frame(DataType = rep(c("GPS", "GL"), 2, each = 2),
#                              Parameter = rep(c("MC", 'rM'), 4),
#                        Quantity = rep(c('Mean', 'Variance'), each = 4),
#                        Value = c(rep(MC14, 4), vars14, vars14.2))
# 
# g.est <- ggplot(est14.all.df, aes(Estimate)) + geom_histogram(bins = 12) +
#   facet_grid(DataType + Parameter ~ Quantity, scales = 'free_x') +
#   geom_vline(aes(xintercept = Value), data = trues14.all.df) +
#   theme_bw() + scale_x_continuous(breaks = c(0.002, 0.003, 0.1, 0.2, 0.3))
# 
# g.est
# 
# 

## ----eval = FALSE-------------------------------------------------------------
# ## Simulation ----
# nBreeding <- 100
# nWintering <- 100
# breedingPos <- matrix(c(rep(seq(-99, -81, 2), each=sqrt(nBreeding)),
#                         rep(seq(49, 31, -2), sqrt(nBreeding))), nBreeding, 2)
# winteringPos <- matrix(c(rep(seq(-79, -61, 2), each=sqrt(nWintering)),
#                          rep(seq(9, -9, -2), sqrt(nWintering))), nWintering, 2)
# head(breedingPos)
# tail(breedingPos)
# head(winteringPos)
# tail(winteringPos)
# 
# breedDist <- distFromPos(breedingPos, 'ellipsoid')
# nonbreedDist <- distFromPos(winteringPos, 'ellipsoid')
# 
# # Breeding Abundance
# breedingN <- rep(5000, nBreeding)
# breedingRelN <- breedingN/sum(breedingN)
# 
# # Set up psi matrix
# o <- optimize(mlogitMC, MC.in = 0.25, origin.dist = breedDist,
#               target.dist = nonbreedDist, origin.abund = breedingRelN,
#               sample.size = sum(breedingN),
#               interval = c(0, 10), tol = .Machine$double.eps^0.5)
# 
# slope <- o$minimum
# psi <- mlogitMat(slope, breedDist)
# 
# # Baseline strength of migratory connectivity
# MC <- calcMC(breedDist, nonbreedDist, breedingRelN, psi, sum(breedingN))
# MC
# 
# # Other basic simulation parameters
# 
# ## Dispersal simulations---
# set.seed(1516)
# 
# nYears <- 15
# 
# nMonths <- 4 # Each season
# 
# Drates <- c(0.02, 0.04, 0.08, 0.16, 0.32, 0.64)    #rates of dispersal
# 
# birdLocDisp <- vector('list', length(Drates))
# 
# Disp.df  <- data.frame(Year=rep(1:nYears, length(Drates)),
#                        Rate=rep(Drates, each = nYears), MC = NA)
# 
# for(i in 1:length(Drates)){
#   cat('Dispersal Rate', Drates[i], '\n')
#   birdLocDisp[[i]] <- simMove(breedingN,
#                               breedDist,
#                               nonbreedDist,
#                               psi,
#                               nYears,
#                               nMonths,
#                               sumDispRate = Drates[i])
#   for(j in 1:nYears){
#     cat('\tYear', j, '\n')
#     temp.results <- calcPsiMC(breedDist,
#                               nonbreedDist,
#                               breedingRelN,
#                               birdLocDisp[[i]]$animalLoc,
#                               years=j,
#                               verbose = FALSE)
# 
#     Disp.df$MC[j + (i - 1) * nYears] <- temp.results$MC
#   }
# } # end i loop
# 
# Disp.df$Year <- Disp.df$Year - 1 #just run once!
# 
# data.frame(Disp.df, roundMC = round(Disp.df$MC, 2),
#            nearZero = Disp.df$MC < 0.01)
# 
# 
# # Convert dispersal rates to probabilities of dispersing at least certain distance
# threshold <- 1000
# probFarDisp <- matrix(NA, nBreeding, length(Drates), dimnames = list(NULL, Drates))
# 
# for (i in 1:length(Drates)) {
#   for (k in 1:nBreeding) {
#     probFarDisp[k, i] <-
#       sum(birdLocDisp[[i]]$natalDispMat[k, which(breedDist[k, ]>= threshold)])
#   }
# }
# 
# summary(probFarDisp)
#