Nothing
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)
#
Try the MigConnectivity package in your browser
Any scripts or data that you put into this service are public.
MigConnectivity documentation built on May 29, 2024, 8:37 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## ----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)
#
Try the MigConnectivity package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.