## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)
## ---- echo = FALSE, warning=F-------------------------------------------------
library(ggplot2, warn.conflicts = F, quietly = T)
library(raster, warn.conflicts = F, quietly = T)
## ----install.package, eval=FALSE----------------------------------------------
# # Install 'devtools' package, if needed
# if(!("devtools" %in% list.files(.libPaths()))) {
# install.packages("devtools", repo = "http://cran.rstudio.com", dep = TRUE)
# }
#
# # Download package, build vignette
# devtools::install_github("wpeterman/ResistanceGA",
# build_vignettes = TRUE)
## ----results='hide',message=FALSE, warning=FALSE------------------------------
library(ResistanceGA)
rm(list = ls())
## ----Plot.trans.demo, fig.height = 4, fig.width = 4---------------------------
Ricker.plot <- Plot.trans(PARM = c(1.5, 200),
Resistance = c(1,10),
transformation="Ricker")
# Change title of plot
Ricker.plot$labels$title <- "Ricker Transformation"
Ricker.plot
# Find original data value that now has maximum resistance
Ricker.plot$data$original[which(Ricker.plot$data$transformed==max(Ricker.plot$data$transformed))]
## ----warning=FALSE, results='hide',message=FALSE, eval=FALSE------------------
#
# if("ResistanceGA_Examples"%in%dir("C:/")==FALSE)
# dir.create(file.path("C:/", "ResistanceGA_Examples"))
#
# # Create a subdirectory for the first example
# dir.create(file.path("C:/ResistanceGA_Examples/","SingleSurface"))
#
# # Directory to write .asc files and results
# write.dir <- "C:/ResistanceGA_Examples/SingleSurface/"
#
# # Give path to CIRCUITSCAPE .exe file
# # Default = '"C:/Program Files/Circuitscape/cs_run.exe"'
# CS.program <- paste('"C:/Program Files/Circuitscape/cs_run.exe"')
## ----load.data, echo = FALSE, message = FALSE, warning=FALSE------------------
data(resistance_surfaces)
continuous <- resistance_surfaces[[2]]
data(samples)
sample.locales <- SpatialPoints(samples[,c(2,3)])
## ----eval=FALSE---------------------------------------------------------------
# data(resistance_surfaces)
# continuous <- resistance_surfaces[[2]]
# writeRaster(continuous,
# filename = paste0(write.dir,"cont.asc"),
# overwrite = TRUE)
## ----eval=FALSE---------------------------------------------------------------
# data(samples)
# write.table(samples,file=paste0(write.dir,"samples.txt"),sep="\t",col.names=F,row.names=F)
#
# # Create a spatial points object for plotting
# sample.locales <- SpatialPoints(samples[,c(2,3)])
## ----single.surface.plot, fig.height = 4, fig.width = 4-----------------------
plot(continuous)
plot(sample.locales, pch = 16, col = "blue", add = TRUE) # Add points
## ----eval=FALSE---------------------------------------------------------------
# # Set the random number seed to reproduce the results presented
# GA.inputs <- GA.prep(ASCII.dir = write.dir,
# max.cat = 500,
# max.cont = 500,
# select.trans = "M",
# method = "LL",
# seed = 555)
#
# CS.inputs <- CS.prep(n.Pops = length(sample.locales),
# CS_Point.File = paste0(write.dir,"samples.txt"),
# CS.program = CS.program)
## ----monomolec.plot, eval = FALSE---------------------------------------------
# r.tran <- Resistance.tran(transformation = "Monomolecular",
# shape = 2,
# max = 275,
# r = continuous)
#
# plot.t <- Plot.trans(PARM = c(2, 275),
# Resistance = continuous,
# transformation = "Monomolecular")
## ----eval=FALSE---------------------------------------------------------------
# # Create the true resistance/response surface
# CS.response <- Run_CS(CS.inputs = CS.inputs,
# GA.inputs = GA.inputs,
# r = r.tran)
## ----eval=FALSE---------------------------------------------------------------
# CS.inputs <- CS.prep(n.Pops = length(sample.locales),
# response = CS.response,
# CS_Point.File = paste0(write.dir,"samples.txt"),
# CS.program = CS.program)
## ----eval=FALSE---------------------------------------------------------------
# SS_RESULTS <- SS_optim(CS.inputs=CS.inputs,
# GA.inputs=GA.inputs)
## ----eval=FALSE---------------------------------------------------------------
# GA.inputs <- GA.prep(ASCII.dir = write.dir)
#
# CS.inputs <- CS.prep(n.Pops = length(sample.locales),
# response = CS.response,
# CS_Point.File = paste0(write.dir,"samples.txt"),
# CS.program = CS.program)
#
# SS_RESULTS <- SS_optim(CS.inputs = CS.inputs,
# GA.inputs = GA.inputs)
## ----gdistance, eval=FALSE----------------------------------------------------
# # Import data
# data(resistance_surfaces)
# continuous <- resistance_surfaces[[2]]
#
# data(samples)
# sample.locales <- SpatialPoints(samples[,c(2,3)])
#
# # Set the random number seed to reproduce the results presented
# # Run in parallel on 4 cores
# GA.inputs <- GA.prep(ASCII.dir=continuous,
# Results.dir=write.dir,
# select.trans = "M",
# max.cat=500,
# max.cont=500,
# seed = 555,
# parallel = 4)
#
#
# gdist.inputs <- gdist.prep(length(sample.locales),
# samples = sample.locales,
# response = CS.response,
# method = 'commuteDistance') ## Optimize using commute distance
#
# # Run optimization
# SS_RESULTS.gdist <- SS_optim(gdist.inputs = gdist.inputs,
# GA.inputs = GA.inputs)
## ----eval=FALSE---------------------------------------------------------------
# Grid.Results <- Grid.Search(shape = seq(1, 3, by = 0.025),
# max = seq(125, 425, by = 50),
# transformation = "Monomolecular",
# Resistance = continuous,
# gdist.inputs = gdist.inputs,
# GA.inputs = GA.inputs)
## ----recreate.grid, eval=FALSE------------------------------------------------
# filled.contour(Grid.Results$Plot.data,
# col = rainbow(14),
# xlab = "Shape parameter",
# ylab = "Maximum value parameter")
## ----warning=FALSE, results='hide',message=FALSE, eval=FALSE------------------
# if("ResistanceGA_Examples"%in%dir("C:/")==FALSE)
# dir.create(file.path("C:/", "ResistanceGA_Examples"))
#
# # Create a subdirectory for the second example
# dir.create(file.path("C:/ResistanceGA_Examples/","MultipleSurfaces"))
#
# # Directory to write .asc files and results
# write.dir <- "C:/ResistanceGA_Examples/MultipleSurfaces/"
## ----multi_surface.sim, warning=FALSE, message=FALSE,results='hide'-----------
data(resistance_surfaces)
data(samples)
sample.locales <- SpatialPoints(samples[ ,c(2,3)])
## ----feature.sim, warning=FALSE,message=FALSE, fig.height = 4, fig.width = 4----
plot(resistance_surfaces[[1]],main = resistance_surfaces[[1]]@data@names)
plot(sample.locales, pch=16, col="blue", add=TRUE)
plot(resistance_surfaces[[2]],main = resistance_surfaces[[2]]@data@names)
plot(sample.locales, pch=16, col="blue", add=TRUE)
plot(resistance_surfaces[[3]],main = resistance_surfaces[[3]]@data@names)
plot(sample.locales, pch=16, col="blue", add=TRUE)
## ----eval=FALSE---------------------------------------------------------------
# # Note that the `resistance_surfaces` is already a RasterStack object.
# # The code below for demonstration of how to make a stack.
# r.stack <- stack(resistance_surfaces$categorical,
# resistance_surfaces$continuous,
# resistance_surfaces$feature)
#
# GA.inputs <- GA.prep(ASCII.dir = r.stack,
# Results.dir = write.dir,
# method = "LL",
# max.cat = 500,
# max.cont = 500,
# seed = 555,
# parallel = 4)
#
# gdist.inputs <- gdist.prep(length(sample.locales),
# samples = sample.locales,
# method = 'commuteDistance') # Optimize using commute distance
#
## ----IR_Mono, eval=FALSE------------------------------------------------------
# plot.t <- Plot.trans(PARM = c(3.5, 150),
# Resistance = continuous,
# transformation = "Inverse-Reverse Monomolecular")
## ----combine.surfaces, eval=FALSE---------------------------------------------
# PARM <- c(1, 250, 75, 1, 3.5, 150, 1, 350)
#
# # PARM<- c(1, # First feature of categorical
# # 250, # Second feature of categorical
# # 75, # Third feature of categorical
# # 1, # Transformation equation for continuous surface = Inverse-Reverse Monomolecular
# # 3.5, # Shape parameter
# # 150, # Scale parameter
# # 1, # First feature of feature surface
# # 350) # Second feature of feature surface
#
# # Combine resistance surfaces
# Resist <- Combine_Surfaces(PARM = PARM,
# gdist.inputs = gdist.inputs,
# GA.inputs = GA.inputs,
# out = NULL,
# rescale = TRUE)
#
# # View combined surface
# plot(Resist, main = "scaled composite resistance")
## ----combine.cs, eval=FALSE---------------------------------------------------
# # Create the true resistance/response surface
# gdist.response <- Run_gdistance(gdist.inputs = gdist.inputs,
# r = Resist)
#
# gdist.inputs <- gdist.prep(n.Pops = length(sample.locales),
# samples = sample.locales,
# response = as.vector(gdist.response),
# method = 'commuteDistance')
## ----eval=FALSE---------------------------------------------------------------
# Multi.Surface_optim <- MS_optim(gdist.inputs = gdist.inputs,
# GA.inputs = GA.inputs)
## ----eval=FALSE---------------------------------------------------------------
# Summary.table <- data.frame(PARM,round(t(Multi.Surface_optim$GA.summary@solution),2))
# colnames(Summary.table)<-c("Truth", "Optimized")
# row.names(Summary.table)<-c("Category1", "Category2", "Category3",
# "Transformation", "Shape", "Max",
# "Feature1", "Feature2")
## ----combined.plots,fig.width=12,fig.height=8, eval=FALSE---------------------
# # Make combined, optimized resistance surface.
# optim.resist <- Combine_Surfaces(PARM = Multi.Surface_optim$GA.summary@solution,
# gdist.inputs = gdist.inputs,
# GA.inputs = GA.inputs,
# rescale = TRUE)
# ms.stack <- stack(Resist, optim.resist)
# names(ms.stack) <- c("Truth", "Optimized")
# plot(ms.stack)
#
# # Correlation between the two surfaces
# pairs(ms.stack)
## ----CS.maps, eval=FALSE------------------------------------------------------
# # Note: You must run `CS.prep` to generate the CS.inputs object for doing this.
# CS.inputs <- CS.prep(n.Pops = length(sample.locales),
# response = gdist.response,
# CS_Point.File = paste0(write.dir,"samples.txt"),
# CS.program = CS.program)
#
# Resist.true <- Run_CS(CS.inputs = CS.inputs,
# GA.inputs = GA.inputs,
# r = Resist,
# CurrentMap = TRUE,
# output = "raster")
#
# Resist.opt <- Run_CS(CS.inputs = CS.inputs,
# GA.inputs = GA.inputs,
# r = optim.resist,
# CurrentMap = TRUE,
# output = "raster")
#
# # We can confirm that, like the resistance surfaces above,
# # the CIRCUITSCAPE current maps are also correlated
# cs.stack <- stack(Resist.true, Resist.opt)
# names(cs.stack) <- c("Truth", "Optimized")
# pairs(cs.stack)
## ----eval = FALSE-------------------------------------------------------------
# data(resistance_surfaces)
# cat <- resistance_surfaces[[1]]
# cat[cat < 2] <- 0
#
# # Make categorical surface binary
# cat[cat == 2] <- 1
#
# # Smooth and visualize
# # The `SCALE` parameter re-scales the surface to 0-10
# cat.smooth <- k.smooth(raster = cat,
# sigma = 1,
# SCALE = TRUE)
# par(mfrow = c(1,2))
# plot(cat, main = "Original")
# plot(cat.smooth, main = "Smoothed, sigma = 1")
# par(mfrow = c(1,1))
## ----k_smooth, eval=FALSE-----------------------------------------------------
# data(samples)
# sample.locales <- SpatialPoints(samples[,c(2,3)])
#
# # Set the random number seed to reproduce the results presented
# # Run in parallel on 4 cores
# # NOTE: `scale = TRUE` to indicate optimization of scaling/smoothing parameter
# GA.inputs <- GA.prep(ASCII.dir = cat,
# Results.dir = write.dir,
# select.trans = "M",
# scale = TRUE,
# max.cat = 500,
# max.cont = 500,
# seed = 321,
# run = 35,
# parallel = 4)
#
# # Optimize using commute distance
# gdist.inputs <- gdist.prep(n.Pops = length(sample.locales),
# samples = sample.locales,
# method = 'commuteDistance')
#
# # Transform resistance surface
# r.tran_smooth <- Resistance.tran(transformation = "Monomolecular",
# shape = 2,
# max = 275,
# r = cat.smooth)
#
# # Create the true resistance/response surface
# gdist.response <- Run_gdistance(gdist.inputs = gdist.inputs,
# r = r.tran_smooth)
#
# # Rerun `gdist.prep` to include response
# gdist.inputs <- gdist.prep(n.Pops = length(sample.locales),
# response = as.vector(gdist.response),
# samples = sample.locales,
# method = 'commuteDistance')
#
# # Run optimization: NOTE use of `SS_optim.scale` to optimize the smoothing parameter
# SS_RESULTS.gdist_scale <- SS_optim.scale(gdist.inputs = gdist.inputs,
# GA.inputs = GA.inputs)
## ----analysis1, eval=FALSE----------------------------------------------------
# data(samples)
# data("resistance_surfaces")
#
# # Create a spatial points object
# sample.locales <- SpatialPoints(samples[, c(2, 3)])
#
# # Run `gdist.prep` & GA.prep
# gdist.inputs <- gdist.prep(n.Pops = length(sample.locales),
# samples = sample.locales,
# method = 'commuteDistance')
#
# # This will be used again later
# # Note: to speed up the analysis, only Monomolecular tranformations will be assessed
# GA.inputs_NoFeature <- GA.prep(method = "LL",
# ASCII.dir = resistance_surfaces[[-3]],
# Results.dir = "C:/ResistanceGA_Examples/run2/",
# max.cat = 500,
# max.cont = 500,
# select.trans = list(NA,
# "M"),
# seed = 123,
# parallel = 6)
#
# # The 'true' resistance surface will be the composite surface
# # Combine resistance surfaces, omitting the feature surface
# # Use an Inverse-Reverse Monomolecular transformation of the continuous surface
# # Inverse-Reverse Monomolecular = 1
# PARM <- c(1, 250, 100, 1, 1.5, 150)
#
# # Setting `p.contribution = TRUE` to see how each surface
# # contributes to the total resistance of the composite surface
# # This is the 'true' resistance surface that the example 'response'
# # data were generated from
# Resist <- Combine_Surfaces(PARM = PARM,
# gdist.inputs = gdist.inputs,
# GA.inputs = GA.inputs_NoFeature,
# out = NULL,
# rescale = TRUE,
# p.contribution = TRUE)
#
# # Assess contribution of each surface
# Resist$percent.contribution
## ----analysis2, eval=FALSE----------------------------------------------------
# # Create a subdirectory for results
# dir.create(file.path("C:/ResistanceGA_Examples/","run1"))
# dir.create(file.path("C:/ResistanceGA_Examples/","run2"))
#
# # Turn response data into vector
# gd.true <- Run_gdistance(gdist.inputs = gdist.inputs,
# r = Resist$combined.surface)
#
# gd.true <- as.vector(gd.true)
#
# # Add some noise to response
# set.seed(321)
# gd.response <- gd.true + rnorm(length(gd.true), 0, 9)
#
# plot(gd.response ~ gd.true)
# ecodist::mantel(gd.response ~ gd.true) # Mantel r = 0.67
#
# # Correlation with distance
# ecodist::mantel(gd.response ~ as.vector(dist(samples[, c(2, 3)]))) # Mantel r = 0.26
# plot(gd.response ~ as.vector(dist(samples[, c(2, 3)])), xlab = "Euclidean distance")
## ----analysis2b, eval=FALSE---------------------------------------------------
# # Re-run `gdist.prep` function
# gdist.inputs <- gdist.prep(n.Pops = length(sample.locales),
# response = gd.response,
# samples = sample.locales)
#
#
# # Re-run GA.prep to include all surfaces
# # Note: to speed up the analysis, only Monomolecular tranformations will be assessed
#
# GA.inputs_All <- GA.prep(method = "LL",
# ASCII.dir = resistance_surfaces,
# Results.dir = "C:/ResistanceGA_Examples/run1/",
# max.cat = 500,
# max.cont = 500,
# select.trans = list(NA,
# "M",
# NA),
# seed = 123,
# parallel = 6)
#
# # First run all single surfaces, Multi-surface is response variable
# SS_RESULTS.gdist <- SS_optim(gdist.inputs = gdist.inputs,
# GA.inputs = GA.inputs_All)
#
# # Run `MS_optim` with all surfaces
# Multi.Surface_optim.gd <- MS_optim(gdist.inputs = gdist.inputs,
# GA.inputs = GA.inputs_All)
#
# # Run `MS_optim` with without Feature surface
# Multi.Surface_optim.gd2 <- MS_optim(gdist.inputs = gdist.inputs,
# GA.inputs = GA.inputs_NoFeature)
## ----analysis3, eval=FALSE----------------------------------------------------
# Multi.Surface_optim.gd$percent.contribution
#
# Multi.Surface_optim.gd2$percent.contribution
## ----boot, eval = F-----------------------------------------------------------
# # Extract relevant components from optimization outputs
# # Make a list of cost/resistance distance matrices
# mat.list <- c(SS_RESULTS.gdist$cd,
# Multi.Surface_optim.gd$cd,
# Multi.Surface_optim.gd2$cd)
#
# k <- rbind(SS_RESULTS.gdist$k,
# Multi.Surface_optim.gd$k,
# Multi.Surface_optim.gd2$k)
#
# # Create square distance matrix for response for use with
# # the bootstrap function
#
# response <- matrix(0, 25, 25)
# response[lower.tri(response)] <- gd.response
#
# # Run bootstrap
# (AIC.boot <- Resist.boot(mod.names = names(mat.list),
# dist.mat = mat.list,
# n.parameters = k[,2],
# sample.prop = 0.75,
# iters = 1000,
# obs = 25,
# genetic.mat = response
# )
# )
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.