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