inst/examples/examples.R
In SMBC-NZP/MigConnectivity: Estimate Migratory Connectivity for Migratory Animals
library(MigConnectivity)
set.seed(75)
###############################################################################
# Utility functions for use in simulations
###############################################################################
# Calculates probability matrix based on exponential decline with distance
mlogitMat <- function(slope, dist) {
preMat <- exp(-slope/mean(dist)*dist)
diag(preMat) <- 0
nr <- nrow(dist)
nc <- ncol(dist)
outMat <- matrix(0, nr, nc)
for (b in 1:nr) {
outMat[b,] <- preMat[b,]/(1+sum(preMat[b, ]))
outMat[b,b] <- 1 - sum(outMat[b, ])
}
return(outMat)
}
# Crude optimizable function for developing MC pattern based on MC strength
mlogitMC <- function(slope, MC.in, origin.dist, target.dist, origin.rel.abund) {
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, psi, origin.rel.abund)
return((MC.in - MC)^2)
}
# rho function for individuals
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] <- originDist2[j,i] <- originDist[locations[i,1,1,1], locations[j,1,1,1]]
targetDist2[i,j] <- targetDist2[j,i] <- targetDist[locations[i,2,1,1], locations[j,2,1,1]]
}
}
return(ncf::mantel.test(originDist2, targetDist2, resamp=resamp, quiet = !verbose))
}
###############################################################################
# rho function for individuals
###############################################################################
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] <- originDist2[j,i] <- originDist[locations[i,1,1,1],
locations[j,1,1,1]]
targetDist2[i,j] <- targetDist2[j,i] <- targetDist[locations[i,2,1,1],
locations[j,2,1,1]]
}
}
return(ncf::mantel.test(originDist2, targetDist2, resamp=resamp, quiet = !verbose))
}
###############################################################################
# Simple approach to estimate psi matrix and MC from simulated (or real) data
# (doesn't include uncertainty)
###############################################################################
calcPsiMC <- function(originDist, targetDist, originRelAbund, locations,
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,,]
targetMat <- locations[i,2,,]
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, psiMat, originRelAbund)
return(list(psi=psiMat, MC=MC))
}
###############################################################################
# Parameters for simulations
###############################################################################
\donttest{
nSeasons <- 2
nYears <- 10
nMonths <- 4 # Each season
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')
breedDist[1:12, 1:12]
breedDist[1:12, c(1,91,100)]
# Breeding Abundance
breedingN <- rep(500, nBreeding)
breedingRelN <- breedingN/sum(breedingN)
# Set up psi matrix
o <- optimize(mlogitMC, MC.in = 0.25, origin.dist = breedDist,
target.dist = nonbreedDist, origin.rel.abund = breedingRelN,
interval = c(0, 10), tol = .Machine$double.eps^0.5)
o
slope <- o$minimum
psi <- mlogitMat(slope, breedDist)
round(psi[1:12, 1:12],5)
rowSums(psi)
# Baseline strength of migratory connectivity
MC <- calcMC(breedDist, nonbreedDist, psi, breedingRelN)
MC
# Simulation
sim <- simMove(breedingN, breedDist, nonbreedDist, psi, nYears, nMonths)
###############################
# Sampling regime 1 of 3
# Researchers divide populations differently than reality
# Delineation of seasonal ranges into regions
#I) Breeding range divided along equal longitudinal breaks into ten regions
#II) Non-breeding range divided along equal longitudinal breaks into ten regions
#III) Breeding and non-breeding ranges divided along equal longitudinal breaks into ten regions
#IV) Breeding range divided along the longitudinal and latitudinal midpoint into four regions
#V) Non-breeding range divided along the longitudinal and latitudinal midpoint into four regions
#VI) Breeding range divided along the longitudinal and latitudinal midpoint into four regions and
#non-breeding range divided along equal longitudinal breaks into ten regions
#VII) Breeding range divided along equal longitudinal breaks into ten regions and
#non-breeding range divided along the longitudinal and latitudinal midpoint into four regions
###############################
#Run functions and parameters above first
set.seed(75)
#each element is for a scenario (see above 1-8)
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:nBreeding, each=10))
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)))
breedingSiteTrans14
lapply(breedingSiteTrans14, matrix, nrow=10, ncol=10)
lapply(winteringSiteTrans14, matrix, nrow=10, ncol=10)
breedingPos14 <- list(breedingPos, rowsum(breedingPos, rep(1:10, each=10))/10, #positions of the human defined populations
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)
breedingPos14
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)
winteringPos14
breedDist14 <- lapply(breedingPos14, distFromPos, surface = 'ellipsoid')
breedDist14
lapply(breedDist14, max)
nonbreedDist14 <- lapply(winteringPos14, distFromPos, surface = 'ellipsoid')
nonbreedDist14
lapply(nonbreedDist14, max)
nBreeding14 <- c(100, 10, 100, 10, 4, 100, 4, 10)
nWintering14 <- c(100, 100, 10, 10, 100, 4, 10, 4)
breedingRelN14 <- lapply(nBreeding14, function(x) rep(1/x, x))
breedingRelN14
nSample14 <- 1000 # Total number sampled per simulation
# for the baseline use the simulation from above, sim
animalLoc14base <- sim$animalLoc
#transferring the simulated bird locations from the true populations to the researcher defined populations
changeLocations <- function(animalLoc, breedingSiteTrans, winteringSiteTrans) {
animalLoc[,1,,] <- breedingSiteTrans[animalLoc[,1,,]]
animalLoc[,2,,] <- winteringSiteTrans[animalLoc[,2,,]]
return(animalLoc)
}
nScenarios14 <- length(breedingSiteTrans14)
animalLoc14 <- vector("list", nScenarios14) #making an empty list to fill
for (i in 1:nScenarios14)
animalLoc14[[i]] <- changeLocations(animalLoc14base, breedingSiteTrans14[[i]], winteringSiteTrans14[[i]])
animalLoc14[[2]][101:110,1,1,1]
animalLoc14[[6]][101:110,1,1,1]
animalLoc14[[6]][101:110,2,1,1]
results14 <- vector("list", nScenarios14)
compare14 <- data.frame(Scenario = c("True", "Base", "Breeding10", "Wintering10",
"Breeding10Wintering10", "Breeding4",
"Wintering4", "Breeding4Wintering10",
"Breeding10Wintering4"),
MC = c(MC, rep(NA, nScenarios14)))
for (i in 1:nScenarios14) {
cat("\nScenario", i, "\n")
results14[[i]] <- calcPsiMC(breedDist14[[i]], nonbreedDist14[[i]],
breedingRelN14[[i]], animalLoc14[[i]], TRUE)
compare14$MC[i+1] <- results14[[i]]$MC
}
compare14 <- transform(compare14, diff=MC - MC[1], prop=MC/MC[1])
compare14
###############################################################################
# Sampling regime 2 of 3
# Researchers divide populations differently than reality PLUS
# Different distributions of sampling animals across breeding range PLUS
# Sample sizes don't always match relative abundances PLUS
# Compare our approach and simple Mantel approach
# 1. Base (10 years, uneven abundances but matches sampling)
# 2. Breeding pops divided into 4 squares, sample across breeding range
# 3. Breeding pops divided into 4 squares, sample at centroid of each square
# 4. Sampling high in low abundance populations plus base
# 5. Scenarios 2 plus 4
# 6. Scenarios 3 plus 4
#May want to create another set because we have varied these two things together
###############################################################################
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)))
breedingSiteTrans15
lapply(breedingSiteTrans15, matrix, nrow=10, ncol=10)
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]])
breedingPos15
winteringPos15 <- rep(list(winteringPos), nScenarios15)
winteringPos15
breedDist15 <- lapply(breedingPos15, distFromPos)
breedDist15[[1]][1,]
nonbreedDist15 <- lapply(winteringPos15, distFromPos)
nBreeding15 <- rep(c(100, 4, 4), 2)
nBreeding15
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
breedingN15base
matrix(breedingN15base, 10, 10)
sum(breedingN15base)
breedingN15 <- lapply(breedingSiteTrans15, rowsum, x=breedingN15base) # For researcher defined populations
breedingRelN15 <- lapply(breedingN15, "/", sum(breedingN15base))
nSample15 <- 1000 # Total number sampled per simulation
# Four sampling regimes/simulations, because repeat a couple of them in different scenarios
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])
lapply(sampleBreeding15, matrix, nrow=10, ncol=10)
lapply(sampleBreeding15, sum)
# Set up psi matrix
o15 <- optimize(mlogitMC, MC.in = 0.25, origin.dist = breedDist15[[1]],
target.dist = nonbreedDist15[[1]],
origin.rel.abund = breedingRelN15[[1]], interval = c(0,10),
tol = .Machine$double.eps^0.5)
o15
slope15 <- o15$minimum
psi15 <- mlogitMat(slope15, breedDist15[[1]])
round(psi15[1:12, 1:12],5)
rowSums(psi15)
# Baseline strength of migratory connectivity
MC15 <- calcMC(breedDist15[[1]], nonbreedDist15[[1]], psi15, breedingRelN15[[1]])
MC15
# 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(sim15[[scenarioToSampleMap15[i]]]$animalLoc,
breedingSiteTrans15[[i]], 1:nWintering15[i])
results15[[i]] <- calcPsiMC(breedDist15[[i]], nonbreedDist15[[i]],
breedingRelN15[[i]], animalLoc15[[i]], 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)
compare15$Mantel[1:nScenarios15 + 1] <- apply(compare15.array[,,'Mantel'], 2,
mean)
compare15 <- transform(compare15, MC.diff=MC - MC[1],
Mantel.diff=Mantel - Mantel[1],
MC.prop=MC/MC[1], Mantel.prop=Mantel/Mantel[1])
compare15
compare15a <- as.matrix(compare15[2:7,c(2,4,3,5)])
rownames(compare15a) <- compare15$Scenario[2:7]
round(compare15a, 3)
round(compare15a, 2)
###############################################################################
# Sampling regime 3 of 3
# Researchers divide populations differently than reality (simulations) PLUS
# Different distributions of sampled animals across breeding range PLUS
# Sample sizes don't always match relative abundances PLUS
# Compare our approach and simple Mantel approach PLUS
# MC not same across subsections of range
# 1. Base (uneven MC (0.15 for NW breeding, 0.3 for SW, 0.45 for NE, and 0.6
# for SE), uneven abundances (lowest in NW, highest in SE), sampling
# proportional to abundance
# 2. Breeding pops divided into 4 squares, sample across breeding range
# 3. Breeding pops divided into 4 squares, sample at centroid of each square
# 4. Sampling high in low abundance populations
# 5. Scenarios 2 plus 4
# 6. Scenarios 3 plus 4
###############################################################################
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))
lapply(breedingRelN16, sum)
# 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.rel.abund = breedingRelN16[[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, ]
}
round(psi16[1:12, 1:12],5)
rowSums(psi16)
# Baseline strength of migratory connectivity
MC16 <- calcMC(breedDist16[[1]], nonbreedDist16[[1]], psi16, breedingRelN16[[1]])
MC16
# Set up sampling regimes (different number than number of scenarios)
nSample16 <- 1000
sampleBreeding16 <- sampleBreeding15
lapply(sampleBreeding16, matrix, nrow=10, ncol=10)
lapply(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")))
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(sim15[[scenarioToSampleMap15[i]]]$animalLoc,
breedingSiteTrans15[[i]], 1:nWintering15[i])
results15[[i]] <- calcPsiMC(breedDist15[[i]], nonbreedDist15[[i]],
breedingRelN15[[i]], animalLoc15[[i]], 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
}
}
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(breedDist16[[i]], nonbreedDist16[[i]],
breedingRelN16[[i]], animalLoc16[[i]], 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)
compare16$Mantel[1:nScenarios16 + 1] <- apply(compare16.array[,,'Mantel'], 2, mean)
compare16 <- transform(compare16, MC.diff=MC - MC[1], Mantel.diff=Mantel - Mantel[1],
MC.prop=MC/MC[1], Mantel.prop=Mantel/Mantel[1])
compare16
compare16a <- as.matrix(compare16[2:7,c(2,4,3,5)])
rownames(compare16a) <- compare16$Scenario[2:7]
round(compare16a, 3)
round(compare16a, 2)
}
SMBC-NZP/MigConnectivity documentation built on March 26, 2024, 4:22 p.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.
library(MigConnectivity)
set.seed(75)
###############################################################################
# Utility functions for use in simulations
###############################################################################
# Calculates probability matrix based on exponential decline with distance
mlogitMat <- function(slope, dist) {
preMat <- exp(-slope/mean(dist)*dist)
diag(preMat) <- 0
nr <- nrow(dist)
nc <- ncol(dist)
outMat <- matrix(0, nr, nc)
for (b in 1:nr) {
outMat[b,] <- preMat[b,]/(1+sum(preMat[b, ]))
outMat[b,b] <- 1 - sum(outMat[b, ])
}
return(outMat)
}
# Crude optimizable function for developing MC pattern based on MC strength
mlogitMC <- function(slope, MC.in, origin.dist, target.dist, origin.rel.abund) {
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, psi, origin.rel.abund)
return((MC.in - MC)^2)
}
# rho function for individuals
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] <- originDist2[j,i] <- originDist[locations[i,1,1,1], locations[j,1,1,1]]
targetDist2[i,j] <- targetDist2[j,i] <- targetDist[locations[i,2,1,1], locations[j,2,1,1]]
}
}
return(ncf::mantel.test(originDist2, targetDist2, resamp=resamp, quiet = !verbose))
}
###############################################################################
# rho function for individuals
###############################################################################
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] <- originDist2[j,i] <- originDist[locations[i,1,1,1],
locations[j,1,1,1]]
targetDist2[i,j] <- targetDist2[j,i] <- targetDist[locations[i,2,1,1],
locations[j,2,1,1]]
}
}
return(ncf::mantel.test(originDist2, targetDist2, resamp=resamp, quiet = !verbose))
}
###############################################################################
# Simple approach to estimate psi matrix and MC from simulated (or real) data
# (doesn't include uncertainty)
###############################################################################
calcPsiMC <- function(originDist, targetDist, originRelAbund, locations,
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,,]
targetMat <- locations[i,2,,]
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, psiMat, originRelAbund)
return(list(psi=psiMat, MC=MC))
}
###############################################################################
# Parameters for simulations
###############################################################################
\donttest{
nSeasons <- 2
nYears <- 10
nMonths <- 4 # Each season
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')
breedDist[1:12, 1:12]
breedDist[1:12, c(1,91,100)]
# Breeding Abundance
breedingN <- rep(500, nBreeding)
breedingRelN <- breedingN/sum(breedingN)
# Set up psi matrix
o <- optimize(mlogitMC, MC.in = 0.25, origin.dist = breedDist,
target.dist = nonbreedDist, origin.rel.abund = breedingRelN,
interval = c(0, 10), tol = .Machine$double.eps^0.5)
o
slope <- o$minimum
psi <- mlogitMat(slope, breedDist)
round(psi[1:12, 1:12],5)
rowSums(psi)
# Baseline strength of migratory connectivity
MC <- calcMC(breedDist, nonbreedDist, psi, breedingRelN)
MC
# Simulation
sim <- simMove(breedingN, breedDist, nonbreedDist, psi, nYears, nMonths)
###############################
# Sampling regime 1 of 3
# Researchers divide populations differently than reality
# Delineation of seasonal ranges into regions
#I) Breeding range divided along equal longitudinal breaks into ten regions
#II) Non-breeding range divided along equal longitudinal breaks into ten regions
#III) Breeding and non-breeding ranges divided along equal longitudinal breaks into ten regions
#IV) Breeding range divided along the longitudinal and latitudinal midpoint into four regions
#V) Non-breeding range divided along the longitudinal and latitudinal midpoint into four regions
#VI) Breeding range divided along the longitudinal and latitudinal midpoint into four regions and
#non-breeding range divided along equal longitudinal breaks into ten regions
#VII) Breeding range divided along equal longitudinal breaks into ten regions and
#non-breeding range divided along the longitudinal and latitudinal midpoint into four regions
###############################
#Run functions and parameters above first
set.seed(75)
#each element is for a scenario (see above 1-8)
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:nBreeding, each=10))
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)))
breedingSiteTrans14
lapply(breedingSiteTrans14, matrix, nrow=10, ncol=10)
lapply(winteringSiteTrans14, matrix, nrow=10, ncol=10)
breedingPos14 <- list(breedingPos, rowsum(breedingPos, rep(1:10, each=10))/10, #positions of the human defined populations
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)
breedingPos14
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)
winteringPos14
breedDist14 <- lapply(breedingPos14, distFromPos, surface = 'ellipsoid')
breedDist14
lapply(breedDist14, max)
nonbreedDist14 <- lapply(winteringPos14, distFromPos, surface = 'ellipsoid')
nonbreedDist14
lapply(nonbreedDist14, max)
nBreeding14 <- c(100, 10, 100, 10, 4, 100, 4, 10)
nWintering14 <- c(100, 100, 10, 10, 100, 4, 10, 4)
breedingRelN14 <- lapply(nBreeding14, function(x) rep(1/x, x))
breedingRelN14
nSample14 <- 1000 # Total number sampled per simulation
# for the baseline use the simulation from above, sim
animalLoc14base <- sim$animalLoc
#transferring the simulated bird locations from the true populations to the researcher defined populations
changeLocations <- function(animalLoc, breedingSiteTrans, winteringSiteTrans) {
animalLoc[,1,,] <- breedingSiteTrans[animalLoc[,1,,]]
animalLoc[,2,,] <- winteringSiteTrans[animalLoc[,2,,]]
return(animalLoc)
}
nScenarios14 <- length(breedingSiteTrans14)
animalLoc14 <- vector("list", nScenarios14) #making an empty list to fill
for (i in 1:nScenarios14)
animalLoc14[[i]] <- changeLocations(animalLoc14base, breedingSiteTrans14[[i]], winteringSiteTrans14[[i]])
animalLoc14[[2]][101:110,1,1,1]
animalLoc14[[6]][101:110,1,1,1]
animalLoc14[[6]][101:110,2,1,1]
results14 <- vector("list", nScenarios14)
compare14 <- data.frame(Scenario = c("True", "Base", "Breeding10", "Wintering10",
"Breeding10Wintering10", "Breeding4",
"Wintering4", "Breeding4Wintering10",
"Breeding10Wintering4"),
MC = c(MC, rep(NA, nScenarios14)))
for (i in 1:nScenarios14) {
cat("\nScenario", i, "\n")
results14[[i]] <- calcPsiMC(breedDist14[[i]], nonbreedDist14[[i]],
breedingRelN14[[i]], animalLoc14[[i]], TRUE)
compare14$MC[i+1] <- results14[[i]]$MC
}
compare14 <- transform(compare14, diff=MC - MC[1], prop=MC/MC[1])
compare14
###############################################################################
# Sampling regime 2 of 3
# Researchers divide populations differently than reality PLUS
# Different distributions of sampling animals across breeding range PLUS
# Sample sizes don't always match relative abundances PLUS
# Compare our approach and simple Mantel approach
# 1. Base (10 years, uneven abundances but matches sampling)
# 2. Breeding pops divided into 4 squares, sample across breeding range
# 3. Breeding pops divided into 4 squares, sample at centroid of each square
# 4. Sampling high in low abundance populations plus base
# 5. Scenarios 2 plus 4
# 6. Scenarios 3 plus 4
#May want to create another set because we have varied these two things together
###############################################################################
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)))
breedingSiteTrans15
lapply(breedingSiteTrans15, matrix, nrow=10, ncol=10)
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]])
breedingPos15
winteringPos15 <- rep(list(winteringPos), nScenarios15)
winteringPos15
breedDist15 <- lapply(breedingPos15, distFromPos)
breedDist15[[1]][1,]
nonbreedDist15 <- lapply(winteringPos15, distFromPos)
nBreeding15 <- rep(c(100, 4, 4), 2)
nBreeding15
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
breedingN15base
matrix(breedingN15base, 10, 10)
sum(breedingN15base)
breedingN15 <- lapply(breedingSiteTrans15, rowsum, x=breedingN15base) # For researcher defined populations
breedingRelN15 <- lapply(breedingN15, "/", sum(breedingN15base))
nSample15 <- 1000 # Total number sampled per simulation
# Four sampling regimes/simulations, because repeat a couple of them in different scenarios
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])
lapply(sampleBreeding15, matrix, nrow=10, ncol=10)
lapply(sampleBreeding15, sum)
# Set up psi matrix
o15 <- optimize(mlogitMC, MC.in = 0.25, origin.dist = breedDist15[[1]],
target.dist = nonbreedDist15[[1]],
origin.rel.abund = breedingRelN15[[1]], interval = c(0,10),
tol = .Machine$double.eps^0.5)
o15
slope15 <- o15$minimum
psi15 <- mlogitMat(slope15, breedDist15[[1]])
round(psi15[1:12, 1:12],5)
rowSums(psi15)
# Baseline strength of migratory connectivity
MC15 <- calcMC(breedDist15[[1]], nonbreedDist15[[1]], psi15, breedingRelN15[[1]])
MC15
# 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(sim15[[scenarioToSampleMap15[i]]]$animalLoc,
breedingSiteTrans15[[i]], 1:nWintering15[i])
results15[[i]] <- calcPsiMC(breedDist15[[i]], nonbreedDist15[[i]],
breedingRelN15[[i]], animalLoc15[[i]], 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)
compare15$Mantel[1:nScenarios15 + 1] <- apply(compare15.array[,,'Mantel'], 2,
mean)
compare15 <- transform(compare15, MC.diff=MC - MC[1],
Mantel.diff=Mantel - Mantel[1],
MC.prop=MC/MC[1], Mantel.prop=Mantel/Mantel[1])
compare15
compare15a <- as.matrix(compare15[2:7,c(2,4,3,5)])
rownames(compare15a) <- compare15$Scenario[2:7]
round(compare15a, 3)
round(compare15a, 2)
###############################################################################
# Sampling regime 3 of 3
# Researchers divide populations differently than reality (simulations) PLUS
# Different distributions of sampled animals across breeding range PLUS
# Sample sizes don't always match relative abundances PLUS
# Compare our approach and simple Mantel approach PLUS
# MC not same across subsections of range
# 1. Base (uneven MC (0.15 for NW breeding, 0.3 for SW, 0.45 for NE, and 0.6
# for SE), uneven abundances (lowest in NW, highest in SE), sampling
# proportional to abundance
# 2. Breeding pops divided into 4 squares, sample across breeding range
# 3. Breeding pops divided into 4 squares, sample at centroid of each square
# 4. Sampling high in low abundance populations
# 5. Scenarios 2 plus 4
# 6. Scenarios 3 plus 4
###############################################################################
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))
lapply(breedingRelN16, sum)
# 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.rel.abund = breedingRelN16[[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, ]
}
round(psi16[1:12, 1:12],5)
rowSums(psi16)
# Baseline strength of migratory connectivity
MC16 <- calcMC(breedDist16[[1]], nonbreedDist16[[1]], psi16, breedingRelN16[[1]])
MC16
# Set up sampling regimes (different number than number of scenarios)
nSample16 <- 1000
sampleBreeding16 <- sampleBreeding15
lapply(sampleBreeding16, matrix, nrow=10, ncol=10)
lapply(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")))
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(sim15[[scenarioToSampleMap15[i]]]$animalLoc,
breedingSiteTrans15[[i]], 1:nWintering15[i])
results15[[i]] <- calcPsiMC(breedDist15[[i]], nonbreedDist15[[i]],
breedingRelN15[[i]], animalLoc15[[i]], 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
}
}
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(breedDist16[[i]], nonbreedDist16[[i]],
breedingRelN16[[i]], animalLoc16[[i]], 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)
compare16$Mantel[1:nScenarios16 + 1] <- apply(compare16.array[,,'Mantel'], 2, mean)
compare16 <- transform(compare16, MC.diff=MC - MC[1], Mantel.diff=Mantel - Mantel[1],
MC.prop=MC/MC[1], Mantel.prop=Mantel/Mantel[1])
compare16
compare16a <- as.matrix(compare16[2:7,c(2,4,3,5)])
rownames(compare16a) <- compare16$Scenario[2:7]
round(compare16a, 3)
round(compare16a, 2)
}
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.